Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 22497)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 22498)
@@ -561,11 +561,11 @@
 			./shared/Threads/PartitionRange.cpp\
 			./shared/Exp/exp.cpp\
-			./shared/TriMesh/AssociateSegmentToElement.cpp\
-			./shared/TriMesh/GridInsideHole.cpp\
-			./shared/TriMesh/OrderSegments.cpp\
-			./shared/TriMesh/SplitMeshForRifts.cpp\
-			./shared/TriMesh/TriMeshUtils.cpp\
-			./modules/TriMeshx/TriMeshx.cpp\
-			./modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp\
+			./shared/Triangle/AssociateSegmentToElement.cpp\
+			./shared/Triangle/GridInsideHole.cpp\
+			./shared/Triangle/OrderSegments.cpp\
+			./shared/Triangle/SplitMeshForRifts.cpp\
+			./shared/Triangle/TriangleUtils.cpp\
+			./modules/Trianglex/Trianglex.cpp\
+			./modules/ProcessRiftsx/ProcessRiftsx.cpp\
 			./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp\
 			./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp\
Index: /issm/trunk-jpl/src/c/main/issm.js
===================================================================
--- /issm/trunk-jpl/src/c/main/issm.js	(revision 22497)
+++ /issm/trunk-jpl/src/c/main/issm.js	(revision 22498)
@@ -11,5 +11,5 @@
 	binHeap.set(new Uint8Array(dbinary.buffer)); var binary=binHeap.byteOffset;
 
-	//Declare TriMesh module: 
+	//Declare module: 
 	issm= Module.cwrap('main','number',['number','number']);
 	
Index: /issm/trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp	(revision 22497)
+++ /issm/trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp	(revision 22498)
@@ -1,4 +1,4 @@
-/*!\file TriMeshx
- * \brief: x code for TriMesh mesher
+/*!\file CoordinateSystemTransformx
+ * \brief: x code for CoordinateSystemTransformx
  */
 
Index: /issm/trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp	(revision 22498)
@@ -0,0 +1,78 @@
+/*!\file:  ProcessRifts.cpp
+ * \brief split a mesh where a rift (or fault) is present
+ */ 
+
+#include "./ProcessRiftsx.h"
+#include "../../classes/RiftStruct.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+
+void ProcessRiftsx(int** pindex, int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct){
+
+	/*Output*/
+	int      numrifts,numrifts0;
+	int     *riftsnumsegments     = NULL;
+	int    **riftssegments        = NULL;
+	int     *riftsnumpairs        = NULL;
+	int    **riftspairs           = NULL;
+	int     *riftstips            = NULL;
+	double **riftspenaltypairs    = NULL;
+	int     *riftsnumpenaltypairs = NULL;
+
+	/*Recover initial mesh*/
+	int     nel            = *pnel;
+	int    *index          = *pindex;
+	double *x              = *px;
+	double *y              = *py;
+	int     nods           = *pnods;
+	int    *segments       = *psegments;
+	int    *segmentmarkers = *psegmentmarkers;
+	int     num_seg        = *pnum_seg;
+
+	/*Intermediary*/
+	int     riftflag;
+
+	/*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie: 
+	 *all the nodes of this element belong to the segments (tends to happen when there are corners: */
+	RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
+
+	/*Figure out if we have rifts, and how many: */
+	IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
+
+	if(!riftflag) _error_("No rift present in mesh");
+
+	/*Split mesh*/
+	SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
+
+	/*Order segments so that their normals point outside the domain: */
+	OrderSegments(&segments,num_seg, index,nel);
+
+	/*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the 
+	 *segmentmarkerlist:*/
+	SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
+
+	/*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
+	PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
+
+	/*Order rifts so that they start from one tip, go to the other tip, and back: */
+	OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
+
+	/*Create penalty pairs, used by Imp: */
+	PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
+
+	/*Create Riftstruct*/
+	RiftStruct* riftstruct = new RiftStruct(numrifts,riftsnumsegments,riftssegments,riftsnumpairs,riftspairs,riftsnumpenaltypairs,riftspenaltypairs,riftstips);
+
+	/*Assign output pointers for mesh*/
+	*pnel            = nel;
+	*pindex          = index;
+	*px              = x;
+	*py              = y;
+	*pnods           = nods;
+	*psegments       = segments;
+	*psegmentmarkers = segmentmarkers;
+	*pnum_seg        = num_seg;
+
+	/*Assign output pointers for rifts*/
+	*priftstruct = riftstruct;
+}
Index: /issm/trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h	(revision 22498)
+++ /issm/trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h	(revision 22498)
@@ -0,0 +1,12 @@
+/*!\file:  ProcessRiftsx.h
+ * \brief header file for ProcessRifts module
+ */ 
+
+#ifndef _PROCESSRIFTX_H
+#define _PROCESSRIFTX_H
+
+class RiftStruct;
+
+void ProcessRiftsx(int** pindex,int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct);
+
+#endif  /* _PROCESSRIFTX_H*/
Index: /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp	(revision 22498)
@@ -0,0 +1,201 @@
+/*!\file Trianglex
+ * \brief: x code for Triangle mesher
+ */
+
+/*Header files: {{{*/
+#include "./Trianglex.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+/*ANSI_DECLARATORS needed to call triangle library: */
+#if defined(_HAVE_TRIANGLE_)
+	#ifndef ANSI_DECLARATORS
+	#define ANSI_DECLARATORS
+	#endif
+	#include "triangle.h"
+	#undef ANSI_DECLARATORS
+#endif
+/*}}}*/
+
+void Trianglex(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnsegs,Contours* domain,Contours* rifts,double area){
+
+#if !defined(_HAVE_TRIANGLE_)
+	_error_("triangle has not been installed");
+#else
+	/*indexing: */
+	int i,j;
+
+	/*output: */
+	int    *index             = NULL;
+	double *x                 = NULL;
+	double *y                 = NULL;
+	int    *segments          = NULL;
+	int    *segmentmarkerlist = NULL;
+
+	/*intermediary: */
+	int counter,counter2,backcounter;
+	Contour<IssmPDouble> *contour = NULL;
+
+	/* Triangle structures needed to call Triangle library routines: */
+	struct triangulateio in,out;
+	char   options[256];
+
+	/*Create initial triangulation to call triangulate(). First number of points:*/
+	in.numberofpoints=0;
+	for (i=0;i<domain->Size();i++){
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+		in.numberofpoints+=contour->nods-1;
+	}
+	for (i=0;i<rifts->Size();i++){
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+		in.numberofpoints+=contour->nods;
+	}
+
+	/*number of point attributes: */
+	in.numberofpointattributes=1;
+
+	/*fill in the point list: */
+	in.pointlist = xNew<REAL>(in.numberofpoints*2);
+
+	counter=0;
+	for (i=0;i<domain->Size();i++){
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+		for (j=0;j<contour->nods-1;j++){
+			in.pointlist[2*counter+0]=contour->x[j];
+			in.pointlist[2*counter+1]=contour->y[j];
+			counter++;
+		}
+	}
+	for (i=0;i<rifts->Size();i++){
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+		for (j=0;j<contour->nods;j++){
+			in.pointlist[2*counter+0]=contour->x[j];
+			in.pointlist[2*counter+1]=contour->y[j];
+			counter++;
+		}
+	}
+
+	/*fill in the point attribute list: */
+	in.pointattributelist = xNew<REAL>(in.numberofpoints*in.numberofpointattributes);
+	for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0;
+
+	/*fill in the point marker list: */
+	in.pointmarkerlist = xNew<int>(in.numberofpoints);
+	for(i=0;i<in.numberofpoints;i++) in.pointmarkerlist[i] = 0;
+
+	/*Build segments. First figure out number of segments: holes and closed outlines have as many segments as vertices: */
+	in.numberofsegments=0;
+	for (i=0;i<domain->Size();i++){
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+		in.numberofsegments+=contour->nods-1;
+	}
+	for(i=0;i<rifts->Size();i++){
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+		/*for rifts, we have one less segment as we have vertices*/
+		in.numberofsegments+=contour->nods-1;
+	}
+
+	in.segmentlist = xNew<int>(in.numberofsegments*2);
+	in.segmentmarkerlist = xNewZeroInit<int>(in.numberofsegments);
+	counter=0;
+	backcounter=0;
+	for (i=0;i<domain->Size();i++){
+		contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
+		for (j=0;j<contour->nods-2;j++){
+			in.segmentlist[2*counter+0]=counter;
+			in.segmentlist[2*counter+1]=counter+1;
+			in.segmentmarkerlist[counter]=0;
+			counter++;
+		}
+		/*Close this profile: */
+		 in.segmentlist[2*counter+0]=counter;
+		 in.segmentlist[2*counter+1]=backcounter;
+		 in.segmentmarkerlist[counter]=0;
+		 counter++;
+		 backcounter=counter;
+	}
+	counter2=counter;
+	for (i=0;i<rifts->Size();i++){
+		contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
+		for (j=0;j<(contour->nods-1);j++){
+			in.segmentlist[2*counter2+0]=counter;
+			in.segmentlist[2*counter2+1]=counter+1;
+			in.segmentmarkerlist[counter2]=2+i;
+			counter2++;
+			counter++;
+		}
+		counter++;
+	}
+
+	/*Build regions: */
+	in.numberofregions = 0;
+
+	/*Build holes: */
+	in.numberofholes = domain->Size()-1; /*everything is a hole, but for the first profile.*/
+	if(in.numberofholes){
+		in.holelist = xNew<REAL>(in.numberofholes*2);
+		for (i=0;i<domain->Size()-1;i++){
+			contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
+			GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
+		}
+	}
+
+	/* Make necessary initializations so that Triangle can return a triangulation in `out': */
+	out.pointlist             = (REAL*)NULL;
+	out.pointattributelist    = (REAL*)NULL;
+	out.pointmarkerlist       = (int *)NULL;
+	out.trianglelist          = (int *)NULL;
+	out.triangleattributelist = (REAL*)NULL;
+	out.neighborlist          = (int *)NULL;
+	out.segmentlist           = (int *)NULL;
+	out.segmentmarkerlist     = (int *)NULL;
+	out.edgelist              = (int *)NULL;
+	out.edgemarkerlist        = (int *)NULL;
+
+	/* Triangulate the points:.  Switches are chosen to read and write a  */
+	/*   PSLG (p), preserve the convex hull (c), number everything from  */
+	/*   zero (z), assign a regional attribute to each element (A), and  */
+	/*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
+	/*   neighbor list (n).                                              */
+	sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
+	triangulate(options, &in, &out, NULL);
+	/*report(&out, 0, 1, 1, 1, 1, 0);*/
+
+	/*Allocate index, x and y: */
+	index=xNew<int>(3*out.numberoftriangles);
+	x=xNew<double>(out.numberofpoints);
+	y=xNew<double>(out.numberofpoints);
+	segments=xNew<int>(3*out.numberofsegments);
+	segmentmarkerlist=xNew<int>(out.numberofsegments);
+
+	for (i = 0; i< out.numberoftriangles; i++) {
+		for (j = 0; j < out.numberofcorners; j++) {
+			index[3*i+j]=(int)out.trianglelist[i*out.numberofcorners+j]+1;
+		}
+	}
+	for (i = 0; i< out.numberofpoints; i++){
+		x[i]=(double)out.pointlist[i*2+0];
+		y[i]=(double)out.pointlist[i*2+1];
+	}
+	for (i = 0; i<out.numberofsegments;i++){
+		segments[3*i+0]=(int)out.segmentlist[i*2+0]+1;
+		segments[3*i+1]=(int)out.segmentlist[i*2+1]+1;
+		segmentmarkerlist[i]=(int)out.segmentmarkerlist[i];
+	}
+
+	/*Associate elements with segments: */
+	AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
+
+	/*Order segments so that their normals point outside the domain: */
+	OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
+
+	/*Output : */
+	*pindex=index;
+	*px=x;
+	*py=y;
+	*psegments=segments;
+	*psegmentmarkerlist=segmentmarkerlist;
+	*pnels=out.numberoftriangles;
+	*pnods=out.numberofpoints;
+	*pnsegs=out.numberofsegments;
+#endif
+}
Index: /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.h	(revision 22498)
+++ /issm/trunk-jpl/src/c/modules/Trianglex/Trianglex.h	(revision 22498)
@@ -0,0 +1,13 @@
+/*!\file:  Trianglex.h
+ * \brief header file for Trianglex module
+ */ 
+
+#ifndef _TRIANGLEX_H_
+#define _TRIANGLEX_H_
+
+#include <string.h>
+#include "../../classes/classes.h"
+
+/* local prototypes: */
+void Trianglex(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnseg,Contours* domain,Contours* rifts,double area);
+#endif  /* _TRIANGLEX_H */
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 22497)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 22498)
@@ -91,6 +91,6 @@
 #include "./SpcNodesx/SpcNodesx.h"
 #include "./SurfaceAreax/SurfaceAreax.h"
