Index: ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp (revision 22497) +++ ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp (nonexistent) @@ -1,201 +0,0 @@ -/*!\file TriMeshx - * \brief: x code for TriMesh mesher - */ - -/*Header files: {{{*/ -#include "./TriMeshx.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 TriMeshx(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 *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;iSize();i++){ - contour=(Contour*)domain->GetObjectByOffset(i); - in.numberofpoints+=contour->nods-1; - } - for (i=0;iSize();i++){ - contour=(Contour*)rifts->GetObjectByOffset(i); - in.numberofpoints+=contour->nods; - } - - /*number of point attributes: */ - in.numberofpointattributes=1; - - /*fill in the point list: */ - in.pointlist = xNew(in.numberofpoints*2); - - counter=0; - for (i=0;iSize();i++){ - contour=(Contour*)domain->GetObjectByOffset(i); - for (j=0;jnods-1;j++){ - in.pointlist[2*counter+0]=contour->x[j]; - in.pointlist[2*counter+1]=contour->y[j]; - counter++; - } - } - for (i=0;iSize();i++){ - contour=(Contour*)rifts->GetObjectByOffset(i); - for (j=0;jnods;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(in.numberofpoints*in.numberofpointattributes); - for (i=0;i(in.numberofpoints); - for(i=0;iSize();i++){ - contour=(Contour*)domain->GetObjectByOffset(i); - in.numberofsegments+=contour->nods-1; - } - for(i=0;iSize();i++){ - contour=(Contour*)rifts->GetObjectByOffset(i); - /*for rifts, we have one less segment as we have vertices*/ - in.numberofsegments+=contour->nods-1; - } - - in.segmentlist = xNew(in.numberofsegments*2); - in.segmentmarkerlist = xNewZeroInit(in.numberofsegments); - counter=0; - backcounter=0; - for (i=0;iSize();i++){ - contour=(Contour*)domain->GetObjectByOffset(i); - for (j=0;jnods-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;iSize();i++){ - contour=(Contour*)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(in.numberofholes*2); - for (i=0;iSize()-1;i++){ - contour=(Contour*)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(3*out.numberoftriangles); - x=xNew(out.numberofpoints); - y=xNew(out.numberofpoints); - segments=xNew(3*out.numberofsegments); - segmentmarkerlist=xNew(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 -#include "../../classes/classes.h" - -/* local prototypes: */ -void TriMeshx(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnseg,Contours* domain,Contours* rifts,double area); -#endif /* _TRIMESHX_H */ Index: ../trunk-jpl/src/c/modules/TriMeshx =================================================================== --- ../trunk-jpl/src/c/modules/TriMeshx (revision 22497) +++ ../trunk-jpl/src/c/modules/TriMeshx (nonexistent) Property changes on: ../trunk-jpl/src/c/modules/TriMeshx ___________________________________________________________________ Deleted: svn:ignore ## -1,2 +0,0 ## -.deps -.dirstamp Index: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h =================================================================== --- ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h (revision 22497) +++ ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h (nonexistent) @@ -1,12 +0,0 @@ -/*!\file: TriMeshProcessRiftsx.h - * \brief header file for TriMeshProcessRifts module - */ - -#ifndef _TRIMESHPROCESSRIFTX_H -#define _TRIMESHPROCESSRIFTX_H - -class RiftStruct; - -void TriMeshProcessRiftsx(int** pindex,int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct); - -#endif /* _TRIMESHPROCESSRIFTX_H*/ Index: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp (revision 22497) +++ ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp (nonexistent) @@ -1,78 +0,0 @@ -/*!\file: TriMeshProcessRifts.cpp - * \brief split a mesh where a rift (or fault) is present - */ - -#include "./TriMeshProcessRiftsx.h" -#include "../../classes/RiftStruct.h" -#include "../../shared/shared.h" -#include "../../toolkits/toolkits.h" - -void TriMeshProcessRiftsx(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: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx =================================================================== --- ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx (revision 22497) +++ ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx (nonexistent) Property changes on: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx ___________________________________________________________________ Deleted: svn:ignore ## -1,2 +0,0 ## -.deps -.dirstamp Index: ../trunk-jpl/src/c/modules/modules.h =================================================================== --- ../trunk-jpl/src/c/modules/modules.h (revision 22497) +++ ../trunk-jpl/src/c/modules/modules.h (revision 22498) @@ -90,8 +90,8 @@ #include "./SystemMatricesx/SystemMatricesx.h" #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" #include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h" Index: ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h =================================================================== --- ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h (nonexistent) +++ ../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: ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp (nonexistent) +++ ../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: ../trunk-jpl/src/c/modules/ProcessRiftsx =================================================================== --- ../trunk-jpl/src/c/modules/ProcessRiftsx (nonexistent) +++ ../trunk-jpl/src/c/modules/ProcessRiftsx (revision 22498) Property changes on: ../trunk-jpl/src/c/modules/ProcessRiftsx ___________________________________________________________________ Added: svn:ignore ## -0,0 +1,2 ## +.deps +.dirstamp Index: ../trunk-jpl/src/c/modules/Trianglex/Trianglex.h =================================================================== --- ../trunk-jpl/src/c/modules/Trianglex/Trianglex.h (nonexistent) +++ ../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 +#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: ../trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp (nonexistent) +++ ../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 *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;iSize();i++){ + contour=(Contour*)domain->GetObjectByOffset(i); + in.numberofpoints+=contour->nods-1; + } + for (i=0;iSize();i++){ + contour=(Contour*)rifts->GetObjectByOffset(i); + in.numberofpoints+=contour->nods; + } + + /*number of point attributes: */ + in.numberofpointattributes=1; + + /*fill in the point list: */ + in.pointlist = xNew(in.numberofpoints*2); + + counter=0; + for (i=0;iSize();i++){ + contour=(Contour*)domain->GetObjectByOffset(i); + for (j=0;jnods-1;j++){ + in.pointlist[2*counter+0]=contour->x[j]; + in.pointlist[2*counter+1]=contour->y[j]; + counter++; + } + } + for (i=0;iSize();i++){ + contour=(Contour*)rifts->GetObjectByOffset(i); + for (j=0;jnods;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(in.numberofpoints*in.numberofpointattributes); + for (i=0;i(in.numberofpoints); + for(i=0;iSize();i++){ + contour=(Contour*)domain->GetObjectByOffset(i); + in.numberofsegments+=contour->nods-1; + } + for(i=0;iSize();i++){ + contour=(Contour*)rifts->GetObjectByOffset(i); + /*for rifts, we have one less segment as we have vertices*/ + in.numberofsegments+=contour->nods-1; + } + + in.segmentlist = xNew(in.numberofsegments*2); + in.segmentmarkerlist = xNewZeroInit(in.numberofsegments); + counter=0; + backcounter=0; + for (i=0;iSize();i++){ + contour=(Contour*)domain->GetObjectByOffset(i); + for (j=0;jnods-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;iSize();i++){ + contour=(Contour*)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(in.numberofholes*2); + for (i=0;iSize()-1;i++){ + contour=(Contour*)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(3*out.numberoftriangles); + x=xNew(out.numberofpoints); + y=xNew(out.numberofpoints); + segments=xNew(3*out.numberofsegments); + segmentmarkerlist=xNew(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(nods); //to make sure we don't split the same nodes twice! - for (i=0;i(x,nods,nods+1); - y=xReNew(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 - -#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: ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp =================================================================== --- ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp (revision 22497) +++ ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp (nonexistent) @@ -1,24 +0,0 @@ -/*!\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 - -#include "./trimesh.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(max_number_elements); - - for (i=0;i(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(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(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(nriftsegs*4); - counter=0; - for (i=0;i(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(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(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(segments, nsegs*3,(nsegs+nriftsegs)*3); - segmentmarkerlist=xReNew(segmentmarkerlist,nsegs,(nsegs+nriftsegs)); - - /*First, update the existing segments to the new nodes :*/ - for (i=0;i(new_numsegs*3); - new_segmentmarkers=xNew(new_numsegs); - - /*Copy new segments info : */ - counter=0; - for (i=0;i(numrifts); - riftssegments=xNew(numrifts); - for (i=0;i(counter*3); - /*Copy new segments info :*/ - counter=0; - for (j=0;j(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(numrifts); - riftspairs=xNew(numrifts); - for (i=0;i(2*numsegs); - for (j=0;j=2 indicates a certain rift: */ - numrifts=0; - for (i=0;imaxmark){ - 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(numrifts*2); - - /*Go through all rifts: */ - for (i=0;i(numsegs*3); - riftpairs_copy=xNew(numsegs*2); - order=xNew(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;j0 && 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(order); - xDelete(riftsegments_copy); - xDelete(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(numrifts); - riftsnumpenaltypairs=xNew(numrifts); - - for(i=0;i((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(x,nods,nods+1); - y=xReNew(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(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 -#include - -//#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_TRIMESH_H */ Index: ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp =================================================================== --- ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp (revision 22497) +++ ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp (nonexistent) @@ -1,46 +0,0 @@ -/* - * 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 + +#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(max_number_elements); + + for (i=0;i(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(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(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(nriftsegs*4); + counter=0; + for (i=0;i(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(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(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(segments, nsegs*3,(nsegs+nriftsegs)*3); + segmentmarkerlist=xReNew(segmentmarkerlist,nsegs,(nsegs+nriftsegs)); + + /*First, update the existing segments to the new nodes :*/ + for (i=0;i(new_numsegs*3); + new_segmentmarkers=xNew(new_numsegs); + + /*Copy new segments info : */ + counter=0; + for (i=0;i(numrifts); + riftssegments=xNew(numrifts); + for (i=0;i(counter*3); + /*Copy new segments info :*/ + counter=0; + for (j=0;j(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(numrifts); + riftspairs=xNew(numrifts); + for (i=0;i(2*numsegs); + for (j=0;j=2 indicates a certain rift: */ + numrifts=0; + for (i=0;imaxmark){ + 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(numrifts*2); + + /*Go through all rifts: */ + for (i=0;i(numsegs*3); + riftpairs_copy=xNew(numsegs*2); + order=xNew(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;j0 && 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(order); + xDelete(riftsegments_copy); + xDelete(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(numrifts); + riftsnumpenaltypairs=xNew(numrifts); + + for(i=0;i((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(x,nods,nods+1); + y=xReNew(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(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(nods); //to make sure we don't split the same nodes twice! + for (i=0;i(x,nods,nods+1); + y=xReNew(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 + +#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: ../trunk-jpl/src/c/shared/Triangle/triangle.h =================================================================== --- ../trunk-jpl/src/c/shared/Triangle/triangle.h (nonexistent) +++ ../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 +#include + +//#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: ../trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp (nonexistent) +++ ../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 -#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__ "TriMeshProcessRifts" - -#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: ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp =================================================================== --- ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 22497) +++ ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp (nonexistent) @@ -1,59 +0,0 @@ -/*!\file: TriMeshProcessRifts.cpp - * \brief split a mesh where a rift (or fault) is present - */ - -#include "./TriMeshProcessRifts.h" - -void TriMeshProcessRiftsUsage(void){/*{{{*/ - _printf_("\n"); - _printf_(" usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(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(TriMeshProcessRifts_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,&TriMeshProcessRiftsUsage); - - /*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*/ - TriMeshProcessRiftsx(&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(index); - xDelete(x); - xDelete(y); - xDelete(segments); - xDelete(segmentmarkers ); - MODULEEND(); -} Index: ../trunk-jpl/src/wrappers/TriMeshProcessRifts =================================================================== --- ../trunk-jpl/src/wrappers/TriMeshProcessRifts (revision 22497) +++ ../trunk-jpl/src/wrappers/TriMeshProcessRifts (nonexistent) Property changes on: ../trunk-jpl/src/wrappers/TriMeshProcessRifts ___________________________________________________________________ Deleted: svn:ignore ## -1,2 +0,0 ## -.deps -.dirstamp Index: ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp =================================================================== --- ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp (revision 22497) +++ ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp (nonexistent) @@ -1,63 +0,0 @@ -/* - * TriMesh: mesh a domain using an .exp file - */ - -#include "./TriMesh.h" - -void TriMeshUsage(void){/*{{{*/ - _printf_("\n"); - _printf_(" usage: [index,x,y,segments,segmentmarkers]=TriMesh(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(TriMesh_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,&TriMeshUsage); - - /*Fetch data needed for meshing: */ - FetchData(&domain,DOMAINOUTLINE); - FetchData(&rifts,RIFTSOUTLINE); - FetchData(&area,AREA); - - /*call x core: */ - TriMeshx(&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(index); - xDelete(x); - xDelete(y); - xDelete(segments); - xDelete(segmentmarkerlist); - - /*end module: */ - MODULEEND(); -} Index: ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h =================================================================== --- ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h (revision 22497) +++ ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h (nonexistent) @@ -1,83 +0,0 @@ -/* - TriMesh.h -*/ - -#ifndef _TRIMESH_H -#define _TRIMESH_H - -#ifdef HAVE_CONFIG_H - #include -#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__ "TriMesh" - -#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 TriMeshModule(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 TriMeshModule.js, not other modules! -#endif - - -/* serial arg counts: */ -#undef NLHS -#define NLHS 5 -#undef NRHS -#define NRHS 3 - -#endif /* _TRIMESH_H */ Index: ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js =================================================================== --- ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js (revision 22497) +++ ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js (nonexistent) @@ -1,82 +0,0 @@ -function TriMesh(md,domain,rifts, area){ -/*TriMesh - usage: var array = TriMesh(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 TriMesh module: - TriMeshModule = Module.cwrap('TriMeshModule','number',['number','number','number','number','number','number','number','number','number','number','number','number']); - - //Call TriMesh module: - TriMeshModule(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: ../trunk-jpl/src/wrappers/TriMesh =================================================================== --- ../trunk-jpl/src/wrappers/TriMesh (revision 22497) +++ ../trunk-jpl/src/wrappers/TriMesh (nonexistent) Property changes on: ../trunk-jpl/src/wrappers/TriMesh ___________________________________________________________________ Deleted: svn:ignore ## -1,2 +0,0 ## -.deps -.dirstamp Index: ../trunk-jpl/src/wrappers/Triangle/Triangle.js =================================================================== --- ../trunk-jpl/src/wrappers/Triangle/Triangle.js (nonexistent) +++ ../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: ../trunk-jpl/src/wrappers/Triangle/Triangle.cpp =================================================================== --- ../trunk-jpl/src/wrappers/Triangle/Triangle.cpp (nonexistent) +++ ../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(index); + xDelete(x); + xDelete(y); + xDelete(segments); + xDelete(segmentmarkerlist); + + /*end module: */ + MODULEEND(); +} Index: ../trunk-jpl/src/wrappers/Triangle/Triangle.h =================================================================== --- ../trunk-jpl/src/wrappers/Triangle/Triangle.h (nonexistent) +++ ../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 +#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: ../trunk-jpl/src/wrappers/Triangle =================================================================== --- ../trunk-jpl/src/wrappers/Triangle (nonexistent) +++ ../trunk-jpl/src/wrappers/Triangle (revision 22498) Property changes on: ../trunk-jpl/src/wrappers/Triangle ___________________________________________________________________ Added: svn:ignore ## -0,0 +1,2 ## +.deps +.dirstamp Index: ../trunk-jpl/src/wrappers/javascript/Makefile.am =================================================================== --- ../trunk-jpl/src/wrappers/javascript/Makefile.am (revision 22497) +++ ../trunk-jpl/src/wrappers/javascript/Makefile.am (revision 22498) @@ -6,7 +6,7 @@ #define prefix (from http://www.gnu.org/software/autoconf/manual/autoconf-2.67/html_node/Defining-Directories.html) 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\ ${ISSM_DIR}/src/wrappers/ElementConnectivity/ElementConnectivity.js\ @@ -80,7 +80,7 @@ libISSMApi_la_LDFLAGS = -static endif -IssmModule_SOURCES = ../TriMesh/TriMesh.cpp \ +IssmModule_SOURCES = ../Triangle/Triangle.cpp \ ../NodeConnectivity/NodeConnectivity.cpp\ ../ContourToMesh/ContourToMesh.cpp\ ../ElementConnectivity/ElementConnectivity.cpp\ @@ -88,6 +88,6 @@ ../IssmConfig/IssmConfig.cpp\ ../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: ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp =================================================================== --- ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp (nonexistent) +++ ../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(index); + xDelete(x); + xDelete(y); + xDelete(segments); + xDelete(segmentmarkers ); + MODULEEND(); +} Index: ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h =================================================================== --- ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h (nonexistent) +++ ../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 +#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: ../trunk-jpl/src/wrappers/ProcessRifts =================================================================== --- ../trunk-jpl/src/wrappers/ProcessRifts (nonexistent) +++ ../trunk-jpl/src/wrappers/ProcessRifts (revision 22498) Property changes on: ../trunk-jpl/src/wrappers/ProcessRifts ___________________________________________________________________ Added: svn:ignore ## -0,0 +1,2 ## +.deps +.dirstamp Index: ../trunk-jpl/src/wrappers/matlab/Makefile.am =================================================================== --- ../trunk-jpl/src/wrappers/matlab/Makefile.am (revision 22497) +++ ../trunk-jpl/src/wrappers/matlab/Makefile.am (revision 22498) @@ -57,8 +57,8 @@ MeshProfileIntersection_matlab.la\ PointCloudFindNeighbors_matlab.la\ PropagateFlagsFromConnectivity_matlab.la\ - TriMesh_matlab.la\ - TriMeshProcessRifts_matlab.la\ + Triangle_matlab.la\ + ProcessRifts_matlab.la\ Scotch_matlab.la if CHACO @@ -232,11 +232,11 @@ ShpRead_matlab_la_CXXFLAGS = ${AM_CXXFLAGS} 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) +Triangle_matlab_la_SOURCES = ../Triangle/Triangle.cpp +Triangle_matlab_la_CXXFLAGS = ${AM_CXXFLAGS} +Triangle_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) +ProcessRifts_matlab_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp +ProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS} +ProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB) #}}} Index: ../trunk-jpl/src/wrappers/python/Makefile.am =================================================================== --- ../trunk-jpl/src/wrappers/python/Makefile.am (revision 22497) +++ ../trunk-jpl/src/wrappers/python/Makefile.am (revision 22498) @@ -40,8 +40,8 @@ IssmConfig_python.la\ MeshProfileIntersection_python.la\ NodeConnectivity_python.la\ - TriMesh_python.la\ - TriMeshProcessRifts_python.la + Triangle_python.la\ + ProcessRifts_python.la endif #}}} #Flags and libraries {{{ @@ -140,11 +140,11 @@ NodeConnectivity_python_la_CXXFLAGS = ${AM_CXXFLAGS} 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) #}}} Index: ../trunk-jpl/src/m/mesh/triangle.js =================================================================== --- ../trunk-jpl/src/m/mesh/triangle.js (revision 22497) +++ ../trunk-jpl/src/m/mesh/triangle.js (revision 22498) @@ -1,7 +1,7 @@ function triangle(md){ //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 // unit). Riftname is an optional argument (Argus domain outline) describing rifts. @@ -35,7 +35,7 @@ var area=Math.pow(resolution,2); //Call mesher: - var return_array=TriMesh(md, domain, rifts, area); + var return_array=Triangle(md, domain, rifts, area); //Plug into md: md.mesh.elements=return_array[0]; Index: ../trunk-jpl/src/m/mesh/triangle2dvertical.m =================================================================== --- ../trunk-jpl/src/m/mesh/triangle2dvertical.m (revision 22497) +++ ../trunk-jpl/src/m/mesh/triangle2dvertical.m (revision 22498) @@ -1,7 +1,7 @@ function md=triangle(md,domainname,resolution) %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 % unit). Riftname is an optional argument (Argus domain outline) describing rifts. @@ -23,8 +23,8 @@ 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 orphan=find(~ismember([1:length(x)],sort(unique(elements(:))))); Index: ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py =================================================================== --- ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 22497) +++ ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 22498) @@ -1,5 +1,5 @@ import numpy as np -from TriMeshProcessRifts import TriMeshProcessRifts +from ProcessRifts import ProcessRifts from ContourToMesh import ContourToMesh from meshprocessoutsiderifts import meshprocessoutsiderifts from GetAreas import GetAreas @@ -21,7 +21,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) md.mesh.elements=md.mesh.elements.astype(int) md.mesh.x=md.mesh.x.reshape(-1) md.mesh.y=md.mesh.y.reshape(-1) @@ -28,7 +28,7 @@ md.mesh.segments=md.mesh.segments.astype(int) 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: numrifts=len(md.rifts.riftstruct) Index: ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m =================================================================== --- ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m (revision 22497) +++ ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m (revision 22498) @@ -24,9 +24,9 @@ end %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 %Fill in rest of fields: Index: ../trunk-jpl/src/m/mesh/argusmesh.m =================================================================== --- ../trunk-jpl/src/m/mesh/argusmesh.m (revision 22497) +++ ../trunk-jpl/src/m/mesh/argusmesh.m (revision 22498) @@ -8,7 +8,7 @@ % md=argusmesh(md,infile) % % Example: -% md=argusmesh(md,'TriMesh.exp') +% md=argusmesh(md,'Domain.exp') %some argument check: if nargin~=2 | nargout~=1, Index: ../trunk-jpl/src/m/mesh/triangle.py =================================================================== --- ../trunk-jpl/src/m/mesh/triangle.py (revision 22497) +++ ../trunk-jpl/src/m/mesh/triangle.py (revision 22498) @@ -1,7 +1,7 @@ import os.path import numpy as np from mesh2d import mesh2d -from TriMesh import TriMesh +from Triangle import Triangle from NodeConnectivity import NodeConnectivity from ElementConnectivity import ElementConnectivity import MatlabFuncs as m @@ -10,7 +10,7 @@ """ 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 unit). Riftname is an optional argument (Argus domain outline) describing rifts. @@ -47,9 +47,9 @@ if not os.path.exists(domainname): 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) md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int) Index: ../trunk-jpl/src/m/mesh/triangle.m =================================================================== --- ../trunk-jpl/src/m/mesh/triangle.m (revision 22497) +++ ../trunk-jpl/src/m/mesh/triangle.m (revision 22498) @@ -1,7 +1,7 @@ function md=triangle(md,domainname,varargin) %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 % unit). Riftname is an optional argument (Argus domain outline) describing rifts. @@ -41,8 +41,8 @@ error(['file "' domainname '" not found']); 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 removeorphans=1; Index: ../trunk-jpl/src/m/modules/TriMesh.m =================================================================== --- ../trunk-jpl/src/m/modules/TriMesh.m (revision 22497) +++ ../trunk-jpl/src/m/modules/TriMesh.m (nonexistent) @@ -1,21 +0,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: ../trunk-jpl/src/m/modules/TriMesh.py =================================================================== --- ../trunk-jpl/src/m/modules/TriMesh.py (revision 22497) +++ ../trunk-jpl/src/m/modules/TriMesh.py (nonexistent) @@ -1,21 +0,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: ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py =================================================================== --- ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py (revision 22497) +++ ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py (nonexistent) @@ -1,16 +0,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: ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m =================================================================== --- ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m (revision 22497) +++ ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m (nonexistent) @@ -1,17 +0,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: ../trunk-jpl/src/m/modules/ProcessRifts.m =================================================================== --- ../trunk-jpl/src/m/modules/ProcessRifts.m (nonexistent) +++ ../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: ../trunk-jpl/src/m/modules/Triangle.py =================================================================== --- ../trunk-jpl/src/m/modules/Triangle.py (nonexistent) +++ ../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: ../trunk-jpl/src/m/modules/Triangle.m =================================================================== --- ../trunk-jpl/src/m/modules/Triangle.m (nonexistent) +++ ../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: ../trunk-jpl/src/m/modules/ProcessRifts.py =================================================================== --- ../trunk-jpl/src/m/modules/ProcessRifts.py (nonexistent) +++ ../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