Changeset 18488

Timestamp:

09/09/14 05:40:32 (11 years ago)

Author:

olason

Message:

CHG: do not include any vertex of the background mesh that may lie outside of the domain

File:

• issm/trunk-jpl/src/c/bamg/Mesh.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/bamg/Mesh.cpp

@@ -2930 +2930 @@
                                 Triangle *tcvj=TriangleFindFromCoord(vj.i,det3);
                                 if (tcvj && !tcvj->link){
+                                        _printf_("While trying to add the following point:\n");
+                                        vj.Echo();
+                                        _printf_("BAMG determined that it was inside the following triangle, which probably lies outside of the geometric domain\n");
                                         tcvj->Echo();
                                         _error_("problem inserting point in InsertNewPoints (tcvj=" << tcvj << " and tcvj->link=" << tcvj->link << ")");
@@ -3146 +3149 @@
                         if (verbose>5) _printf_("         Inserting initial mesh points\n");
                         bool pointsoutside = false;
-                        for (i=0;i<Bh.nbv;i++){ 
+                        for(i=0;i<Bh.nbv;i++){ 
                                 BamgVertex &bv=Bh[i];
                                 /*Do not insert if the point is outside*/
                                 long long det3[3];
                                 Triangle* tcvj=TriangleFindFromCoord(bv.i,det3);
-                                if(tcvj->det<0){
+                                if(tcvj->det<0 || !tcvj->link){
                                         pointsoutside = true;
                                         continue;
                                 }
-                                if (!bv.GeomEdgeHook){
+                                IssmDouble area_1=((bv.r.x -(*tcvj)(2)->r.x)*((*tcvj)(1)->r.y-(*tcvj)(2)->r.y) 
+                                                - (bv.r.y -(*tcvj)(2)->r.y)*((*tcvj)(1)->r.x-(*tcvj)(2)->r.x));
+                                IssmDouble area_2=(((*tcvj)(0)->r.x -(*tcvj)(2)->r.x)*(bv.r.y -(*tcvj)(2)->r.y) 
+                                                - ((*tcvj)(0)->r.y -(*tcvj)(2)->r.y)*(bv.r.x -(*tcvj)(2)->r.x));
+                                IssmDouble area_3 =((bv.r.x -(*tcvj)(1)->r.x)*((*tcvj)(0)->r.y-(*tcvj)(1)->r.y)
+                                                - (bv.r.y -(*tcvj)(1)->r.y)*((*tcvj)(0)->r.x-(*tcvj)(1)->r.x));
+                                if(area_1<0 || area_2<0 || area_3<0){
+                                        pointsoutside = true;
+                                        continue;
+                                }
+                                if(!bv.GeomEdgeHook){
                                         vertices[nbv].r              = bv.r;
                                         vertices[nbv].PreviousNumber = i+1;
@@ -3164 +3177 @@
                         Bh.CreateSingleVertexToTriangleConnectivity();     
                         InsertNewPoints(nbvold,NbTSwap,bamgopts->random);
-                }  
+                }
                 else Bh.CreateSingleVertexToTriangleConnectivity();