-#include "./TriMeshx/TriMeshx.h"
-#include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h"
+#include "./Trianglex/Trianglex.h"
+#include "./ProcessRiftsx/ProcessRiftsx.h"
 #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
 #include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h"
Index: /issm/trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp	(revision 22498)
@@ -0,0 +1,24 @@
+/*!\file:  AssociateSegmentToElement.cpp
+ * \brief for each segment, look for the corresponding element.
+ */ 
+
+#include "./trimesh.h"
+
+int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel){
+
+	/*node indices: */
+	int A,B;
+
+	/*Recover segments: */
+	int* segments=*psegments;
+
+	for(int i=0;i<nseg;i++){
+		A=segments[3*i+0];
+		B=segments[3*i+1];
+		segments[3*i+2]=FindElement(A,B,index,nel)+1; //matlab indexing.
+	}
+
+	/*Assign output pointers: */
+	*psegments=segments;
+	return 1;
+}
Index: /issm/trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp	(revision 22498)
@@ -0,0 +1,51 @@
+/*
+ * GridInsideHole.c:
+ * from a convex set of points, figure out a point that for sure lies inside the profile.
+ */
+
+#include <math.h>
+
+#include "./trimesh.h"
+#include "../Exp/exp.h"
+
+#undef M_PI
+#define M_PI 3.141592653589793238462643
+
+int GridInsideHole(double* px0,double* py0,int n,double* x,double* y){
+
+	double flag=0.0;
+	double xA,xB,xC,xD,xE;
+	double yA,yB,yC,yD,yE;
+
+	/*Take first and last vertices: */
+	xA=x[0];
+	yA=y[0];
+	xB=x[n-1];
+	yB=y[n-1];
+
+	/*Figure out middle of segment [A B]: */
+	xC=(xA+xB)/2;
+	yC=(yA+yB)/2;
+
+	/*D and E are on each side of segment [A B], on the median line between segment [A  B], 
+	 *at an angle of 10 degree (less than the minimum 30 enforced by the quality of the mesh: */
+	xD=xC+tan(10./180.*M_PI)*(yC-yA);
+	yD=yC+tan(10./180.*M_PI)*(xA-xC);
+	xE=xC-tan(10./180.*M_PI)*(yC-yA);
+	yE=yC-tan(10./180.*M_PI)*(xA-xC);
+
+	/*Either E or D is inside profile (x,y): */
+	IsInPolySerial(&flag,&xD,&yD,1,x,y,n,2);
+	/*FIXME: used to be 'flag' and not '!flag', check*/
+	if(!flag){
+		/*D is inside the poly: */
+		*px0=xD;
+		*py0=yD;
+	}
+	else{
+		/*E is inside the poly: */
+		*px0=xE;
+		*py0=yE;
+	}
+	return 1;
+}
Index: /issm/trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp	(revision 22498)
@@ -0,0 +1,46 @@
+/*
+ * OrderSegments.c: 
+ * reorder segments so that their normals point outside the domain outline.
+ */
+#include "./trimesh.h"
+
+int OrderSegments(int** psegments,int nseg,int* index,int nel){
+
+	/*vertex indices: */
+	int A,B;
+
+	/*element index*/
+	int el;
+
+	/*Recover segments: */
+	int* segments=*psegments;
+
+	for(int i=0;i<nseg;i++){
+		A=segments[3*i+0];
+		B=segments[3*i+1];
+		el=segments[3*i+2]-1; //after AssociateSegmentToElement, el was a matlab index, we need the c index now.
+
+		if (index[3*el+0]==A){
+			if (index[3*el+2]==B){
+				segments[3*i+0]=B;
+				segments[3*i+1]=A;
+			}
+		}
+		else if (index[3*el+1]==A){
+			if (index[3*el+0]==B){
+				segments[3*i+0]=B;
+				segments[3*i+1]=A;
+			}
+		}
+		else{
+			if (index[3*el+1]==B){
+				segments[3*i+0]=B;
+				segments[3*i+1]=A;
+			}
+		}
+	}
+
+	/*Assign output pointers: */
+	*psegments=segments;
+	return 1;
+}
Index: /issm/trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp	(revision 22498)
@@ -0,0 +1,96 @@
+/*
+ * SplitMeshForRifts.c:
+ */
+#include "./trimesh.h"
+#include "../MemOps/MemOps.h"
+
+int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist){
+
+	/*Some notes on dimensions: 
+	index  of size nelx3
+	x and y of size nodsx1
+	segments of size nsegsx3*/
+
+	int i,j,k,l;
+	int node;
+	int el;
+	int  nriftsegs;
+	int* riftsegments=NULL; 
+	int* flags=NULL;
+	int  NumGridElementListOnOneSideOfRift;
+	int* GridElementListOnOneSideOfRift=NULL;
+
+	/*Recover input: */
+	int     nel               = *pnel;
+	int    *index             = *pindex;
+	int     nods              = *pnods;
+	double *x                 = *px;
+	double *y                 = *py;
+	int     nsegs             = *pnsegs;
+	int    *segments          = *psegments;
+	int    *segmentmarkerlist = *psegmentmarkerlist;
+
+	/*Establish list of segments that belong to a rift: */
+	/*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
+	RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);
+
+	/*Go through all nodes of the rift segments, and start splitting the mesh: */
+	flags=xNewZeroInit<int>(nods); //to make sure we don't split the same nodes twice!
+	for (i=0;i<nriftsegs;i++){
+		for (j=0;j<2;j++){
+
+			node=riftsegments[4*i+j+2];
+			if(flags[node-1]){
+				/*This node was already split, skip:*/
+				continue;
+			}
+			else{
+				flags[node-1]=1;
+			}
+
+			if(IsGridOnRift(riftsegments,nriftsegs,node)){
+
+				DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
+
+				/*Summary: we have for node, a list of elements
+				 * (GridElementListOnOneSideOfRift, of size
+				 * NumGridElementListOnOneSideOfRift) that all contain node 
+				 *and that are on the same side of the rift. For all these
+				 elements, we clone node into another node, and we swap all
+				 instances of node in the triangulation *for those elements, to the
+				 new node.*/
+
+				//create new node
+				x=xReNew<double>(x,nods,nods+1);
+				y=xReNew<double>(y,nods,nods+1);
+				x[nods]=x[node-1]; //matlab indexing
+				y[nods]=y[node-1]; //matlab indexing
+
+				//augment number of nodes 
+				nods++;
+
+				//change elements owning this node
+				for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
+					el=GridElementListOnOneSideOfRift[k];
+					for (l=0;l<3;l++){
+						if (index[3*el+l]==node) index[3*el+l]=nods; //again, matlab indexing.
+					}
+				}
+			}// if(IsGridOnRift(riftsegments,nriftsegs,node))
+		} //for(j=0;j<2;j++)
+	} //for (i=0;i<nriftsegs;i++)
+
+	/*update segments: they got modified completely by adding new nodes.*/
+	UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
+
+	/*Assign output pointers: */
+	*pnel=nel;
+	*pindex=index;
+	*pnods=nods;
+	*px=x;
+	*py=y;
+	*pnsegs=nsegs;
+	*psegments=segments;
+	*psegmentmarkerlist=segmentmarkerlist;
+	return 1;
+}
Index: /issm/trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp	(revision 22498)
@@ -0,0 +1,912 @@
+/*
+ * TriangleUtils: mesh manipulation routines: 
+ */
+
+#include <stdio.h>
+
+#include "./triangle.h"
+#include "../Exceptions/exceptions.h"
+#include "../MemOps/MemOps.h"
+
+#define RIFTPENALTYPAIRSWIDTH 8
+int IsGridOnRift(int* riftsegments, int nriftsegs, int node){/*{{{*/
+
+	/*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift, 
+	 *if it belongs to 2 elements, it is on the tip of a rift, or it has already been split across the rift (see below).*/
+
+	int i;
+	int j;
+	int count;
+
+	count=0;
+	for (i=0;i<nriftsegs;i++){
+		for (j=0;j<2;j++){
+			if ((*(riftsegments+4*i+2+j))==node) count++;
+		}
+	}
+	if (count==2){
+		return 1;
+	}
+	else{
+		return 0;
+	}
+}/*}}}*/
+int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){/*{{{*/
+
+	/*From a node, recover all the elements that are connected to it: */
+	int i,j;
+	int noerr=1;
+
+	int max_number_elements=12;
+	int current_size;
+	int NumGridElements;
+	int* GridElements=NULL;
+	int* GridElementsRealloc=NULL;
+
+	/*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own 
+	 * the node. We start by allocating GridElements with that size, and realloc 
+	 * more if needed.*/
+
+	current_size=max_number_elements;
+	NumGridElements=0;
+	GridElements=xNew<int>(max_number_elements);
+
+	for (i=0;i<nel;i++){
+		for (j=0;j<3;j++){
+			if (index[3*i+j]==node){
+				if (NumGridElements<=(current_size-1)){
+					GridElements[NumGridElements]=i;
+					NumGridElements++;
+					break;
+				}
+				else{
+					/*Reallocate another max_number_elements slots in the GridElements: */
+					GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
+					if (!GridElementsRealloc){
+						noerr=0;
+						goto cleanup_and_return;
+					}
+					current_size+=max_number_elements;
+					GridElements=GridElementsRealloc;
+					GridElements[NumGridElements]=i;
+					NumGridElements++;
+					break;
+				}
+			}
+		}
+	}
+	cleanup_and_return:
+	if(!noerr){
+		xDelete<int>(GridElements);
+	}
+	/*Allocate return pointers: */
+	*pGridElements=GridElements;
+	*pNumGridElements=NumGridElements;
+	return noerr;
+}/*}}}*/
+int IsNeighbor(int el1,int el2,int* index){/*{{{*/
+	/*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
+	int i,j;
+	int count=0;
+	for (i=0;i<3;i++){
+		for (j=0;j<3;j++){
+			if (index[3*el1+i]==index[3*el2+j])count++;
+		}
+	}
+	if (count==2){
+		return 1;
+	}
+	else{
+		return 0;
+	}
+}/*}}}*/
+int IsOnRift(int el,int nriftsegs,int* riftsegments){/*{{{*/
+	/*From a list of elements segments, figure out if el belongs to it: */
+	int i;
+	for (i=0;i<nriftsegs;i++){
+		if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
+			return 1;
+		}
+	}
+	return 0;
+}/*}}}*/
+void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){/*{{{*/
+
+	int i,counter;
+	int el,el2;
+
+	int  nriftsegs;
+	int* riftsegments=NULL;
+	int* riftsegments_uncompressed=NULL; 
+	int  element_nodes[3];
+
+	/*Allocate segmentflags: */
+	riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
+
+	/*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary 
+	 *or a hole: */
+	nriftsegs=0;
+	for (i=0;i<nsegs;i++){
+		el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
+		/*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and 
+		 *then  proceed to find another element that owns the segment. If we don't find it, we know 
+		 *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
+		element_nodes[0]=*(index+3*el+0);
+		element_nodes[1]=*(index+3*el+1);
+		element_nodes[2]=*(index+3*el+2);
+
+		index[3*el+0]=-1;
+		index[3*el+1]=-1;
+		index[3*el+2]=-1;
+
+		el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel); 
+
+		/*Restore index: */
+		index[3*el+0]=element_nodes[0];
+		index[3*el+1]=element_nodes[1];
+		index[3*el+2]=element_nodes[2];
+
+		if (el2!=-1){
+			/*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
+		    riftsegments_uncompressed[5*i+0]=1;
+		    riftsegments_uncompressed[5*i+1]=el;
+		    riftsegments_uncompressed[5*i+2]=el2;
+		    riftsegments_uncompressed[5*i+3]=segments[3*i+0];
+			 riftsegments_uncompressed[5*i+4]=segments[3*i+1];
+			 nriftsegs++;
+		}
+	}
+
+	/*Compress riftsegments_uncompressed:*/
+	riftsegments=xNew<int>(nriftsegs*4);
+	counter=0;
+	for (i=0;i<nsegs;i++){
+		if (riftsegments_uncompressed[5*i+0]){
+			riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
+			riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
+			riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
+			riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
+			counter++;
+		}
+	}
+	xDelete<int>(riftsegments_uncompressed);
+
+	/*Assign output pointers: */
+	*priftsegments=riftsegments;
+	*pnriftsegs=nriftsegs;
+}/*}}}*/
+int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){/*{{{*/
+
+	int noerr=1;
+	int k,l,counter;
+	int newel;
+
+	int* GridElements=NULL;
+	int  NumGridElements;
+
+	/*Output: */
+	int NumGridElementListOnOneSideOfRift;
+	int* GridElementListOnOneSideOfRift=NULL;
+
+	/*Build a list of all the elements connected to this node: */
+	GridElementsList(&GridElements,&NumGridElements,node,index,nel);
+
+	/*Figure out the list of elements  that are on the same side of the rift. To do so, we start from one 
+	 * side of the rift and keep rotating in the same direction:*/
+	GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
+	//bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
+	GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there 
+															   for a rotation direction, we 'll take it out when we are 
+															   done rotating*/
+	GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
+	counter=1;
+	for (;;){
+		/*Find neighbour of element GridElementListOnOneSideOfRift[counter], not 
+		 * equal to GridElementListOnOneSideOfRift[counter-1]*/
+		for (k=0;k<NumGridElements;k++){
+			if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
+				/*Verify this element is not already in our list of element on the same side of the rift: */
+				newel=1;
+				for (l=0;l<=counter;l++){
+					if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
+						newel=0;
+						break;
+					}
+				}
+				if (newel){
+					counter++;
+					GridElementListOnOneSideOfRift[counter]=GridElements[k];
+					if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
+						break;
+					}
+					k=-1;
+				}
+			}
+		}
+		/*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
+		NumGridElementListOnOneSideOfRift=counter;
+		for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
+			GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
+		}
+		break;
+	}// for (;;)
+
+	/*Free ressources: */
+	xDelete<int>(GridElements);
+	/*Assign output pointers: */
+	*pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
+	*pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
+	return noerr;
+}/*}}}*/
+int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){/*{{{*/
+
+	int noerr=1;
+	int i,j,k;
+	int el1,el2;
+
+	int *segments          = NULL;
+	int *segmentmarkerlist = NULL;
+	int  nsegs;
+
+	/*Recover input: */
+	segments          = *psegments;
+	segmentmarkerlist = *psegmentmarkerlist;
+	nsegs             = *pnsegs;
+
+	/*Reallocate segments: */
+	segments         =xReNew<int>(segments,         nsegs*3,(nsegs+nriftsegs)*3);
+	segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
+
+	/*First, update the existing segments to the new nodes :*/
+	for (i=0;i<nriftsegs;i++){
+		el1=riftsegments[4*i+0];
+		el2=riftsegments[4*i+1];
+		for (j=0;j<nsegs;j++){
+			if (segments[3*j+2]==(el1+1)){
+				/*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
+				 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts, 
+				 *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
+				for (k=0;k<3;k++){
+					if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
+						*(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
+						break;
+					}
+				}
+				for (k=0;k<3;k++){
+					if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1])  && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
+						*(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
+						break;
+					}
+				}
+				/*Deal with el2: */
+				*(segments+3*(nsegs+i)+2)=el2+1;
+				*(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
+				for (k=0;k<3;k++){
+					if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
+						*(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
+						break;
+					}
+				}
+				for (k=0;k<3;k++){
+					if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
+						*(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
+						break;
+					}
+				}
+			}
+			if (*(segments+3*j+2)==(el2+1)){
+				/*segment j is the same as rift segment i.*/
+				/*Let's update segments[j][:] using  element el2 and the corresponding rift segment: */
+				for (k=0;k<3;k++){
+					if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
+						*(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
+						break;
+					}
+				}
+				for (k=0;k<3;k++){
+					if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
+						*(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
+						break;
+					}
+				}
+				/*Deal with el1: */
+				*(segments+3*(nsegs+i)+2)=el1+1;
+				*(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
+				for (k=0;k<3;k++){
+					if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
+						*(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
+						break;
+					}
+				}
+				for (k=0;k<3;k++){
+					if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
+						*(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
+						break;
+					}
+				}
+			}
+		}
+	}
+	nsegs+=nriftsegs;
+
+	/*Assign output pointers: */
+	*psegments=segments;
+	*psegmentmarkerlist=segmentmarkerlist;
+	*pnsegs=nsegs;
+
+	return noerr;
+}/*}}}*/
+int FindElement(int A,int B,int* index,int nel){/*{{{*/
+
+	int el=-1;
+	for (int n=0;n<nel;n++){
+		if(((index[3*n+0]==A) || (index[3*n+1]==A) || (index[3*n+2]==A)) && ((index[3*n+0]==B) || (index[3*n+1]==B) || (index[3*n+2]==B))){
+			el=n;
+			break;
+		}
+	}
+	return el;
+}/*}}}*/
+int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){/*{{{*/
+
+	/*Using segment markers, wring out the rift segments from the segments. Rift markers are 
+	 *of the form 2+i where i=0 to number of rifts */
+
+	int noerr=1;
+	int i,j,counter;
+
+	/*input: */
+	int *segments          = NULL;
+	int *segmentmarkerlist = NULL;
+	int numsegs;
+
+	/*output: */
+	int   new_numsegs;
+	int  *riftsnumsegs       = NULL;
+	int **riftssegments      = NULL;
+	int  *new_segments       = NULL;
+	int  *new_segmentmarkers = NULL;
+
+	/*intermediary: */
+	int* riftsegment=NULL;
+
+	/*Recover input arguments: */
+	segments          = *psegments;
+	numsegs           = *pnumsegs;
+	segmentmarkerlist = *psegmentmarkerlist;
+
+	/*First, figure out  how many segments will be left in 'segments': */
+	counter=0;
+	for (i=0;i<numsegs;i++){
+		if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
+	}
+	/*Allocate new segments: */
+	new_numsegs=counter;
+	new_segments=xNew<int>(new_numsegs*3);
+	new_segmentmarkers=xNew<int>(new_numsegs);
+
+	/*Copy new segments info : */
+	counter=0;
+	for (i=0;i<numsegs;i++){
+		if (segmentmarkerlist[i]==1){
+			new_segments[3*counter+0]=segments[3*i+0];
+			new_segments[3*counter+1]=segments[3*i+1];
+			new_segments[3*counter+2]=segments[3*i+2];
+			new_segmentmarkers[counter]=segmentmarkerlist[i];
+			counter++;
+		}
+	}
+
+	/*Now deal with rift segments: */
+	riftsnumsegs=xNew<int>(numrifts);
+	riftssegments=xNew<int*>(numrifts);
+	for (i=0;i<numrifts;i++){
+		/*Figure out how many segments for rift i: */
+		counter=0;
+		for (j=0;j<numsegs;j++){
+			if (segmentmarkerlist[j]==2+i)counter++;
+		}
+		riftsnumsegs[i]=counter;
+		riftsegment=xNew<int>(counter*3);
+		/*Copy new segments info :*/
+		counter=0;
+		for (j=0;j<numsegs;j++){
+			if (segmentmarkerlist[j]==(2+i)){
+				riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods+1);
+				riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods+1);
+				riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel+1);
+				counter++;
+			}
+		}
+		*(riftssegments+i)=riftsegment;
+	}
+
+	/*Free ressources: */
+	xDelete<int>(segments);
+
+	/*Assign output pointers: */
+	*psegments=new_segments;
+	*psegmentmarkerlist=new_segmentmarkers;
+	*pnumsegs=new_numsegs;
+	*pnumrifts=numrifts;
+	*priftssegments=riftssegments;
+	*priftsnumsegs=riftsnumsegs;
+	return noerr;
+}/*}}}*/
+int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){/*{{{*/
+
+	int noerr=1;
+	int i,j,k;
+
+	/*output: */
+	int  *riftsnumpairs = NULL;
+	int **riftspairs    = NULL;
+
+	/*intermediary :*/
+	int  numsegs;
+	int* segments=NULL;
+	int* pairs=NULL;
+	int  node1,node2,node3,node4;
+
+	riftsnumpairs=xNew<int>(numrifts);
+	riftspairs=xNew<int*>(numrifts);
+	for (i=0;i<numrifts;i++){
+		segments=riftssegments[i];
+		numsegs =riftsnumsegments[i];
+		riftsnumpairs[i]=numsegs;
+		pairs=xNew<int>(2*numsegs);
+		for (j=0;j<numsegs;j++){
+			pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
+			node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
+			/*Find element facing on other side of rift: */
+			for (k=0;k<numsegs;k++){
+				if (k==j)continue;
+				node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
+				/*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
+				if (   ((x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]))
+				    || ((x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1]))  ){
+					/*We found the corresponding element: */
+					pairs[2*j+1]=segments[3*k+2];
+					break;
+				}
+			}
+		}
+		riftspairs[i]=pairs;
+	}
+
+	/*Assign output pointers: */
+	*priftsnumpairs=riftsnumpairs;
+	*priftspairs=riftspairs;
+	return noerr;
+}/*}}}*/
+int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){/*{{{*/
+
+	int i;
+	int noerr=1;
+
+	/*output: */
+	int riftflag=0;
+	int numrifts=0;
+
+	int maxmark=1; //default marker for regular segments
+
+	/*Any marker >=2 indicates a certain rift: */
+	numrifts=0;
+	for (i=0;i<nsegs;i++){
+		if (segmentmarkerlist[i]>maxmark){
+			numrifts++;
+			maxmark=segmentmarkerlist[i];
+		}
+	}
+	if(numrifts)riftflag=1;
+
+	/*Assign output pointers:*/
+	*priftflag=riftflag;
+	*pnumrifts=numrifts;
+	return noerr;
+}/*}}}*/
+int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){/*{{{*/
+
+	int noerr=1;
+	int i,j,k,counter;
+
+	/*intermediary: */
+	int *riftsegments = NULL;
+	int *riftpairs    = NULL;
+	int numsegs;
+
+	/*ordering and copy: */
+	int *order             = NULL;
+	int *riftsegments_copy = NULL;
+	int *riftpairs_copy    = NULL;
+
+	/*node and element manipulation: */
+	int node1,node2,node3,node4,temp_node,tip1,tip2,node;
+	int el2;
+	int already_ordered=0;
+
+	/*output: */
+	int* riftstips=NULL;
+
+	/*Allocate byproduct of this routine, riftstips: */
+	riftstips=xNew<int>(numrifts*2);
+
+	/*Go through all rifts: */
+	for (i=0;i<numrifts;i++){
+		riftsegments = riftssegments[i];
+		riftpairs    = riftspairs[i];
+		numsegs      = riftsnumsegments[i];
+
+		/*Allocate copy of riftsegments and riftpairs, 
+		 *as well as ordering vector: */
+		riftsegments_copy=xNew<int>(numsegs*3);
+		riftpairs_copy=xNew<int>(numsegs*2);
+		order=xNew<int>(numsegs);
+
+		/*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
+		tip1=-1;
+		tip2=-1;
+
+		for (j=0;j<numsegs;j++){
+			el2=*(riftpairs+2*j+1);
+			node1=*(riftsegments+3*j+0);
+			node2=*(riftsegments+3*j+1);
+			/*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and 
+			 *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
+			for (k=0;k<numsegs;k++){
+				if (*(riftsegments+3*k+2)==el2){
+					node3=*(riftsegments+3*k+0);
+					node4=*(riftsegments+3*k+1);
+					break;
+				}
+			}
+			/* Make sure node3 faces node1 and node4 faces node2: */
+			_assert_(node1<nods+1 && node4<nods+1);
+			_assert_(node1>0 && node4>0);
+			if ((x[node1-1]==x[node4-1]) && (y[node1-1]==y[node4-1])){
+				/*Swap node3 and node4:*/
+				temp_node=node3;
+				node3=node4;
+				node4=temp_node;
+			}
+
+			/*Figure out if a tip is on this element: */
+			if (node3==node1){
+				/*node1 is a tip*/
+				if (tip1==-1) {
+					tip1=node1;
+					continue;
+				}
+				if ((tip2==-1) && (node1!=tip1)){
+					tip2=node1;
+					break;
+				}
+			}
+
+			if (node4==node2){
+				/*node2 is a tip*/
+				if (tip1==-1){
+					tip1=node2;
+					continue;
+				}
+				if ((tip2==-1) && (node2!=tip1)){
+					tip2=node2;
+					break;
+				}
+			}
+		}
+
+		/*Record tips in riftstips: */
+		*(riftstips+2*i+0)=tip1;
+		*(riftstips+2*i+1)=tip2;
+
+		/*We have the two tips for this rift.  Go from tip1 to tip2, and figure out the order in which segments are sequential. 
+		 *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
+		node=tip1;
+		for (counter=0;counter<numsegs;counter++){
+			for (j=0;j<numsegs;j++){
+				node1=*(riftsegments+3*j+0);
+				node2=*(riftsegments+3*j+1);
+
+				if ((node1==node) || (node2==node)){
+					/*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
+					already_ordered=0;
+					for (k=0;k<counter;k++){
+						if(order[k]==j){
+							already_ordered=1;
+							break;
+						}
+					}
+					if (!already_ordered){
+						order[counter]=j;
+						if(node1==node){
+							node=node2;
+						}
+						else if(node2==node){
+							node=node1;
+						}
+						break;
+					}
+				}
+			}
+		}
+
+		/*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
+		for (j=0;j<numsegs;j++){
+			_assert_(order[j]<numsegs);
+			*(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
+			*(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
+			*(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
+			*(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
+			*(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
+		}
+
+		for (j=0;j<numsegs;j++){
+			*(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
+			*(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
+			*(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
+			*(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
+			*(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
+		}
+
+		xDelete<int>(order);
+		xDelete<int>(riftsegments_copy);
+		xDelete<int>(riftpairs_copy);
+
+	}
+
+	/*Assign output pointer:*/
+	*priftstips=riftstips;
+	return noerr;
+}/*}}}*/
+int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,/*{{{*/
+		int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
+
+	int noerr=1;
+	int i,j,k,k0;
+
+	double el1,el2,node1,node2,node3,node4;
+	double temp_node;
+
+	/*output: */
+	double **riftspenaltypairs    = NULL;
+	double  *riftpenaltypairs     = NULL;
+	int     *riftsnumpenaltypairs = NULL;
+
+	/*intermediary: */
+	int numsegs;
+	int* riftsegments=NULL;
+	int* riftpairs=NULL;
+	int counter;
+	double normal[2];
+	double length;
+	int    k1,k2;
+
+	/*Allocate: */
+	riftspenaltypairs=xNew<double*>(numrifts);
+	riftsnumpenaltypairs=xNew<int>(numrifts);
+
+	for(i=0;i<numrifts;i++){
+		numsegs=riftsnumsegs[i];
+		riftsegments=riftssegments[i];
+		riftpairs=riftspairs[i];
+
+		/*allocate riftpenaltypairs, and riftnumpenaltypairs: */
+		if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
+
+		/*Go through only one flank of the rifts, not counting the tips: */
+		counter=0;
+		for(j=0;j<(numsegs/2);j++){
+			el1=*(riftpairs+2*j+0);
+			el2=*(riftpairs+2*j+1);
+			node1=*(riftsegments+3*j+0);
+			node2=*(riftsegments+3*j+1);
+			/*Find segment index to recover node3 and node4, facing node1 and node2: */
+			k0=-1;
+			for(k=0;k<numsegs;k++){
+				if(*(riftsegments+3*k+2)==el2){
+					k0=k;
+					break;
+				}
+			}
+			node3=*(riftsegments+3*k0+0);
+			node4=*(riftsegments+3*k0+1);
+
+			/* Make sure node3 faces node1 and node4 faces node2: */
+			if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
+				/*Swap node3 and node4:*/
+				temp_node=node3;
+				node3=node4;
+				node4=temp_node;
+			}	
+			/*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to 
+			 *this segment, and its length: */
+			normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
+			normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
+			length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
+
+			/*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
+			 * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4, 
+			 * only once. We'll add the normals and the lengths : */
+
+			if(node1!=node3){ //exclude tips from loads
+				k1=-1;
+				for(k=0;k<counter;k++){
+					if( (*(riftpenaltypairs+k*7+0))==node1){
+						k1=k; 
+						break;
+					}
+				}
+				if(k1==-1){
+					*(riftpenaltypairs+counter*7+0)=node1;
+					*(riftpenaltypairs+counter*7+1)=node3;
+					*(riftpenaltypairs+counter*7+2)=el1;
+					*(riftpenaltypairs+counter*7+3)=el2;
+					*(riftpenaltypairs+counter*7+4)=normal[0];
+					*(riftpenaltypairs+counter*7+5)=normal[1];
+					*(riftpenaltypairs+counter*7+6)=length/2;
+					counter++;
+				}
+				else{
+					*(riftpenaltypairs+k1*7+4)+=normal[0];
+					*(riftpenaltypairs+k1*7+5)+=normal[1];
+					*(riftpenaltypairs+k1*7+6)+=length/2;
+				}
+			}
+			if(node2!=node4){
+				k2=-1;
+				for(k=0;k<counter;k++){
+					if( (*(riftpenaltypairs+k*7+0))==node2){
+						k2=k;
+						break;
+					}
+				}
+				if(k2==-1){
+					*(riftpenaltypairs+counter*7+0)=node2;
+					*(riftpenaltypairs+counter*7+1)=node4;
+					*(riftpenaltypairs+counter*7+2)=el1;
+					*(riftpenaltypairs+counter*7+3)=el2;
+					*(riftpenaltypairs+counter*7+4)=normal[0];
+					*(riftpenaltypairs+counter*7+5)=normal[1];
+					*(riftpenaltypairs+counter*7+6)=length/2;
+					counter++;
+				}
+				else{
+					*(riftpenaltypairs+k2*7+4)+=normal[0];
+					*(riftpenaltypairs+k2*7+5)+=normal[1];
+					*(riftpenaltypairs+k2*7+6)+=length/2;
+				}
+			}
+		}
+		/*Renormalize normals: */
+		for(j=0;j<counter;j++){
+			double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
+			*(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
+			*(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
+		}
+
+		riftspenaltypairs[i]=riftpenaltypairs;
+		riftsnumpenaltypairs[i]=(numsegs/2-1);
+	}
+
+	/*Assign output pointers: */
+	*priftspenaltypairs=riftspenaltypairs;
+	*priftsnumpenaltypairs=riftsnumpenaltypairs;
+	return noerr;
+}/*}}}*/
+int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
+
+	int noerr=1;
+	int i,j,k;
+	int node1,node2,node3;
+	int el;
+	double  pair[2];
+	int     pair_count=0;
+	int     triple=0;
+
+	/*Recover input: */
+	int    *index = *pindex;
+	int     nel   = *pnel;
+	double *x     = *px;
+	double *y     = *py;
+	int     nods  = *pnods;
+
+	for (i=0;i<num_seg;i++){
+		node1=*(segments+3*i+0);
+		node2=*(segments+3*i+1);
+		/*Find all elements connected to [node1 node2]: */
+		pair_count=0;
+		for (j=0;j<nel;j++){
+			if (*(index+3*j+0)==node1){
+				if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
+					pair[pair_count]=j;
+					pair_count++;
+				}
+			}
+			if (*(index+3*j+1)==node1){
+				if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
+					pair[pair_count]=j;
+					pair_count++;
+				}
+			}
+			if (*(index+3*j+2)==node1){
+				if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
+					pair[pair_count]=j;
+					pair_count++;
+				}
+			}
+		}
+		/*Ok, we have pair_count elements connected to this segment. For each of these elements, 
+		 *figure out if the third node also belongs to a segment: */
+		if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to  2 elements
+			continue;
+		}
+		else{
+			for (j=0;j<pair_count;j++){
+				el=(int)pair[j];
+				triple=0;
+				/*First find node3: */
+				if (*(index+3*el+0)==node1){
+					if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
+					else node3=*(index+3*el+1);
+				}
+				if (*(index+3*el+1)==node1){
+					if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
+					else node3=*(index+3*el+0);
+				}
+				if (*(index+3*el+2)==node1){
+					if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
+					else node3=*(index+3*el+0);
+				}
+				/*Ok, we have node3. Does node3 belong to a segment? : */
+				for (k=0;k<num_seg;k++){
+					if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
+						triple=1;
+						break;
+					}
+				}
+				if(triple==1){
+					/*el is a corner element: we need to split it in 3 triangles: */
+					x=xReNew<double>(x,nods,nods+1);
+					y=xReNew<double>(y,nods,nods+1);
+					x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
+					y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
+					index=xReNew<int>(index,nel*3,(nel+2*3));
+					/*First, reassign element el: */
+					*(index+3*el+0)=node1;
+					*(index+3*el+1)=node2;
+					*(index+3*el+2)=nods+1;
+					/*Other two elements: */
+					*(index+3*nel+0)=node2;
+					*(index+3*nel+1)=node3;
+					*(index+3*nel+2)=nods+1;
+
+					*(index+3*(nel+1)+0)=node3;
+					*(index+3*(nel+1)+1)=node1;
+					*(index+3*(nel+1)+2)=nods+1;
+					/*we need  to change the segment elements corresponding to el: */
+					for (k=0;k<num_seg;k++){
+						if (*(segments+3*k+2)==(el+1)){
+							if ( ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node2)) || ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node1))) *(segments+3*k+2)=el+1;
+							if ( ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node3)) || ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node2))) *(segments+3*k+2)=nel+1;
+							if ( ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node1)) || ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node3))) *(segments+3*k+2)=nel+2;
+						}
+					}
+
+					nods=nods+1;
+					nel=nel+2;
+					i=0;
+					break;
+				}
+			} //for (j=0;j<pair_count;j++)
+		}
+	}// for (i=0;i<num_seg;i++)
+
+	/*Assign output pointers: */
+	*pindex=index;
+	*pnel=nel;
+	*px=x;
+	*py=y;
+	*pnods=nods;
+	return noerr;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/shared/Triangle/triangle.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Triangle/triangle.h	(revision 22498)
+++ /issm/trunk-jpl/src/c/shared/Triangle/triangle.h	(revision 22498)
@@ -0,0 +1,33 @@
+/*!\file:  triangle.h
+ * \brief
+ */ 
+
+#ifndef _SHARED_TRIANGLE_H
+#define _SHARED_TRIANGLE_H
+
+#include <stdio.h>
+#include <math.h>
+
+//#define REAL double //took  it out because it may conflict with stdlib.h defines. put back if necessary
+int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel);
+int OrderSegments(int** psegments,int nseg, int* index,int nel);
+int GridInsideHole(double* px0,double* py0,int n,double* x,double* y);
+int FindElement(int A,int B,int* index,int nel);
+int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist);
+int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
+int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
+int IsNeighbor(int el1,int el2,int* index);
+int IsOnRift(int el,int nriftsegs,int* riftsegments);
+void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index, int nsegs,int* segments);
+int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift,int segmentnumber, int nriftsegs,int* riftsegments, int node,int* index,int nel);
+int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
+int FindElement(double A,double B,int* index,int nel);
+int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
+int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels);
+int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels);
+int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int**  riftssegments,
+		int* riftsnumsegments,int** riftspairs,int* riftstips,double* x,double* y);
+int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg);
+int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y);
+
+#endif  /* _SHARED_TRIANGLE_H */
Index: /issm/trunk-jpl/src/c/shared/shared.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/shared.h	(revision 22497)
+++ /issm/trunk-jpl/src/c/shared/shared.h	(revision 22498)
@@ -19,5 +19,5 @@
 #include "./String/sharedstring.h"
 #include "./Threads/issm_threads.h"
