Index: /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp	(revision 11880)
+++ /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp	(revision 11880)
@@ -0,0 +1,169 @@
+/*!\file TriMeshx
+ * \brief: x code for TriMesh mesher
+ */
+
+/*Header files: {{{*/
+#include "./TriMeshx.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+/*}}}*/
+
+
+void TriMeshx(Matrix** pindex,Vector** px,Vector** py,Matrix** psegments,Vector** psegmentmarkerlist,DataSet* domain,double area,bool order){
+
+	/*indexing: */
+	int i,j;
+
+	/*output: */
+	double* index=NULL;
+	double* x=NULL;
+	double* y=NULL;
+	double* segments=NULL;
+	double* segmentmarkerlist=NULL;
+
+	/*intermediary: */
+	int      counter,backcounter;
+	Contour* 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*)domain->GetObjectByOffset(i);
+		in.numberofpoints+=contour->nods;
+	}
+	/*number of point attributes: */
+	in.numberofpointattributes=1;
+
+	/*fill in the point list: */
+	in.pointlist = (REAL *) xmalloc(in.numberofpoints * 2 * sizeof(REAL));
+
+	counter=0;
+	for (i=0;i<domain->Size();i++){
+		contour=(Contour*)domain->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 = (REAL *) xmalloc(in.numberofpoints * in.numberofpointattributes * sizeof(REAL));
+	for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0;
+	
+	/*fill in the point marker list: */
+	in.pointmarkerlist = (int *) xmalloc(in.numberofpoints * sizeof(int));
+	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*)domain->GetObjectByOffset(i);
+		in.numberofsegments+=contour->nods;
+	}
+	
+	in.segmentlist = (int *) xmalloc(in.numberofsegments * 2 * sizeof(int));
+	in.segmentmarkerlist = (int *) xcalloc(in.numberofsegments,sizeof(int));
+	counter=0;
+	backcounter=0;
+	for (i=0;i<domain->Size();i++){
+		contour=(Contour*)domain->GetObjectByOffset(i);
+		for (j=0;j<contour->nods-1;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;
+	}
+	
+	/*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 = (REAL *) xmalloc(in.numberofholes * 2 * sizeof(REAL));
+		for (i=0;i<domain->Size()-1;i++){
+			contour=(Contour*)domain->GetObjectByOffset(i+1);
+			GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods,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=(double*)xmalloc(3*out.numberoftriangles*sizeof(double));
+	x=(double*)xmalloc(out.numberofpoints*sizeof(double));
+	y=(double*)xmalloc(out.numberofpoints*sizeof(double));
+	segments=(double*)xmalloc(3*out.numberofsegments*sizeof(double));
+	segmentmarkerlist=(double*)xmalloc(out.numberofsegments*sizeof(double));
+
+	for (i = 0; i < out.numberoftriangles; i++) {
+		for (j = 0; j < out.numberofcorners; j++) {
+			*(index+3*i+j)=(double)out.trianglelist[i * out.numberofcorners + j]+1;
+		}
+	}
+	for (i = 0; i < out.numberofpoints; i++) {
+		x[i]=out.pointlist[i * 2 + 0];
+		y[i]=out.pointlist[i * 2 + 1];
+	}
+	
+	for (i = 0; i < out.numberofsegments; i++) {
+		segments[3*i+0]=(double)out.segmentlist[i*2+0]+1;
+		segments[3*i+1]=(double)out.segmentlist[i*2+1]+1;
+		segmentmarkerlist[i]=(double)out.segmentmarkerlist[i];
+	}
+
+	/*Associate elements with segments: */
+	AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
+
+	/*Order segments so that their normals point outside the domain: */
+	if(order){
+		OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
+	}
+
+	printf("%i %i %i\n",out.numberoftriangles,out.numberofpoints,out.numberofsegments);
+
+	/*Output : */
+	*pindex=new Matrix(index,out.numberoftriangles,3,1);
+	*px=new Vector(x,out.numberofpoints);
+	*py=new Vector(y,out.numberofpoints);
+	*psegments=new Matrix(segments,out.numberofsegments,3,1);
+	*psegmentmarkerlist=new Vector(segmentmarkerlist,out.numberofsegments);
+}
Index: /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp.bak
===================================================================
--- /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp.bak	(revision 11880)
+++ /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp.bak	(revision 11880)
@@ -0,0 +1,260 @@
+/*
+ * TriMesh: mesh a domain using an .exp file
+ */
+
+#include "./TriMesh.h"
+
+void mexFunction(	int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] ){
+
+	/* output: */
+	double* index=NULL;
+	double* x=NULL;
+	double* y=NULL;
+	double* segments=NULL;
+	double* segmentmarkerlist=NULL;
+
+	/* input: */
+	char*  domainname=NULL;
+	double area;
+	char*  order=NULL;
+	
+	/*Domain outline variables: */
+	int      nprof;
+	int*     profnvertices=NULL;
+	double** pprofx=NULL;
+	double** pprofy=NULL;
+	double*  xprof=NULL;
+	double*  yprof=NULL;
+	int      numberofpoints;
+
+	/* Triangle structures: */
+	struct triangulateio in,out;
+	char   options[256];
+	
+	/*Boot module: */
+	MODULEBOOT();
+	
+	/*checks on arguments on the matlab side: */
+	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&TriMeshUsage);
+
+	/*Fetch data needed for meshing: */
+	FetchMatlabData(&domainname,DOMAINOUTLINE);
+	FetchMatlabData(&area,AREA);
+	FetchMatlabData(&order,ORDER);
+
+	/*Read domain outline: */
+	domain=DomainOutlineRead(domainname,false);
+
+	/*Create initial triangulation to call triangulate():*/
+	numberofpoints=0;
+	for (i=0;i<nprof;i++){
+		numberofpoints+=profnvertices[i];
+	}
+	if (riftname){
+		for (i=0;i<numrifts;i++){
+			numberofpoints+=riftsnumvertices[i];
+		}
+	}
+	in.numberofpoints=numberofpoints;
+
+	in.numberofpointattributes=1;
+	in.pointlist = (REAL *) mxMalloc(in.numberofpoints * 2 * sizeof(REAL));
+
+	counter=0;
+	for (i=0;i<nprof;i++){
+		xprof=pprofx[i];
+		yprof=pprofy[i];
+		for (j=0;j<profnvertices[i];j++){
+			in.pointlist[2*counter+0]=xprof[j];
+			in.pointlist[2*counter+1]=yprof[j];
+			counter++;
+		}
+	}
+	if(riftname){
+		for (i=0;i<numrifts;i++){
+			xprof=riftsverticesx[i];
+			yprof=riftsverticesy[i];
+			for (j=0;j<riftsnumvertices[i];j++){
+				in.pointlist[2*counter+0]=xprof[j];
+				in.pointlist[2*counter+1]=yprof[j];
+				counter++;
+			}
+		}
+	}
+	
+	in.pointattributelist = (REAL *) mxMalloc(in.numberofpoints *
+										  in.numberofpointattributes *
+										  sizeof(REAL));
+	for (i=0;i<in.numberofpoints;i++){
+		in.pointattributelist[i] = 0.0;
+	}
+	in.pointmarkerlist = (int *) mxMalloc(in.numberofpoints * sizeof(int));
+	for(i=0;i<in.numberofpoints;i++){
+		in.pointmarkerlist[i] = 0;
+	}
+
+	/*Build segments: */
+	/*Figure out number of segments: holes and closed outlines have as many segments as vertices, 
+	 *for rifts, we have one less segment as we have vertices*/
+	in.numberofsegments=0;
+	for (i=0;i<nprof;i++){
+		in.numberofsegments+=profnvertices[i];
+	}
+	if (riftname){
+		for (i=0;i<numrifts;i++){
+			in.numberofsegments+=riftsnumvertices[i]-1;
+		}
+	}
+	
+	in.segmentlist = (int *) mxMalloc(in.numberofsegments * 2 * sizeof(int));
+	in.segmentmarkerlist = (int *) mxCalloc(in.numberofsegments,sizeof(int));
+	counter=0;
+	backcounter=0;
+	for (i=0;i<nprof;i++){
+		for (j=0;j<(profnvertices[i]-1);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;
+	if(riftname){
+		for (i=0;i<numrifts;i++){
+			for (j=0;j<(riftsnumvertices[i]-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 = nprof-1; /*everything is a hole, but for the first profile.*/
+	in.holelist = (REAL *) mxMalloc(in.numberofholes * 2 * sizeof(REAL));
+	for (i=0;i<nprof-1;i++){
+		/*We are looking for a vertex that lies inside the hole: */
+		GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],profnvertices[i+1],pprofx[i+1],pprofy[i+1]);
+	}
+
+	/* 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=(double*)mxMalloc(3*out.numberoftriangles*sizeof(double));
+	x=(double*)mxMalloc(out.numberofpoints*sizeof(double));
+	y=(double*)mxMalloc(out.numberofpoints*sizeof(double));
+	segments=(double*)mxMalloc(3*out.numberofsegments*sizeof(double));
+	segmentmarkerlist=(double*)mxMalloc(out.numberofsegments*sizeof(double));
+
+	for (i = 0; i < out.numberoftriangles; i++) {
+		for (j = 0; j < out.numberofcorners; j++) {
+			*(index+3*i+j)=(double)out.trianglelist[i * out.numberofcorners + j]+1;
+		}
+	}
+	for (i = 0; i < out.numberofpoints; i++) {
+		x[i]=out.pointlist[i * 2 + 0];
+		y[i]=out.pointlist[i * 2 + 1];
+	}
+	
+	for (i = 0; i < out.numberofsegments; i++) {
+		segments[3*i+0]=(double)out.segmentlist[i*2+0]+1;
+		segments[3*i+1]=(double)out.segmentlist[i*2+1]+1;
+		segmentmarkerlist[i]=(double)out.segmentmarkerlist[i];
+	}
+
+	/*Associate elements with segments: */
+	AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
+
+	/*Order segments so that their normals point outside the domain: */
+	if(!strcmp(order,"yes")){
+		OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
+	}
+
+	/*Output : */
+	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(pmxa_array,3);
+	mxSetN(pmxa_array,out.numberoftriangles);
+	mxSetPr(pmxa_array,index);
+	mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "transpose");
+
+	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(pmxa_array,1);
+	mxSetN(pmxa_array,out.numberofpoints);
+	mxSetPr(pmxa_array,x);
+	mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "transpose");
+
+	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(pmxa_array,1);
+	mxSetN(pmxa_array,out.numberofpoints);
+	mxSetPr(pmxa_array,y);
+	mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "transpose");
+
+	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(pmxa_array,3);
+	mxSetN(pmxa_array,out.numberofsegments);
+	mxSetPr(pmxa_array,segments);
+	mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "transpose");
+
+	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+	mxSetM(pmxa_array,1);
+	mxSetN(pmxa_array,out.numberofsegments);
+	mxSetPr(pmxa_array,segmentmarkerlist);
+	mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "transpose");
+	
+	return;
+}
+
+void TriMeshUsage(void)
+{
+	printf("\n");
+	printf("   usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,riftsoutlinename,area,ordered) \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("      (if rifts are present, markers >=2 flag them ), outlinefilename an Argus domain outline file.\n");
+	printf("      riftsoutlinename is an Argus domain file, defining rifts (ie: open profiles), \n");
+	printf("      area is the maximum area desired for any element of the resulting mesh. \n");
+	printf("      and ordered is a string ('yes' or 'no') that determines whether segments are output in the \n");
+	printf("      order they are made by Triangle (ie none), or ordered counter clockwise around the domain outline.\n");
+	printf("      riftsoutlinename and ordered are optional arguments.\n");
+	printf("\n");
+}
+
+
+
+
Index: /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h	(revision 11880)
+++ /issm/trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h	(revision 11880)
@@ -0,0 +1,27 @@
+/*!\file:  TriMeshx.h
+ * \brief header file for TriMeshx module
+ */ 
+
+#ifndef _TRIMESHX_H_
+#define _TRIMESHX_H_
+
+
+
+/*ANSI_DECLARATORS needed to call triangle library: */
+#ifndef ANSI_DECLARATORS
+#define ANSI_DECLARATORS
+#include "triangle.h"
+#undef ANSI_DECLARATORS
+#else
+#include "triangle.h"
+#endif
+
+#include "string.h"
+
+#include "../../Container/Container.h"
+#include "../../objects/objects.h"
+
+/* local prototypes: */
+void TriMeshx(Matrix** pindex,Vector** px,Vector** py,Matrix** psegments,Vector** psegmentmarkerlist,DataSet* domain,double area,bool order);
+
+#endif  /* _TRIMESHX_H */