Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 18487)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 18488)
@@ -2930,4 +2930,7 @@
 				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,14 +3149,24 @@
 			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,5 +3177,5 @@
 			Bh.CreateSingleVertexToTriangleConnectivity();     
 			InsertNewPoints(nbvold,NbTSwap,bamgopts->random);
-		}  
+		}
 		else Bh.CreateSingleVertexToTriangleConnectivity();