-#include "./TriMesh/trimesh.h"
+#include "./Triangle/triangle.h"
 #include "./LatLong/latlong.h"
 
Index: /issm/trunk-jpl/src/m/mesh/argusmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/argusmesh.m	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/argusmesh.m	(revision 22498)
@@ -9,5 +9,5 @@
 %
 %   Example:
-%     md=argusmesh(md,'TriMesh.exp')
+%     md=argusmesh(md,'Domain.exp')
 
 %some argument check: 
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m	(revision 22498)
@@ -25,7 +25,7 @@
 
 %Call MEX file
-[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers);
+[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=ProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers);
 if ~isstruct(md.rifts.riftstruct),
-	error('TriMeshProcessRifts did not find any rift');
+	error('ProcessRifts did not find any rift');
 end
 
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py	(revision 22498)
@@ -1,4 +1,4 @@
 import numpy as np
-from TriMeshProcessRifts import TriMeshProcessRifts
+from ProcessRifts import ProcessRifts
 from ContourToMesh import ContourToMesh
 from meshprocessoutsiderifts import meshprocessoutsiderifts
@@ -22,5 +22,5 @@
 
 	#Call MEX file
-	md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers)
+	md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct=ProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers)
 	md.mesh.elements=md.mesh.elements.astype(int)
 	md.mesh.x=md.mesh.x.reshape(-1)
@@ -29,5 +29,5 @@
 	md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
 	if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
-		raise RuntimeError("TriMeshProcessRifts did not find any rift")
+		raise RuntimeError("ProcessRifts did not find any rift")
 
 	#Fill in rest of fields:
Index: /issm/trunk-jpl/src/m/mesh/triangle.js
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.js	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/triangle.js	(revision 22498)
@@ -2,5 +2,5 @@
 //TRIANGLE - create model mesh using the triangle package
 //
-//   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
+//   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
 //   where md is a @model object, domainname is the name of an Argus domain outline file, 
 //   and resolution is a characteristic length for the mesh (same unit as the domain outline
@@ -36,5 +36,5 @@
 
 	//Call mesher: 
-	var return_array=TriMesh(md, domain, rifts, area); 
+	var return_array=Triangle(md, domain, rifts, area); 
 
 	//Plug into md:
Index: /issm/trunk-jpl/src/m/mesh/triangle.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.m	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/triangle.m	(revision 22498)
@@ -2,5 +2,5 @@
 %TRIANGLE - create model mesh using the triangle package
 %
-%   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
+%   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
 %   where md is a @model object, domainname is the name of an Argus domain outline file, 
 %   and resolution is a characteristic length for the mesh (same unit as the domain outline
@@ -42,6 +42,6 @@
 end
 
-%Mesh using TriMesh
-[elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area);
+%Mesh using Triangle
+[elements,x,y,segments,segmentmarkers]=Triangle(domainname,riftname,area);
 
 %check that all the created nodes belong to at least one element
Index: /issm/trunk-jpl/src/m/mesh/triangle.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 22498)
@@ -2,5 +2,5 @@
 import numpy as np
 from mesh2d import mesh2d
-from TriMesh import TriMesh
+from Triangle import Triangle
 from NodeConnectivity import NodeConnectivity
 from ElementConnectivity import ElementConnectivity
@@ -11,5 +11,5 @@
 	TRIANGLE - create model mesh using the triangle package
 
-	   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
+	   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
 	   where md is a @model object, domainname is the name of an Argus domain outline file, 
 	   and resolution is a characteristic length for the mesh (same unit as the domain outline
@@ -48,7 +48,7 @@
 		raise IOError("file '%s' not found" % domainname)
 
-	#Mesh using TriMesh
+	#Mesh using Triangle
 	md.mesh=mesh2d()
-	md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=TriMesh(domainname,riftname,area)
+	md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=Triangle(domainname,riftname,area)
 	md.mesh.elements=md.mesh.elements.astype(int)
 	md.mesh.segments=md.mesh.segments.astype(int)
Index: /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m	(revision 22497)
+++ /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m	(revision 22498)
@@ -2,5 +2,5 @@
 %TRIANGLE - create model mesh using the triangle package
 %
-%   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
+%   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
 %   where md is a @model object, domainname is the name of an Argus domain outline file, 
 %   and resolution is a characteristic length for the mesh (same unit as the domain outline
@@ -24,6 +24,6 @@
 area=resolution^2;
 
-%Mesh using TriMesh
-[elements,x,z,segments,segmentmarkers]=TriMesh(domainname,'',area);
+%Mesh using Triangle
+[elements,x,z,segments,segmentmarkers]=Triangle(domainname,'',area);
 
 %check that all the created nodes belong to at least one element
Index: /issm/trunk-jpl/src/m/modules/ProcessRifts.m
===================================================================
--- /issm/trunk-jpl/src/m/modules/ProcessRifts.m	(revision 22498)
+++ /issm/trunk-jpl/src/m/modules/ProcessRifts.m	(revision 22498)
@@ -0,0 +1,17 @@
+function [index2,x2,y2,segments2,segmentmarkers2,rifts2] = ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
+%TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
+%
+%   Usage: 
+%      [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1); 
+%   
+%   (index1,x1,y1,segments1,segmentmarkers1):	An initial triangulation.
+%   [index2,x2,y2,segments2,segmentmarkers2,rifts2]:	The resulting triangulation where rifts have been processed.
+
+% Check usage
+if nargin~=5 && nargout~=6
+	help ProcessRifts
+	error('Wrong usage (see above)');
+end
+
+% Call mex module
+[index2,x2,y2,segments2,segmentmarkers2,rifts2] = ProcessRifts_matlab(index1,x1,y1,segments1,segmentmarkers1);
Index: /issm/trunk-jpl/src/m/modules/ProcessRifts.py
===================================================================
--- /issm/trunk-jpl/src/m/modules/ProcessRifts.py	(revision 22498)
+++ /issm/trunk-jpl/src/m/modules/ProcessRifts.py	(revision 22498)
@@ -0,0 +1,16 @@
+from ProcessRifts_python import ProcessRifts_python
+
+def ProcessRifts(index1,x1,y1,segments1,segmentmarkers1):
+	"""
+	TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
+
+	   Usage: 
+		   [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
+
+	   (index1,x1,y1,segments1,segmentmarkers1):	An initial triangulation.
+	   [index2,x2,y2,segments2,segmentmarkers2,rifts2]:	The resulting triangulation where rifts have been processed.
+	"""
+	# Call mex module
+	index2,x2,y2,segments2,segmentmarkers2,rifts2 = ProcessRifts_python(index1,x1,y1,segments1,segmentmarkers1)
+	# Return
+	return index2,x2,y2,segments2,segmentmarkers2,rifts2
Index: sm/trunk-jpl/src/m/modules/TriMesh.m
===================================================================
--- /issm/trunk-jpl/src/m/modules/TriMesh.m	(revision 22497)
+++ 	(revision )
@@ -1,21 +1,0 @@
-function [index,x,y,segments,segmentmarkers] = TriMesh(domainoutlinefilename,rifts,mesh_area);
-%TRIMESH - Mesh a domain using an .exp file
-%
-%   Usage: 
-%     [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,mesh_area);
-%	      
-%   index,x,y:	Defines a triangulation 
-%   segments:	Array made of exterior segments to the mesh domain outline 
-%   segmentmarkers:	Array flagging each segment
-%
-%   domainoutlinefilename:	Argus domain outline file
-%   mesh_area:	Maximum area desired for any element of the resulting mesh
-
-% Check usage
-if nargin~=3 && nargout~=5
-	help TriMesh
-	error('Wrong usage (see above)');
-end
-
-% Call mex module
-[index,x,y,segments,segmentmarkers]=TriMesh_matlab(domainoutlinefilename,rifts,mesh_area);
Index: sm/trunk-jpl/src/m/modules/TriMesh.py
===================================================================
--- /issm/trunk-jpl/src/m/modules/TriMesh.py	(revision 22497)
+++ 	(revision )
@@ -1,21 +1,0 @@
-from TriMesh_python import TriMesh_python
-
-def TriMesh(domainoutlinefilename,rifts,mesh_area):
-	"""
-	TRIMESH - Mesh a domain using an .exp file
-
-	   Usage: 
-			[index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,mesh_area); 
-
-	   index,x,y: defines a triangulation 
-		segments: An array made of exterior segments to the mesh domain outline 
-		segmentmarkers: An array flagging each segment
-
-	   domainoutlinefilename: an Argus domain outline file
-		mesh_area: The maximum area desired for any element of the resulting mesh
-	"""
-	# Call mex module
-	index,x,y,segments,segmentmarkers=TriMesh_python(domainoutlinefilename,rifts,mesh_area)
-	# Return
-	return index,x,y,segments,segmentmarkers
-
Index: sm/trunk-jpl/src/m/modules/TriMeshProcessRifts.m
===================================================================
--- /issm/trunk-jpl/src/m/modules/TriMeshProcessRifts.m	(revision 22497)
+++ 	(revision )
@@ -1,17 +1,0 @@
-function [index2,x2,y2,segments2,segmentmarkers2,rifts2] = TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
-%TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
-%
-%   Usage: 
-%      [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1); 
-%   
-%   (index1,x1,y1,segments1,segmentmarkers1):	An initial triangulation.
-%   [index2,x2,y2,segments2,segmentmarkers2,rifts2]:	The resulting triangulation where rifts have been processed.
-
-% Check usage
-if nargin~=5 && nargout~=6
-	help TriMeshProcessRifts
-	error('Wrong usage (see above)');
-end
-
-% Call mex module
-[index2,x2,y2,segments2,segmentmarkers2,rifts2] = TriMeshProcessRifts_matlab(index1,x1,y1,segments1,segmentmarkers1);
Index: sm/trunk-jpl/src/m/modules/TriMeshProcessRifts.py
===================================================================
--- /issm/trunk-jpl/src/m/modules/TriMeshProcessRifts.py	(revision 22497)
+++ 	(revision )
@@ -1,16 +1,0 @@
-from TriMeshProcessRifts_python import TriMeshProcessRifts_python
-
-def TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1):
-	"""
-	TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
-
-	   Usage: 
-		   [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
-
-	   (index1,x1,y1,segments1,segmentmarkers1):	An initial triangulation.
-	   [index2,x2,y2,segments2,segmentmarkers2,rifts2]:	The resulting triangulation where rifts have been processed.
-	"""
-	# Call mex module
-	index2,x2,y2,segments2,segmentmarkers2,rifts2 = TriMeshProcessRifts_python(index1,x1,y1,segments1,segmentmarkers1)
-	# Return
-	return index2,x2,y2,segments2,segmentmarkers2,rifts2
Index: /issm/trunk-jpl/src/m/modules/Triangle.m
===================================================================
--- /issm/trunk-jpl/src/m/modules/Triangle.m	(revision 22498)
+++ /issm/trunk-jpl/src/m/modules/Triangle.m	(revision 22498)
@@ -0,0 +1,21 @@
+function [index,x,y,segments,segmentmarkers] = Triangle(domainoutlinefilename,rifts,mesh_area);
+%TRIMESH - Mesh a domain using an .exp file
+%
+%   Usage: 
+%     [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,mesh_area);
+%	      
+%   index,x,y:	Defines a triangulation 
+%   segments:	Array made of exterior segments to the mesh domain outline 
+%   segmentmarkers:	Array flagging each segment
+%
+%   domainoutlinefilename:	Argus domain outline file
+%   mesh_area:	Maximum area desired for any element of the resulting mesh
+
+% Check usage
+if nargin~=3 && nargout~=5
+	help Triangle
+	error('Wrong usage (see above)');
+end
+
+% Call mex module
+[index,x,y,segments,segmentmarkers]=Triangle_matlab(domainoutlinefilename,rifts,mesh_area);
Index: /issm/trunk-jpl/src/m/modules/Triangle.py
===================================================================
--- /issm/trunk-jpl/src/m/modules/Triangle.py	(revision 22498)
+++ /issm/trunk-jpl/src/m/modules/Triangle.py	(revision 22498)
@@ -0,0 +1,21 @@
+from Triangle_python import Triangle_python
+
+def Triangle(domainoutlinefilename,rifts,mesh_area):
+	"""
+	TRIMESH - Mesh a domain using an .exp file
+
+	   Usage: 
+			[index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,mesh_area); 
+
+	   index,x,y: defines a triangulation 
+		segments: An array made of exterior segments to the mesh domain outline 
+		segmentmarkers: An array flagging each segment
+
+	   domainoutlinefilename: an Argus domain outline file
+		mesh_area: The maximum area desired for any element of the resulting mesh
+	"""
+	# Call mex module
+	index,x,y,segments,segmentmarkers=Triangle_python(domainoutlinefilename,rifts,mesh_area)
+	# Return
+	return index,x,y,segments,segmentmarkers
+
Index: /issm/trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp	(revision 22498)
@@ -0,0 +1,59 @@
+/*!\file:  ProcessRifts.cpp
+ * \brief split a mesh where a rift (or fault) is present
+ */ 
+
+#include "./ProcessRifts.h"
+
+void ProcessRiftsUsage(void){/*{{{*/
+	_printf_("\n");
+	_printf_("   usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1) \n");
+	_printf_("      where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
+	_printf_("      index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
+}/*}}}*/
+WRAPPER(ProcessRifts_python){
+
+	/* returned quantities: */
+	RiftStruct *riftstruct = NULL;
+
+	/* input: */
+	int     nel,nods;
+	int    *index          = NULL;
+	double *x              = NULL;
+	double *y              = NULL;
+	int    *segments       = NULL;
+	int    *segmentmarkers = NULL;
+	int     num_seg;
+
+	/*Boot module*/
+	MODULEBOOT();
+
+	/*checks on arguments on the matlab side: */
+	CHECKARGUMENTS(NLHS,NRHS,&ProcessRiftsUsage);
+
+	/*Fetch data */
+	FetchData(&index,&nel,NULL,INDEXIN);
+	FetchData(&x,&nods,XIN);
+	FetchData(&y,NULL,YIN);
+	FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
+	FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
+
+	/*call x layer*/
+	ProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,&riftstruct);
+
+	/*Output : */
+	WriteData(INDEXOUT,index,nel,3);
+	WriteData(XOUT,x,nods,1);
+	WriteData(YOUT,y,nods,1);
+	WriteData(SEGMENTSOUT,segments,num_seg,3);
+	WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
+	WriteData(RIFTSTRUCT,riftstruct);
+
+	/*end module: */
+	delete riftstruct;
+	xDelete<int>(index);
+	xDelete<double>(x);
+	xDelete<double>(y);
+	xDelete<int>(segments);
+	xDelete<int>(segmentmarkers );
+	MODULEEND();
+}
Index: /issm/trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h
===================================================================
--- /issm/trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h	(revision 22498)
+++ /issm/trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h	(revision 22498)
@@ -0,0 +1,65 @@
+/*
+ * ProcessRifts.h
+ */ 
+
+#ifndef _PROCESSRIFTS_H_
+#define _PROCESSRIFTS_H_
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+/*For python modules: needs to come before header files inclusion*/
+#ifdef _HAVE_PYTHON_
+#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
+#endif
+
+#include "../bindings.h"
+#include "../../c/main/globals.h"
+#include "../../c/modules/modules.h"
+#include "../../c/shared/shared.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__  "ProcessRifts"
+
+#ifdef _HAVE_MATLAB_MODULES_
+/* serial input macros: */
+#define INDEXIN          prhs[0]
+#define XIN              prhs[1]
+#define YIN              prhs[2]
+#define SEGMENTSIN       prhs[3]
+#define SEGMENTMARKERSIN prhs[4]
+/* serial output macros: */
+#define INDEXOUT          (mxArray**)&plhs[0]
+#define XOUT              (mxArray**)&plhs[1]
+#define YOUT              (mxArray**)&plhs[2]
+#define SEGMENTSOUT       (mxArray**)&plhs[3]
+#define SEGMENTMARKERSOUT (mxArray**)&plhs[4]
+#define RIFTSTRUCT        (mxArray**)&plhs[5]
+#endif
+
+#ifdef _HAVE_PYTHON_MODULES_
+/* serial input macros: */
+#define INDEXIN          PyTuple_GetItem(args,0)
+#define XIN              PyTuple_GetItem(args,1)
+#define YIN              PyTuple_GetItem(args,2)
+#define SEGMENTSIN       PyTuple_GetItem(args,3)
+#define SEGMENTMARKERSIN PyTuple_GetItem(args,4)
+/* serial output macros: */
+#define INDEXOUT          output,0
+#define XOUT              output,1
+#define YOUT              output,2
+#define SEGMENTSOUT       output,3
+#define SEGMENTMARKERSOUT output,4
+#define RIFTSTRUCT        output,5
+#endif
+
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  6
+#undef NRHS
+#define NRHS  5
+
+#endif
Index: /issm/trunk-jpl/src/wrappers/Triangle/Triangle.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/Triangle/Triangle.cpp	(revision 22498)
+++ /issm/trunk-jpl/src/wrappers/Triangle/Triangle.cpp	(revision 22498)
@@ -0,0 +1,63 @@
+/*
+ * Triangle: mesh a domain using an .exp file
+ */
+
+#include "./Triangle.h"
+
+void TriangleUsage(void){/*{{{*/
+	_printf_("\n");
+	_printf_("   usage: [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,area) \n");
+	_printf_("      where: index,x,y defines a triangulation, segments is an array made \n");
+	_printf_("      of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment, \n");
+	_printf_("      outlinefilename an Argus domain outline file, \n");
+	_printf_("      area is the maximum area desired for any element of the resulting mesh, \n");
+	_printf_("\n");
+}/*}}}*/
+WRAPPER(Triangle_python){
+	
+	/*intermediary: */
+	double    area;
+	Contours *domain = NULL;
+	Contours *rifts  = NULL;
+
+	/* output: */
+	int    *index             = NULL;
+	double *x                 = NULL;
+	double *y                 = NULL;
+	int    *segments          = NULL;
+	int    *segmentmarkerlist = NULL;
+	int     nel,nods,nsegs;
+
+	/*Boot module: */
+	MODULEBOOT();
+
+	/*checks on arguments: */
+	CHECKARGUMENTS(NLHS,NRHS,&TriangleUsage);
+
+	/*Fetch data needed for meshing: */
+	FetchData(&domain,DOMAINOUTLINE);
+	FetchData(&rifts,RIFTSOUTLINE);
+	FetchData(&area,AREA);
+
+	/*call x core: */
+	Trianglex(&index,&x,&y,&segments,&segmentmarkerlist,&nel,&nods,&nsegs,domain,rifts,area);
+
+	/*write outputs: */
+	WriteData(INDEX,index,nel,3);
+	WriteData(X,x,nods);
+	WriteData(Y,y,nods);
+	WriteData(SEGMENTS,segments,nsegs,3);
+	WriteData(SEGMENTMARKERLIST,segmentmarkerlist,nsegs);
+
+	/*free ressources: */
+	delete domain;
+	delete rifts;
+	xDelete<int>(index);
+	xDelete<double>(x);
+	xDelete<double>(y);
+	xDelete<int>(segments);
+	xDelete<int>(segmentmarkerlist);
+
+	/*end module: */
+	MODULEEND();
+}
Index: /issm/trunk-jpl/src/wrappers/Triangle/Triangle.h
===================================================================
--- /issm/trunk-jpl/src/wrappers/Triangle/Triangle.h	(revision 22498)
+++ /issm/trunk-jpl/src/wrappers/Triangle/Triangle.h	(revision 22498)
@@ -0,0 +1,83 @@
+/*
+	Triangle.h
+*/
+
+#ifndef _TRIANGLE_H
+#define _TRIANGLE_H
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+	#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+/*For python modules: needs to come before header files inclusion*/
+#ifdef _HAVE_PYTHON_
+#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
+#endif
+
+#ifdef _HAVE_JAVASCRIPT_MODULES_
+#undef _DO_NOT_LOAD_GLOBALS_ /*only module where this needs to be undefined, so as to 
+							   not include IssmComm several times in the javascript Modle construct.*/
+#endif
+
+/*Header files: */
+#include "../bindings.h"
+#include "../../c/main/globals.h"
+#include "../../c/toolkits/toolkits.h"
+#include "../../c/modules/modules.h"
+#include "../../c/shared/shared.h"
+#include "../../c/shared/io/io.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__  "Triangle"
+
+#ifdef _HAVE_MATLAB_MODULES_
+/* serial input macros: */
+#define DOMAINOUTLINE  prhs[0]
+#define RIFTSOUTLINE   prhs[1]
+#define AREA           prhs[2]
+/* serial output macros: */
+#define INDEX             (mxArray**)&plhs[0]
+#define X                 (mxArray**)&plhs[1]
+#define Y                 (mxArray**)&plhs[2]
+#define SEGMENTS          (mxArray**)&plhs[3]
+#define SEGMENTMARKERLIST (mxArray**)&plhs[4]
+#endif
+
+#ifdef _HAVE_PYTHON_MODULES_
+/* serial input macros: */
+#define DOMAINOUTLINE PyTuple_GetItem(args,0)
+#define RIFTSOUTLINE  PyTuple_GetItem(args,1)
+#define AREA          PyTuple_GetItem(args,2)
+/* serial output macros: */
+#define INDEX             output,0
+#define X                 output,1
+#define Y                 output,2
+#define SEGMENTS          output,3
+#define SEGMENTMARKERLIST output,4
+#endif
+
+#ifdef _HAVE_JAVASCRIPT_MODULES_
+/* serial input macros: */
+#define DOMAINOUTLINE domainx,domainy,domainnods
+#define RIFTSOUTLINE  NULL,NULL,0
+#define AREA          areain
+/* serial output macros: */
+#define INDEX             pindex,pnel
+#define X                 px,pnods
+#define Y                 py,pnods
+#define SEGMENTS          psegments,pnsegs
+#define SEGMENTMARKERLIST psegmentmarkers,pnsegs
+#define WRAPPER(modulename) extern "C" { int  TriangleModule(double** pindex, double** px, double** py, int* pnel, int* pnods, double** psegments, double** psegmentmarkers, int* pnsegs, double* domainx, double* domainy, int domainnods, double areain)
+#define _DO_NOT_LOAD_GLOBALS_//we only load globals for TriangleModule.js, not other modules!
+#endif
+
+
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  5
+#undef NRHS
+#define NRHS  3
+
+#endif  /* _TRIANGLE_H */
Index: /issm/trunk-jpl/src/wrappers/Triangle/Triangle.js
===================================================================
--- /issm/trunk-jpl/src/wrappers/Triangle/Triangle.js	(revision 22498)
+++ /issm/trunk-jpl/src/wrappers/Triangle/Triangle.js	(revision 22498)
@@ -0,0 +1,82 @@
+function Triangle(md,domain,rifts, area){
+/*Triangle 
+	   usage: var array = Triangle(domain,rifts,area);
+	      where: array is made of [index,x,y,segments,segmentmarkers]
+		  and index,x,y defines a triangulation, segments is an array made 
+	      of exterior segments to the mesh domain outline, segmentmarkers is an array 
+		  flagging each segment, domain a js array defining the domain outline  (sames for 
+		  rifts) and area is the maximum area desired for any element of the resulting mesh.
+
+		  Ok, for now, we are not dealing with rifts. Also, the domain is made of only one 
+		  profile, this to avoid passing a double** pointer to js. 
+*/
+
+	//Dynamic allocations: {{{
+	//Retrieve domain arrays, and allocate on Module heap: 
+	
+	//input
+	var dx=new Float64Array(domain['x']); var nx=dx.length * dx.BYTES_PER_ELEMENT;
+	var dxPtr= Module._malloc(nx); var domainxHeap = new Uint8Array(Module.HEAPU8.buffer,dxPtr,nx);
+	domainxHeap.set(new Uint8Array(dx.buffer)); var domainx=domainxHeap.byteOffset;
+
+	var dy=new Float64Array(domain['y']); var ny=dy.length * dy.BYTES_PER_ELEMENT;
+	var dyPtr = Module._malloc(ny); var domainyHeap = new Uint8Array(Module.HEAPU8.buffer,dyPtr,ny);
+	domainyHeap.set(new Uint8Array(dy.buffer)); var domainy=domainyHeap.byteOffset;
+	
+	//output
+	var nel,indexlinear,index,nods,x,y;
+	var pnel= Module._malloc(4); 
+	var pindex= Module._malloc(4); 
+	var pnods= Module._malloc(4); 
+	var px= Module._malloc(4); 
+	var py= Module._malloc(4); 
+	var psegments= Module._malloc(4); 
+	var psegmentmarkers= Module._malloc(4); 
+	var pnsegs= Module._malloc(4); 
+	//}}}
+
+	//Declare Triangle module: 
+	TriangleModule = Module.cwrap('TriangleModule','number',['number','number','number','number','number','number','number','number','number','number','number','number']);
+	
+	//Call Triangle module: 
+	TriangleModule(pindex,px,py,pnel,pnods,psegments,psegmentmarkers,pnsegs, domainx,domainy,dx.length,area);
+	
+	/*Dynamic copying from heap: {{{*/
+	//recover mesh: 
+	nel = Module.getValue(pnel, 'i32');
+	var indexptr = Module.getValue(pindex,'i32');
+	indexlinear = Module.HEAPF64.slice(indexptr /8, indexptr/8 + nel*3);
+	index = ListToMatrix(indexlinear,3);
+
+	nods = Module.getValue(pnods, 'i32');
+	var xptr = Module.getValue(px,'i32');
+	var yptr = Module.getValue(py,'i32');
+	x = Module.HEAPF64.slice(xptr /8, xptr/8 + nods);
+	y = Module.HEAPF64.slice(yptr /8, yptr/8 + nods);
+	
+	nsegs = Module.getValue(pnsegs, 'i32');
+	var segmentsptr = Module.getValue(psegments,'i32');
+	segmentslinear = Module.HEAPF64.slice(segmentsptr /8, segmentsptr/8 + nsegs*3);
+	segments = ListToMatrix(segmentslinear,3);
+	
+	var segmentmarkersptr = Module.getValue(psegmentmarkers,'i32');
+	segmentmarkers = Module.HEAPF64.slice(segmentmarkersptr /8, segmentmarkersptr/8 + nsegs);
+	/*}}}*/
+
+	var return_array=[index,x,y,segments,segmentmarkers];
+
+	/*Free ressources: */
+	Module._free(pindex); 
+	Module._free(indexlinear); 
+	Module._free(px); 
+	Module._free(x); 
+	Module._free(py); 
+	Module._free(y); 
+	Module._free(pnel); 
+	Module._free(pnods); 
+	Module._free(psegments); 
+	Module._free(psegmentmarkers); 
+	Module._free(pnsegs); 
+
+	return return_array;
+}
Index: /issm/trunk-jpl/src/wrappers/javascript/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/wrappers/javascript/Makefile.am	(revision 22497)
+++ /issm/trunk-jpl/src/wrappers/javascript/Makefile.am	(revision 22498)
@@ -7,5 +7,5 @@
 AM_CPPFLAGS+=  -DISSM_PREFIX='"$(prefix)"'
 
-js_scripts = ${ISSM_DIR}/src/wrappers/TriMesh/TriMesh.js  \
+js_scripts = ${ISSM_DIR}/src/wrappers/Triangle/Triangle.js  \
 			 ${ISSM_DIR}/src/wrappers/NodeConnectivity/NodeConnectivity.js\
 			 ${ISSM_DIR}/src/wrappers/ContourToMesh/ContourToMesh.js\
@@ -81,5 +81,5 @@
 endif
 
-IssmModule_SOURCES = ../TriMesh/TriMesh.cpp \
+IssmModule_SOURCES = ../Triangle/Triangle.cpp \
 					 ../NodeConnectivity/NodeConnectivity.cpp\
 					 ../ContourToMesh/ContourToMesh.cpp\
@@ -89,5 +89,5 @@
 					 ../Issm/issm.cpp
 
-IssmModule_CXXFLAGS= -fPIC -D_DO_NOT_LOAD_GLOBALS_  --memory-init-file 0 $(AM_CXXFLAGS) $(CXXFLAGS) $(CXXOPTFLAGS) $(COPTFLAGS) -s EXPORTED_FUNCTIONS="['_TriMeshModule','_NodeConnectivityModule','_ContourToMeshModule','_ElementConnectivityModule','_InterpFromMeshToMesh2dModule','_IssmConfigModule','_IssmModule']"  -s DISABLE_EXCEPTION_CATCHING=0 -s ALLOW_MEMORY_GROWTH=1 -s INVOKE_RUN=0
+IssmModule_CXXFLAGS= -fPIC -D_DO_NOT_LOAD_GLOBALS_  --memory-init-file 0 $(AM_CXXFLAGS) $(CXXFLAGS) $(CXXOPTFLAGS) $(COPTFLAGS) -s EXPORTED_FUNCTIONS="['_TriangleModule','_NodeConnectivityModule','_ContourToMeshModule','_ElementConnectivityModule','_InterpFromMeshToMesh2dModule','_IssmConfigModule','_IssmModule']"  -s DISABLE_EXCEPTION_CATCHING=0 -s ALLOW_MEMORY_GROWTH=1 -s INVOKE_RUN=0
 IssmModule_LDADD = ${deps} $(TRIANGLELIB)  $(GSLLIB)
 #}}}
Index: /issm/trunk-jpl/src/wrappers/matlab/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/wrappers/matlab/Makefile.am	(revision 22497)
+++ /issm/trunk-jpl/src/wrappers/matlab/Makefile.am	(revision 22498)
@@ -58,6 +58,6 @@
 						 PointCloudFindNeighbors_matlab.la\
 						 PropagateFlagsFromConnectivity_matlab.la\
-						 TriMesh_matlab.la\
-						 TriMeshProcessRifts_matlab.la\
+						 Triangle_matlab.la\
+						 ProcessRifts_matlab.la\
 						 Scotch_matlab.la
 
@@ -233,10 +233,10 @@
 ShpRead_matlab_la_LIBADD = ${deps} $(SHAPELIBLIB) $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
 
-TriMesh_matlab_la_SOURCES = ../TriMesh/TriMesh.cpp
-TriMesh_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
-TriMesh_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(TRIANGLELIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
-
-TriMeshProcessRifts_matlab_la_SOURCES = ../TriMeshProcessRifts/TriMeshProcessRifts.cpp
-TriMeshProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
-TriMeshProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
-#}}}
+Triangle_matlab_la_SOURCES = ../Triangle/Triangle.cpp
+Triangle_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
+Triangle_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(TRIANGLELIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
+
+ProcessRifts_matlab_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp
+ProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
+ProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
+#}}}
Index: /issm/trunk-jpl/src/wrappers/python/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/wrappers/python/Makefile.am	(revision 22497)
+++ /issm/trunk-jpl/src/wrappers/python/Makefile.am	(revision 22498)
@@ -41,6 +41,6 @@
 						MeshProfileIntersection_python.la\
 						NodeConnectivity_python.la\
-						TriMesh_python.la\
-						TriMeshProcessRifts_python.la
+						Triangle_python.la\
+						ProcessRifts_python.la
 endif 
 #}}}
@@ -141,10 +141,10 @@
 NodeConnectivity_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
 
-TriMesh_python_la_SOURCES = ../TriMesh/TriMesh.cpp
-TriMesh_python_la_CXXFLAGS = ${AM_CXXFLAGS}
-TriMesh_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)
+Triangle_python_la_SOURCES = ../Triangle/Triangle.cpp
+Triangle_python_la_CXXFLAGS = ${AM_CXXFLAGS}
+Triangle_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)
 
-TriMeshProcessRifts_python_la_SOURCES = ../TriMeshProcessRifts/TriMeshProcessRifts.cpp
-TriMeshProcessRifts_python_la_CXXFLAGS = ${AM_CXXFLAGS}
-TriMeshProcessRifts_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
+ProcessRifts_python_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp
+ProcessRifts_python_la_CXXFLAGS = ${AM_CXXFLAGS}
+ProcessRifts_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
 #}}}