source: issm/oecreview/Archive/21724-22754/ISSM-22497-22498.diff@ 22755

Last change on this file since 22755 was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 150.8 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp (revision 22497)
4+++ ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp (nonexistent)
5@@ -1,201 +0,0 @@
6-/*!\file TriMeshx
7- * \brief: x code for TriMesh mesher
8- */
9-
10-/*Header files: {{{*/
11-#include "./TriMeshx.h"
12-#include "../../shared/shared.h"
13-#include "../../toolkits/toolkits.h"
14-/*ANSI_DECLARATORS needed to call triangle library: */
15-#if defined(_HAVE_TRIANGLE_)
16- #ifndef ANSI_DECLARATORS
17- #define ANSI_DECLARATORS
18- #endif
19- #include "triangle.h"
20- #undef ANSI_DECLARATORS
21-#endif
22-/*}}}*/
23-
24-void TriMeshx(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnsegs,Contours* domain,Contours* rifts,double area){
25-
26-#if !defined(_HAVE_TRIANGLE_)
27- _error_("triangle has not been installed");
28-#else
29- /*indexing: */
30- int i,j;
31-
32- /*output: */
33- int *index = NULL;
34- double *x = NULL;
35- double *y = NULL;
36- int *segments = NULL;
37- int *segmentmarkerlist = NULL;
38-
39- /*intermediary: */
40- int counter,counter2,backcounter;
41- Contour<IssmPDouble> *contour = NULL;
42-
43- /* Triangle structures needed to call Triangle library routines: */
44- struct triangulateio in,out;
45- char options[256];
46-
47- /*Create initial triangulation to call triangulate(). First number of points:*/
48- in.numberofpoints=0;
49- for (i=0;i<domain->Size();i++){
50- contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
51- in.numberofpoints+=contour->nods-1;
52- }
53- for (i=0;i<rifts->Size();i++){
54- contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
55- in.numberofpoints+=contour->nods;
56- }
57-
58- /*number of point attributes: */
59- in.numberofpointattributes=1;
60-
61- /*fill in the point list: */
62- in.pointlist = xNew<REAL>(in.numberofpoints*2);
63-
64- counter=0;
65- for (i=0;i<domain->Size();i++){
66- contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
67- for (j=0;j<contour->nods-1;j++){
68- in.pointlist[2*counter+0]=contour->x[j];
69- in.pointlist[2*counter+1]=contour->y[j];
70- counter++;
71- }
72- }
73- for (i=0;i<rifts->Size();i++){
74- contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
75- for (j=0;j<contour->nods;j++){
76- in.pointlist[2*counter+0]=contour->x[j];
77- in.pointlist[2*counter+1]=contour->y[j];
78- counter++;
79- }
80- }
81-
82- /*fill in the point attribute list: */
83- in.pointattributelist = xNew<REAL>(in.numberofpoints*in.numberofpointattributes);
84- for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0;
85-
86- /*fill in the point marker list: */
87- in.pointmarkerlist = xNew<int>(in.numberofpoints);
88- for(i=0;i<in.numberofpoints;i++) in.pointmarkerlist[i] = 0;
89-
90- /*Build segments. First figure out number of segments: holes and closed outlines have as many segments as vertices: */
91- in.numberofsegments=0;
92- for (i=0;i<domain->Size();i++){
93- contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
94- in.numberofsegments+=contour->nods-1;
95- }
96- for(i=0;i<rifts->Size();i++){
97- contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
98- /*for rifts, we have one less segment as we have vertices*/
99- in.numberofsegments+=contour->nods-1;
100- }
101-
102- in.segmentlist = xNew<int>(in.numberofsegments*2);
103- in.segmentmarkerlist = xNewZeroInit<int>(in.numberofsegments);
104- counter=0;
105- backcounter=0;
106- for (i=0;i<domain->Size();i++){
107- contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
108- for (j=0;j<contour->nods-2;j++){
109- in.segmentlist[2*counter+0]=counter;
110- in.segmentlist[2*counter+1]=counter+1;
111- in.segmentmarkerlist[counter]=0;
112- counter++;
113- }
114- /*Close this profile: */
115- in.segmentlist[2*counter+0]=counter;
116- in.segmentlist[2*counter+1]=backcounter;
117- in.segmentmarkerlist[counter]=0;
118- counter++;
119- backcounter=counter;
120- }
121- counter2=counter;
122- for (i=0;i<rifts->Size();i++){
123- contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
124- for (j=0;j<(contour->nods-1);j++){
125- in.segmentlist[2*counter2+0]=counter;
126- in.segmentlist[2*counter2+1]=counter+1;
127- in.segmentmarkerlist[counter2]=2+i;
128- counter2++;
129- counter++;
130- }
131- counter++;
132- }
133-
134- /*Build regions: */
135- in.numberofregions = 0;
136-
137- /*Build holes: */
138- in.numberofholes = domain->Size()-1; /*everything is a hole, but for the first profile.*/
139- if(in.numberofholes){
140- in.holelist = xNew<REAL>(in.numberofholes*2);
141- for (i=0;i<domain->Size()-1;i++){
142- contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
143- GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
144- }
145- }
146-
147- /* Make necessary initializations so that Triangle can return a triangulation in `out': */
148- out.pointlist = (REAL*)NULL;
149- out.pointattributelist = (REAL*)NULL;
150- out.pointmarkerlist = (int *)NULL;
151- out.trianglelist = (int *)NULL;
152- out.triangleattributelist = (REAL*)NULL;
153- out.neighborlist = (int *)NULL;
154- out.segmentlist = (int *)NULL;
155- out.segmentmarkerlist = (int *)NULL;
156- out.edgelist = (int *)NULL;
157- out.edgemarkerlist = (int *)NULL;
158-
159- /* Triangulate the points:. Switches are chosen to read and write a */
160- /* PSLG (p), preserve the convex hull (c), number everything from */
161- /* zero (z), assign a regional attribute to each element (A), and */
162- /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
163- /* neighbor list (n). */
164- sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
165- triangulate(options, &in, &out, NULL);
166- /*report(&out, 0, 1, 1, 1, 1, 0);*/
167-
168- /*Allocate index, x and y: */
169- index=xNew<int>(3*out.numberoftriangles);
170- x=xNew<double>(out.numberofpoints);
171- y=xNew<double>(out.numberofpoints);
172- segments=xNew<int>(3*out.numberofsegments);
173- segmentmarkerlist=xNew<int>(out.numberofsegments);
174-
175- for (i = 0; i< out.numberoftriangles; i++) {
176- for (j = 0; j < out.numberofcorners; j++) {
177- index[3*i+j]=(int)out.trianglelist[i*out.numberofcorners+j]+1;
178- }
179- }
180- for (i = 0; i< out.numberofpoints; i++){
181- x[i]=(double)out.pointlist[i*2+0];
182- y[i]=(double)out.pointlist[i*2+1];
183- }
184- for (i = 0; i<out.numberofsegments;i++){
185- segments[3*i+0]=(int)out.segmentlist[i*2+0]+1;
186- segments[3*i+1]=(int)out.segmentlist[i*2+1]+1;
187- segmentmarkerlist[i]=(int)out.segmentmarkerlist[i];
188- }
189-
190- /*Associate elements with segments: */
191- AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
192-
193- /*Order segments so that their normals point outside the domain: */
194- OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
195-
196- /*Output : */
197- *pindex=index;
198- *px=x;
199- *py=y;
200- *psegments=segments;
201- *psegmentmarkerlist=segmentmarkerlist;
202- *pnels=out.numberoftriangles;
203- *pnods=out.numberofpoints;
204- *pnsegs=out.numberofsegments;
205-#endif
206-}
207Index: ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h
208===================================================================
209--- ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h (revision 22497)
210+++ ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h (nonexistent)
211@@ -1,13 +0,0 @@
212-/*!\file: TriMeshx.h
213- * \brief header file for TriMeshx module
214- */
215-
216-#ifndef _TRIMESHX_H_
217-#define _TRIMESHX_H_
218-
219-#include <string.h>
220-#include "../../classes/classes.h"
221-
222-/* local prototypes: */
223-void TriMeshx(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnseg,Contours* domain,Contours* rifts,double area);
224-#endif /* _TRIMESHX_H */
225Index: ../trunk-jpl/src/c/modules/TriMeshx
226===================================================================
227--- ../trunk-jpl/src/c/modules/TriMeshx (revision 22497)
228+++ ../trunk-jpl/src/c/modules/TriMeshx (nonexistent)
229
230Property changes on: ../trunk-jpl/src/c/modules/TriMeshx
231___________________________________________________________________
232Deleted: svn:ignore
233## -1,2 +0,0 ##
234-.deps
235-.dirstamp
236Index: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h
237===================================================================
238--- ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h (revision 22497)
239+++ ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h (nonexistent)
240@@ -1,12 +0,0 @@
241-/*!\file: TriMeshProcessRiftsx.h
242- * \brief header file for TriMeshProcessRifts module
243- */
244-
245-#ifndef _TRIMESHPROCESSRIFTX_H
246-#define _TRIMESHPROCESSRIFTX_H
247-
248-class RiftStruct;
249-
250-void TriMeshProcessRiftsx(int** pindex,int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct);
251-
252-#endif /* _TRIMESHPROCESSRIFTX_H*/
253Index: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp
254===================================================================
255--- ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp (revision 22497)
256+++ ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp (nonexistent)
257@@ -1,78 +0,0 @@
258-/*!\file: TriMeshProcessRifts.cpp
259- * \brief split a mesh where a rift (or fault) is present
260- */
261-
262-#include "./TriMeshProcessRiftsx.h"
263-#include "../../classes/RiftStruct.h"
264-#include "../../shared/shared.h"
265-#include "../../toolkits/toolkits.h"
266-
267-void TriMeshProcessRiftsx(int** pindex, int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct){
268-
269- /*Output*/
270- int numrifts,numrifts0;
271- int *riftsnumsegments = NULL;
272- int **riftssegments = NULL;
273- int *riftsnumpairs = NULL;
274- int **riftspairs = NULL;
275- int *riftstips = NULL;
276- double **riftspenaltypairs = NULL;
277- int *riftsnumpenaltypairs = NULL;
278-
279- /*Recover initial mesh*/
280- int nel = *pnel;
281- int *index = *pindex;
282- double *x = *px;
283- double *y = *py;
284- int nods = *pnods;
285- int *segments = *psegments;
286- int *segmentmarkers = *psegmentmarkers;
287- int num_seg = *pnum_seg;
288-
289- /*Intermediary*/
290- int riftflag;
291-
292- /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
293- *all the nodes of this element belong to the segments (tends to happen when there are corners: */
294- RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
295-
296- /*Figure out if we have rifts, and how many: */
297- IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
298-
299- if(!riftflag) _error_("No rift present in mesh");
300-
301- /*Split mesh*/
302- SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
303-
304- /*Order segments so that their normals point outside the domain: */
305- OrderSegments(&segments,num_seg, index,nel);
306-
307- /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
308- *segmentmarkerlist:*/
309- SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
310-
311- /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
312- PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
313-
314- /*Order rifts so that they start from one tip, go to the other tip, and back: */
315- OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
316-
317- /*Create penalty pairs, used by Imp: */
318- PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
319-
320- /*Create Riftstruct*/
321- RiftStruct* riftstruct = new RiftStruct(numrifts,riftsnumsegments,riftssegments,riftsnumpairs,riftspairs,riftsnumpenaltypairs,riftspenaltypairs,riftstips);
322-
323- /*Assign output pointers for mesh*/
324- *pnel = nel;
325- *pindex = index;
326- *px = x;
327- *py = y;
328- *pnods = nods;
329- *psegments = segments;
330- *psegmentmarkers = segmentmarkers;
331- *pnum_seg = num_seg;
332-
333- /*Assign output pointers for rifts*/
334- *priftstruct = riftstruct;
335-}
336Index: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx
337===================================================================
338--- ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx (revision 22497)
339+++ ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx (nonexistent)
340
341Property changes on: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx
342___________________________________________________________________
343Deleted: svn:ignore
344## -1,2 +0,0 ##
345-.deps
346-.dirstamp
347Index: ../trunk-jpl/src/c/modules/modules.h
348===================================================================
349--- ../trunk-jpl/src/c/modules/modules.h (revision 22497)
350+++ ../trunk-jpl/src/c/modules/modules.h (revision 22498)
351@@ -90,8 +90,8 @@
352 #include "./SystemMatricesx/SystemMatricesx.h"
353 #include "./SpcNodesx/SpcNodesx.h"
354 #include "./SurfaceAreax/SurfaceAreax.h"
355-#include "./TriMeshx/TriMeshx.h"
356-#include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h"
357+#include "./Trianglex/Trianglex.h"
358+#include "./ProcessRiftsx/ProcessRiftsx.h"
359 #include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
360 #include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h"
361 #include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"
362Index: ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h
363===================================================================
364--- ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h (nonexistent)
365+++ ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h (revision 22498)
366@@ -0,0 +1,12 @@
367+/*!\file: ProcessRiftsx.h
368+ * \brief header file for ProcessRifts module
369+ */
370+
371+#ifndef _PROCESSRIFTX_H
372+#define _PROCESSRIFTX_H
373+
374+class RiftStruct;
375+
376+void ProcessRiftsx(int** pindex,int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct);
377+
378+#endif /* _PROCESSRIFTX_H*/
379Index: ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp
380===================================================================
381--- ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp (nonexistent)
382+++ ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp (revision 22498)
383@@ -0,0 +1,78 @@
384+/*!\file: ProcessRifts.cpp
385+ * \brief split a mesh where a rift (or fault) is present
386+ */
387+
388+#include "./ProcessRiftsx.h"
389+#include "../../classes/RiftStruct.h"
390+#include "../../shared/shared.h"
391+#include "../../toolkits/toolkits.h"
392+
393+void ProcessRiftsx(int** pindex, int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct){
394+
395+ /*Output*/
396+ int numrifts,numrifts0;
397+ int *riftsnumsegments = NULL;
398+ int **riftssegments = NULL;
399+ int *riftsnumpairs = NULL;
400+ int **riftspairs = NULL;
401+ int *riftstips = NULL;
402+ double **riftspenaltypairs = NULL;
403+ int *riftsnumpenaltypairs = NULL;
404+
405+ /*Recover initial mesh*/
406+ int nel = *pnel;
407+ int *index = *pindex;
408+ double *x = *px;
409+ double *y = *py;
410+ int nods = *pnods;
411+ int *segments = *psegments;
412+ int *segmentmarkers = *psegmentmarkers;
413+ int num_seg = *pnum_seg;
414+
415+ /*Intermediary*/
416+ int riftflag;
417+
418+ /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
419+ *all the nodes of this element belong to the segments (tends to happen when there are corners: */
420+ RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
421+
422+ /*Figure out if we have rifts, and how many: */
423+ IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
424+
425+ if(!riftflag) _error_("No rift present in mesh");
426+
427+ /*Split mesh*/
428+ SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
429+
430+ /*Order segments so that their normals point outside the domain: */
431+ OrderSegments(&segments,num_seg, index,nel);
432+
433+ /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
434+ *segmentmarkerlist:*/
435+ SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
436+
437+ /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
438+ PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
439+
440+ /*Order rifts so that they start from one tip, go to the other tip, and back: */
441+ OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
442+
443+ /*Create penalty pairs, used by Imp: */
444+ PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
445+
446+ /*Create Riftstruct*/
447+ RiftStruct* riftstruct = new RiftStruct(numrifts,riftsnumsegments,riftssegments,riftsnumpairs,riftspairs,riftsnumpenaltypairs,riftspenaltypairs,riftstips);
448+
449+ /*Assign output pointers for mesh*/
450+ *pnel = nel;
451+ *pindex = index;
452+ *px = x;
453+ *py = y;
454+ *pnods = nods;
455+ *psegments = segments;
456+ *psegmentmarkers = segmentmarkers;
457+ *pnum_seg = num_seg;
458+
459+ /*Assign output pointers for rifts*/
460+ *priftstruct = riftstruct;
461+}
462Index: ../trunk-jpl/src/c/modules/ProcessRiftsx
463===================================================================
464--- ../trunk-jpl/src/c/modules/ProcessRiftsx (nonexistent)
465+++ ../trunk-jpl/src/c/modules/ProcessRiftsx (revision 22498)
466
467Property changes on: ../trunk-jpl/src/c/modules/ProcessRiftsx
468___________________________________________________________________
469Added: svn:ignore
470## -0,0 +1,2 ##
471+.deps
472+.dirstamp
473Index: ../trunk-jpl/src/c/modules/Trianglex/Trianglex.h
474===================================================================
475--- ../trunk-jpl/src/c/modules/Trianglex/Trianglex.h (nonexistent)
476+++ ../trunk-jpl/src/c/modules/Trianglex/Trianglex.h (revision 22498)
477@@ -0,0 +1,13 @@
478+/*!\file: Trianglex.h
479+ * \brief header file for Trianglex module
480+ */
481+
482+#ifndef _TRIANGLEX_H_
483+#define _TRIANGLEX_H_
484+
485+#include <string.h>
486+#include "../../classes/classes.h"
487+
488+/* local prototypes: */
489+void Trianglex(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnseg,Contours* domain,Contours* rifts,double area);
490+#endif /* _TRIANGLEX_H */
491Index: ../trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp
492===================================================================
493--- ../trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp (nonexistent)
494+++ ../trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp (revision 22498)
495@@ -0,0 +1,201 @@
496+/*!\file Trianglex
497+ * \brief: x code for Triangle mesher
498+ */
499+
500+/*Header files: {{{*/
501+#include "./Trianglex.h"
502+#include "../../shared/shared.h"
503+#include "../../toolkits/toolkits.h"
504+/*ANSI_DECLARATORS needed to call triangle library: */
505+#if defined(_HAVE_TRIANGLE_)
506+ #ifndef ANSI_DECLARATORS
507+ #define ANSI_DECLARATORS
508+ #endif
509+ #include "triangle.h"
510+ #undef ANSI_DECLARATORS
511+#endif
512+/*}}}*/
513+
514+void Trianglex(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnsegs,Contours* domain,Contours* rifts,double area){
515+
516+#if !defined(_HAVE_TRIANGLE_)
517+ _error_("triangle has not been installed");
518+#else
519+ /*indexing: */
520+ int i,j;
521+
522+ /*output: */
523+ int *index = NULL;
524+ double *x = NULL;
525+ double *y = NULL;
526+ int *segments = NULL;
527+ int *segmentmarkerlist = NULL;
528+
529+ /*intermediary: */
530+ int counter,counter2,backcounter;
531+ Contour<IssmPDouble> *contour = NULL;
532+
533+ /* Triangle structures needed to call Triangle library routines: */
534+ struct triangulateio in,out;
535+ char options[256];
536+
537+ /*Create initial triangulation to call triangulate(). First number of points:*/
538+ in.numberofpoints=0;
539+ for (i=0;i<domain->Size();i++){
540+ contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
541+ in.numberofpoints+=contour->nods-1;
542+ }
543+ for (i=0;i<rifts->Size();i++){
544+ contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
545+ in.numberofpoints+=contour->nods;
546+ }
547+
548+ /*number of point attributes: */
549+ in.numberofpointattributes=1;
550+
551+ /*fill in the point list: */
552+ in.pointlist = xNew<REAL>(in.numberofpoints*2);
553+
554+ counter=0;
555+ for (i=0;i<domain->Size();i++){
556+ contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
557+ for (j=0;j<contour->nods-1;j++){
558+ in.pointlist[2*counter+0]=contour->x[j];
559+ in.pointlist[2*counter+1]=contour->y[j];
560+ counter++;
561+ }
562+ }
563+ for (i=0;i<rifts->Size();i++){
564+ contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
565+ for (j=0;j<contour->nods;j++){
566+ in.pointlist[2*counter+0]=contour->x[j];
567+ in.pointlist[2*counter+1]=contour->y[j];
568+ counter++;
569+ }
570+ }
571+
572+ /*fill in the point attribute list: */
573+ in.pointattributelist = xNew<REAL>(in.numberofpoints*in.numberofpointattributes);
574+ for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0;
575+
576+ /*fill in the point marker list: */
577+ in.pointmarkerlist = xNew<int>(in.numberofpoints);
578+ for(i=0;i<in.numberofpoints;i++) in.pointmarkerlist[i] = 0;
579+
580+ /*Build segments. First figure out number of segments: holes and closed outlines have as many segments as vertices: */
581+ in.numberofsegments=0;
582+ for (i=0;i<domain->Size();i++){
583+ contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
584+ in.numberofsegments+=contour->nods-1;
585+ }
586+ for(i=0;i<rifts->Size();i++){
587+ contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
588+ /*for rifts, we have one less segment as we have vertices*/
589+ in.numberofsegments+=contour->nods-1;
590+ }
591+
592+ in.segmentlist = xNew<int>(in.numberofsegments*2);
593+ in.segmentmarkerlist = xNewZeroInit<int>(in.numberofsegments);
594+ counter=0;
595+ backcounter=0;
596+ for (i=0;i<domain->Size();i++){
597+ contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
598+ for (j=0;j<contour->nods-2;j++){
599+ in.segmentlist[2*counter+0]=counter;
600+ in.segmentlist[2*counter+1]=counter+1;
601+ in.segmentmarkerlist[counter]=0;
602+ counter++;
603+ }
604+ /*Close this profile: */
605+ in.segmentlist[2*counter+0]=counter;
606+ in.segmentlist[2*counter+1]=backcounter;
607+ in.segmentmarkerlist[counter]=0;
608+ counter++;
609+ backcounter=counter;
610+ }
611+ counter2=counter;
612+ for (i=0;i<rifts->Size();i++){
613+ contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
614+ for (j=0;j<(contour->nods-1);j++){
615+ in.segmentlist[2*counter2+0]=counter;
616+ in.segmentlist[2*counter2+1]=counter+1;
617+ in.segmentmarkerlist[counter2]=2+i;
618+ counter2++;
619+ counter++;
620+ }
621+ counter++;
622+ }
623+
624+ /*Build regions: */
625+ in.numberofregions = 0;
626+
627+ /*Build holes: */
628+ in.numberofholes = domain->Size()-1; /*everything is a hole, but for the first profile.*/
629+ if(in.numberofholes){
630+ in.holelist = xNew<REAL>(in.numberofholes*2);
631+ for (i=0;i<domain->Size()-1;i++){
632+ contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
633+ GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
634+ }
635+ }
636+
637+ /* Make necessary initializations so that Triangle can return a triangulation in `out': */
638+ out.pointlist = (REAL*)NULL;
639+ out.pointattributelist = (REAL*)NULL;
640+ out.pointmarkerlist = (int *)NULL;
641+ out.trianglelist = (int *)NULL;
642+ out.triangleattributelist = (REAL*)NULL;
643+ out.neighborlist = (int *)NULL;
644+ out.segmentlist = (int *)NULL;
645+ out.segmentmarkerlist = (int *)NULL;
646+ out.edgelist = (int *)NULL;
647+ out.edgemarkerlist = (int *)NULL;
648+
649+ /* Triangulate the points:. Switches are chosen to read and write a */
650+ /* PSLG (p), preserve the convex hull (c), number everything from */
651+ /* zero (z), assign a regional attribute to each element (A), and */
652+ /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
653+ /* neighbor list (n). */
654+ sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
655+ triangulate(options, &in, &out, NULL);
656+ /*report(&out, 0, 1, 1, 1, 1, 0);*/
657+
658+ /*Allocate index, x and y: */
659+ index=xNew<int>(3*out.numberoftriangles);
660+ x=xNew<double>(out.numberofpoints);
661+ y=xNew<double>(out.numberofpoints);
662+ segments=xNew<int>(3*out.numberofsegments);
663+ segmentmarkerlist=xNew<int>(out.numberofsegments);
664+
665+ for (i = 0; i< out.numberoftriangles; i++) {
666+ for (j = 0; j < out.numberofcorners; j++) {
667+ index[3*i+j]=(int)out.trianglelist[i*out.numberofcorners+j]+1;
668+ }
669+ }
670+ for (i = 0; i< out.numberofpoints; i++){
671+ x[i]=(double)out.pointlist[i*2+0];
672+ y[i]=(double)out.pointlist[i*2+1];
673+ }
674+ for (i = 0; i<out.numberofsegments;i++){
675+ segments[3*i+0]=(int)out.segmentlist[i*2+0]+1;
676+ segments[3*i+1]=(int)out.segmentlist[i*2+1]+1;
677+ segmentmarkerlist[i]=(int)out.segmentmarkerlist[i];
678+ }
679+
680+ /*Associate elements with segments: */
681+ AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
682+
683+ /*Order segments so that their normals point outside the domain: */
684+ OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
685+
686+ /*Output : */
687+ *pindex=index;
688+ *px=x;
689+ *py=y;
690+ *psegments=segments;
691+ *psegmentmarkerlist=segmentmarkerlist;
692+ *pnels=out.numberoftriangles;
693+ *pnods=out.numberofpoints;
694+ *pnsegs=out.numberofsegments;
695+#endif
696+}
697Index: ../trunk-jpl/src/c/modules/Trianglex
698===================================================================
699--- ../trunk-jpl/src/c/modules/Trianglex (nonexistent)
700+++ ../trunk-jpl/src/c/modules/Trianglex (revision 22498)
701
702Property changes on: ../trunk-jpl/src/c/modules/Trianglex
703___________________________________________________________________
704Added: svn:ignore
705## -0,0 +1,2 ##
706+.deps
707+.dirstamp
708Index: ../trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp
709===================================================================
710--- ../trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp (revision 22497)
711+++ ../trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp (revision 22498)
712@@ -1,5 +1,5 @@
713-/*!\file TriMeshx
714- * \brief: x code for TriMesh mesher
715+/*!\file CoordinateSystemTransformx
716+ * \brief: x code for CoordinateSystemTransformx
717 */
718
719 /*Header files*/
720Index: ../trunk-jpl/src/c/Makefile.am
721===================================================================
722--- ../trunk-jpl/src/c/Makefile.am (revision 22497)
723+++ ../trunk-jpl/src/c/Makefile.am (revision 22498)
724@@ -560,13 +560,13 @@
725 modules_sources= ./shared/Threads/LaunchThread.cpp\
726 ./shared/Threads/PartitionRange.cpp\
727 ./shared/Exp/exp.cpp\
728- ./shared/TriMesh/AssociateSegmentToElement.cpp\
729- ./shared/TriMesh/GridInsideHole.cpp\
730- ./shared/TriMesh/OrderSegments.cpp\
731- ./shared/TriMesh/SplitMeshForRifts.cpp\
732- ./shared/TriMesh/TriMeshUtils.cpp\
733- ./modules/TriMeshx/TriMeshx.cpp\
734- ./modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp\
735+ ./shared/Triangle/AssociateSegmentToElement.cpp\
736+ ./shared/Triangle/GridInsideHole.cpp\
737+ ./shared/Triangle/OrderSegments.cpp\
738+ ./shared/Triangle/SplitMeshForRifts.cpp\
739+ ./shared/Triangle/TriangleUtils.cpp\
740+ ./modules/Trianglex/Trianglex.cpp\
741+ ./modules/ProcessRiftsx/ProcessRiftsx.cpp\
742 ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp\
743 ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp\
744 ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
745Index: ../trunk-jpl/src/c/main/issm.js
746===================================================================
747--- ../trunk-jpl/src/c/main/issm.js (revision 22497)
748+++ ../trunk-jpl/src/c/main/issm.js (revision 22498)
749@@ -10,7 +10,7 @@
750 var dbinaryPtr= Module._malloc(nb); var binHeap = new Uint8Array(Module.HEAPU8.buffer,dbinaryPtr,nb);
751 binHeap.set(new Uint8Array(dbinary.buffer)); var binary=binHeap.byteOffset;
752
753- //Declare TriMesh module:
754+ //Declare module:
755 issm= Module.cwrap('main','number',['number','number']);
756
757 //Call issm:
758Index: ../trunk-jpl/src/c/shared/TriMesh/SplitMeshForRifts.cpp
759===================================================================
760--- ../trunk-jpl/src/c/shared/TriMesh/SplitMeshForRifts.cpp (revision 22497)
761+++ ../trunk-jpl/src/c/shared/TriMesh/SplitMeshForRifts.cpp (nonexistent)
762@@ -1,96 +0,0 @@
763-/*
764- * SplitMeshForRifts.c:
765- */
766-#include "./trimesh.h"
767-#include "../MemOps/MemOps.h"
768-
769-int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist){
770-
771- /*Some notes on dimensions:
772- index of size nelx3
773- x and y of size nodsx1
774- segments of size nsegsx3*/
775-
776- int i,j,k,l;
777- int node;
778- int el;
779- int nriftsegs;
780- int* riftsegments=NULL;
781- int* flags=NULL;
782- int NumGridElementListOnOneSideOfRift;
783- int* GridElementListOnOneSideOfRift=NULL;
784-
785- /*Recover input: */
786- int nel = *pnel;
787- int *index = *pindex;
788- int nods = *pnods;
789- double *x = *px;
790- double *y = *py;
791- int nsegs = *pnsegs;
792- int *segments = *psegments;
793- int *segmentmarkerlist = *psegmentmarkerlist;
794-
795- /*Establish list of segments that belong to a rift: */
796- /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
797- RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);
798-
799- /*Go through all nodes of the rift segments, and start splitting the mesh: */
800- flags=xNewZeroInit<int>(nods); //to make sure we don't split the same nodes twice!
801- for (i=0;i<nriftsegs;i++){
802- for (j=0;j<2;j++){
803-
804- node=riftsegments[4*i+j+2];
805- if(flags[node-1]){
806- /*This node was already split, skip:*/
807- continue;
808- }
809- else{
810- flags[node-1]=1;
811- }
812-
813- if(IsGridOnRift(riftsegments,nriftsegs,node)){
814-
815- DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
816-
817- /*Summary: we have for node, a list of elements
818- * (GridElementListOnOneSideOfRift, of size
819- * NumGridElementListOnOneSideOfRift) that all contain node
820- *and that are on the same side of the rift. For all these
821- elements, we clone node into another node, and we swap all
822- instances of node in the triangulation *for those elements, to the
823- new node.*/
824-
825- //create new node
826- x=xReNew<double>(x,nods,nods+1);
827- y=xReNew<double>(y,nods,nods+1);
828- x[nods]=x[node-1]; //matlab indexing
829- y[nods]=y[node-1]; //matlab indexing
830-
831- //augment number of nodes
832- nods++;
833-
834- //change elements owning this node
835- for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
836- el=GridElementListOnOneSideOfRift[k];
837- for (l=0;l<3;l++){
838- if (index[3*el+l]==node) index[3*el+l]=nods; //again, matlab indexing.
839- }
840- }
841- }// if(IsGridOnRift(riftsegments,nriftsegs,node))
842- } //for(j=0;j<2;j++)
843- } //for (i=0;i<nriftsegs;i++)
844-
845- /*update segments: they got modified completely by adding new nodes.*/
846- UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
847-
848- /*Assign output pointers: */
849- *pnel=nel;
850- *pindex=index;
851- *pnods=nods;
852- *px=x;
853- *py=y;
854- *pnsegs=nsegs;
855- *psegments=segments;
856- *psegmentmarkerlist=segmentmarkerlist;
857- return 1;
858-}
859Index: ../trunk-jpl/src/c/shared/TriMesh/GridInsideHole.cpp
860===================================================================
861--- ../trunk-jpl/src/c/shared/TriMesh/GridInsideHole.cpp (revision 22497)
862+++ ../trunk-jpl/src/c/shared/TriMesh/GridInsideHole.cpp (nonexistent)
863@@ -1,51 +0,0 @@
864-/*
865- * GridInsideHole.c:
866- * from a convex set of points, figure out a point that for sure lies inside the profile.
867- */
868-
869-#include <math.h>
870-
871-#include "./trimesh.h"
872-#include "../Exp/exp.h"
873-
874-#undef M_PI
875-#define M_PI 3.141592653589793238462643
876-
877-int GridInsideHole(double* px0,double* py0,int n,double* x,double* y){
878-
879- double flag=0.0;
880- double xA,xB,xC,xD,xE;
881- double yA,yB,yC,yD,yE;
882-
883- /*Take first and last vertices: */
884- xA=x[0];
885- yA=y[0];
886- xB=x[n-1];
887- yB=y[n-1];
888-
889- /*Figure out middle of segment [A B]: */
890- xC=(xA+xB)/2;
891- yC=(yA+yB)/2;
892-
893- /*D and E are on each side of segment [A B], on the median line between segment [A B],
894- *at an angle of 10 degree (less than the minimum 30 enforced by the quality of the mesh: */
895- xD=xC+tan(10./180.*M_PI)*(yC-yA);
896- yD=yC+tan(10./180.*M_PI)*(xA-xC);
897- xE=xC-tan(10./180.*M_PI)*(yC-yA);
898- yE=yC-tan(10./180.*M_PI)*(xA-xC);
899-
900- /*Either E or D is inside profile (x,y): */
901- IsInPolySerial(&flag,&xD,&yD,1,x,y,n,2);
902- /*FIXME: used to be 'flag' and not '!flag', check*/
903- if(!flag){
904- /*D is inside the poly: */
905- *px0=xD;
906- *py0=yD;
907- }
908- else{
909- /*E is inside the poly: */
910- *px0=xE;
911- *py0=yE;
912- }
913- return 1;
914-}
915Index: ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp
916===================================================================
917--- ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp (revision 22497)
918+++ ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp (nonexistent)
919@@ -1,24 +0,0 @@
920-/*!\file: AssociateSegmentToElement.cpp
921- * \brief for each segment, look for the corresponding element.
922- */
923-
924-#include "./trimesh.h"
925-
926-int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel){
927-
928- /*node indices: */
929- int A,B;
930-
931- /*Recover segments: */
932- int* segments=*psegments;
933-
934- for(int i=0;i<nseg;i++){
935- A=segments[3*i+0];
936- B=segments[3*i+1];
937- segments[3*i+2]=FindElement(A,B,index,nel)+1; //matlab indexing.
938- }
939-
940- /*Assign output pointers: */
941- *psegments=segments;
942- return 1;
943-}
944Index: ../trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp
945===================================================================
946--- ../trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp (revision 22497)
947+++ ../trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp (nonexistent)
948@@ -1,912 +0,0 @@
949-/*
950- * TriMeshUtils: mesh manipulation routines:
951- */
952-
953-#include <stdio.h>
954-
955-#include "./trimesh.h"
956-#include "../Exceptions/exceptions.h"
957-#include "../MemOps/MemOps.h"
958-
959-#define RIFTPENALTYPAIRSWIDTH 8
960-int IsGridOnRift(int* riftsegments, int nriftsegs, int node){/*{{{*/
961-
962- /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
963- *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).*/
964-
965- int i;
966- int j;
967- int count;
968-
969- count=0;
970- for (i=0;i<nriftsegs;i++){
971- for (j=0;j<2;j++){
972- if ((*(riftsegments+4*i+2+j))==node) count++;
973- }
974- }
975- if (count==2){
976- return 1;
977- }
978- else{
979- return 0;
980- }
981-}/*}}}*/
982-int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){/*{{{*/
983-
984- /*From a node, recover all the elements that are connected to it: */
985- int i,j;
986- int noerr=1;
987-
988- int max_number_elements=12;
989- int current_size;
990- int NumGridElements;
991- int* GridElements=NULL;
992- int* GridElementsRealloc=NULL;
993-
994- /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own
995- * the node. We start by allocating GridElements with that size, and realloc
996- * more if needed.*/
997-
998- current_size=max_number_elements;
999- NumGridElements=0;
1000- GridElements=xNew<int>(max_number_elements);
1001-
1002- for (i=0;i<nel;i++){
1003- for (j=0;j<3;j++){
1004- if (index[3*i+j]==node){
1005- if (NumGridElements<=(current_size-1)){
1006- GridElements[NumGridElements]=i;
1007- NumGridElements++;
1008- break;
1009- }
1010- else{
1011- /*Reallocate another max_number_elements slots in the GridElements: */
1012- GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
1013- if (!GridElementsRealloc){
1014- noerr=0;
1015- goto cleanup_and_return;
1016- }
1017- current_size+=max_number_elements;
1018- GridElements=GridElementsRealloc;
1019- GridElements[NumGridElements]=i;
1020- NumGridElements++;
1021- break;
1022- }
1023- }
1024- }
1025- }
1026- cleanup_and_return:
1027- if(!noerr){
1028- xDelete<int>(GridElements);
1029- }
1030- /*Allocate return pointers: */
1031- *pGridElements=GridElements;
1032- *pNumGridElements=NumGridElements;
1033- return noerr;
1034-}/*}}}*/
1035-int IsNeighbor(int el1,int el2,int* index){/*{{{*/
1036- /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
1037- int i,j;
1038- int count=0;
1039- for (i=0;i<3;i++){
1040- for (j=0;j<3;j++){
1041- if (index[3*el1+i]==index[3*el2+j])count++;
1042- }
1043- }
1044- if (count==2){
1045- return 1;
1046- }
1047- else{
1048- return 0;
1049- }
1050-}/*}}}*/
1051-int IsOnRift(int el,int nriftsegs,int* riftsegments){/*{{{*/
1052- /*From a list of elements segments, figure out if el belongs to it: */
1053- int i;
1054- for (i=0;i<nriftsegs;i++){
1055- if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
1056- return 1;
1057- }
1058- }
1059- return 0;
1060-}/*}}}*/
1061-void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){/*{{{*/
1062-
1063- int i,counter;
1064- int el,el2;
1065-
1066- int nriftsegs;
1067- int* riftsegments=NULL;
1068- int* riftsegments_uncompressed=NULL;
1069- int element_nodes[3];
1070-
1071- /*Allocate segmentflags: */
1072- riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
1073-
1074- /*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary
1075- *or a hole: */
1076- nriftsegs=0;
1077- for (i=0;i<nsegs;i++){
1078- el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
1079- /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
1080- *then proceed to find another element that owns the segment. If we don't find it, we know
1081- *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
1082- element_nodes[0]=*(index+3*el+0);
1083- element_nodes[1]=*(index+3*el+1);
1084- element_nodes[2]=*(index+3*el+2);
1085-
1086- index[3*el+0]=-1;
1087- index[3*el+1]=-1;
1088- index[3*el+2]=-1;
1089-
1090- el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
1091-
1092- /*Restore index: */
1093- index[3*el+0]=element_nodes[0];
1094- index[3*el+1]=element_nodes[1];
1095- index[3*el+2]=element_nodes[2];
1096-
1097- if (el2!=-1){
1098- /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
1099- riftsegments_uncompressed[5*i+0]=1;
1100- riftsegments_uncompressed[5*i+1]=el;
1101- riftsegments_uncompressed[5*i+2]=el2;
1102- riftsegments_uncompressed[5*i+3]=segments[3*i+0];
1103- riftsegments_uncompressed[5*i+4]=segments[3*i+1];
1104- nriftsegs++;
1105- }
1106- }
1107-
1108- /*Compress riftsegments_uncompressed:*/
1109- riftsegments=xNew<int>(nriftsegs*4);
1110- counter=0;
1111- for (i=0;i<nsegs;i++){
1112- if (riftsegments_uncompressed[5*i+0]){
1113- riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
1114- riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
1115- riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
1116- riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
1117- counter++;
1118- }
1119- }
1120- xDelete<int>(riftsegments_uncompressed);
1121-
1122- /*Assign output pointers: */
1123- *priftsegments=riftsegments;
1124- *pnriftsegs=nriftsegs;
1125-}/*}}}*/
1126-int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){/*{{{*/
1127-
1128- int noerr=1;
1129- int k,l,counter;
1130- int newel;
1131-
1132- int* GridElements=NULL;
1133- int NumGridElements;
1134-
1135- /*Output: */
1136- int NumGridElementListOnOneSideOfRift;
1137- int* GridElementListOnOneSideOfRift=NULL;
1138-
1139- /*Build a list of all the elements connected to this node: */
1140- GridElementsList(&GridElements,&NumGridElements,node,index,nel);
1141-
1142- /*Figure out the list of elements that are on the same side of the rift. To do so, we start from one
1143- * side of the rift and keep rotating in the same direction:*/
1144- GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
1145- //bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
1146- GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there
1147- for a rotation direction, we 'll take it out when we are
1148- done rotating*/
1149- GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
1150- counter=1;
1151- for (;;){
1152- /*Find neighbour of element GridElementListOnOneSideOfRift[counter], not
1153- * equal to GridElementListOnOneSideOfRift[counter-1]*/
1154- for (k=0;k<NumGridElements;k++){
1155- if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
1156- /*Verify this element is not already in our list of element on the same side of the rift: */
1157- newel=1;
1158- for (l=0;l<=counter;l++){
1159- if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
1160- newel=0;
1161- break;
1162- }
1163- }
1164- if (newel){
1165- counter++;
1166- GridElementListOnOneSideOfRift[counter]=GridElements[k];
1167- if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
1168- break;
1169- }
1170- k=-1;
1171- }
1172- }
1173- }
1174- /*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
1175- NumGridElementListOnOneSideOfRift=counter;
1176- for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
1177- GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
1178- }
1179- break;
1180- }// for (;;)
1181-
1182- /*Free ressources: */
1183- xDelete<int>(GridElements);
1184- /*Assign output pointers: */
1185- *pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
1186- *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
1187- return noerr;
1188-}/*}}}*/
1189-int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){/*{{{*/
1190-
1191- int noerr=1;
1192- int i,j,k;
1193- int el1,el2;
1194-
1195- int *segments = NULL;
1196- int *segmentmarkerlist = NULL;
1197- int nsegs;
1198-
1199- /*Recover input: */
1200- segments = *psegments;
1201- segmentmarkerlist = *psegmentmarkerlist;
1202- nsegs = *pnsegs;
1203-
1204- /*Reallocate segments: */
1205- segments =xReNew<int>(segments, nsegs*3,(nsegs+nriftsegs)*3);
1206- segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
1207-
1208- /*First, update the existing segments to the new nodes :*/
1209- for (i=0;i<nriftsegs;i++){
1210- el1=riftsegments[4*i+0];
1211- el2=riftsegments[4*i+1];
1212- for (j=0;j<nsegs;j++){
1213- if (segments[3*j+2]==(el1+1)){
1214- /*segment j is the same as rift segment i.Let's update segments[j][:] using element el1 and the corresponding rift segment.
1215- *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
1216- *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
1217- for (k=0;k<3;k++){
1218- if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
1219- *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
1220- break;
1221- }
1222- }
1223- for (k=0;k<3;k++){
1224- if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
1225- *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
1226- break;
1227- }
1228- }
1229- /*Deal with el2: */
1230- *(segments+3*(nsegs+i)+2)=el2+1;
1231- *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
1232- for (k=0;k<3;k++){
1233- if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
1234- *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
1235- break;
1236- }
1237- }
1238- for (k=0;k<3;k++){
1239- if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
1240- *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
1241- break;
1242- }
1243- }
1244- }
1245- if (*(segments+3*j+2)==(el2+1)){
1246- /*segment j is the same as rift segment i.*/
1247- /*Let's update segments[j][:] using element el2 and the corresponding rift segment: */
1248- for (k=0;k<3;k++){
1249- if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
1250- *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
1251- break;
1252- }
1253- }
1254- for (k=0;k<3;k++){
1255- if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
1256- *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
1257- break;
1258- }
1259- }
1260- /*Deal with el1: */
1261- *(segments+3*(nsegs+i)+2)=el1+1;
1262- *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
1263- for (k=0;k<3;k++){
1264- if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
1265- *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
1266- break;
1267- }
1268- }
1269- for (k=0;k<3;k++){
1270- if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
1271- *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
1272- break;
1273- }
1274- }
1275- }
1276- }
1277- }
1278- nsegs+=nriftsegs;
1279-
1280- /*Assign output pointers: */
1281- *psegments=segments;
1282- *psegmentmarkerlist=segmentmarkerlist;
1283- *pnsegs=nsegs;
1284-
1285- return noerr;
1286-}/*}}}*/
1287-int FindElement(int A,int B,int* index,int nel){/*{{{*/
1288-
1289- int el=-1;
1290- for (int n=0;n<nel;n++){
1291- if(((index[3*n+0]==A) || (index[3*n+1]==A) || (index[3*n+2]==A)) && ((index[3*n+0]==B) || (index[3*n+1]==B) || (index[3*n+2]==B))){
1292- el=n;
1293- break;
1294- }
1295- }
1296- return el;
1297-}/*}}}*/
1298-int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){/*{{{*/
1299-
1300- /*Using segment markers, wring out the rift segments from the segments. Rift markers are
1301- *of the form 2+i where i=0 to number of rifts */
1302-
1303- int noerr=1;
1304- int i,j,counter;
1305-
1306- /*input: */
1307- int *segments = NULL;
1308- int *segmentmarkerlist = NULL;
1309- int numsegs;
1310-
1311- /*output: */
1312- int new_numsegs;
1313- int *riftsnumsegs = NULL;
1314- int **riftssegments = NULL;
1315- int *new_segments = NULL;
1316- int *new_segmentmarkers = NULL;
1317-
1318- /*intermediary: */
1319- int* riftsegment=NULL;
1320-
1321- /*Recover input arguments: */
1322- segments = *psegments;
1323- numsegs = *pnumsegs;
1324- segmentmarkerlist = *psegmentmarkerlist;
1325-
1326- /*First, figure out how many segments will be left in 'segments': */
1327- counter=0;
1328- for (i=0;i<numsegs;i++){
1329- if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
1330- }
1331- /*Allocate new segments: */
1332- new_numsegs=counter;
1333- new_segments=xNew<int>(new_numsegs*3);
1334- new_segmentmarkers=xNew<int>(new_numsegs);
1335-
1336- /*Copy new segments info : */
1337- counter=0;
1338- for (i=0;i<numsegs;i++){
1339- if (segmentmarkerlist[i]==1){
1340- new_segments[3*counter+0]=segments[3*i+0];
1341- new_segments[3*counter+1]=segments[3*i+1];
1342- new_segments[3*counter+2]=segments[3*i+2];
1343- new_segmentmarkers[counter]=segmentmarkerlist[i];
1344- counter++;
1345- }
1346- }
1347-
1348- /*Now deal with rift segments: */
1349- riftsnumsegs=xNew<int>(numrifts);
1350- riftssegments=xNew<int*>(numrifts);
1351- for (i=0;i<numrifts;i++){
1352- /*Figure out how many segments for rift i: */
1353- counter=0;
1354- for (j=0;j<numsegs;j++){
1355- if (segmentmarkerlist[j]==2+i)counter++;
1356- }
1357- riftsnumsegs[i]=counter;
1358- riftsegment=xNew<int>(counter*3);
1359- /*Copy new segments info :*/
1360- counter=0;
1361- for (j=0;j<numsegs;j++){
1362- if (segmentmarkerlist[j]==(2+i)){
1363- riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods+1);
1364- riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods+1);
1365- riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel+1);
1366- counter++;
1367- }
1368- }
1369- *(riftssegments+i)=riftsegment;
1370- }
1371-
1372- /*Free ressources: */
1373- xDelete<int>(segments);
1374-
1375- /*Assign output pointers: */
1376- *psegments=new_segments;
1377- *psegmentmarkerlist=new_segmentmarkers;
1378- *pnumsegs=new_numsegs;
1379- *pnumrifts=numrifts;
1380- *priftssegments=riftssegments;
1381- *priftsnumsegs=riftsnumsegs;
1382- return noerr;
1383-}/*}}}*/
1384-int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){/*{{{*/
1385-
1386- int noerr=1;
1387- int i,j,k;
1388-
1389- /*output: */
1390- int *riftsnumpairs = NULL;
1391- int **riftspairs = NULL;
1392-
1393- /*intermediary :*/
1394- int numsegs;
1395- int* segments=NULL;
1396- int* pairs=NULL;
1397- int node1,node2,node3,node4;
1398-
1399- riftsnumpairs=xNew<int>(numrifts);
1400- riftspairs=xNew<int*>(numrifts);
1401- for (i=0;i<numrifts;i++){
1402- segments=riftssegments[i];
1403- numsegs =riftsnumsegments[i];
1404- riftsnumpairs[i]=numsegs;
1405- pairs=xNew<int>(2*numsegs);
1406- for (j=0;j<numsegs;j++){
1407- pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
1408- node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
1409- /*Find element facing on other side of rift: */
1410- for (k=0;k<numsegs;k++){
1411- if (k==j)continue;
1412- node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
1413- /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
1414- if ( ((x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]))
1415- || ((x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1])) ){
1416- /*We found the corresponding element: */
1417- pairs[2*j+1]=segments[3*k+2];
1418- break;
1419- }
1420- }
1421- }
1422- riftspairs[i]=pairs;
1423- }
1424-
1425- /*Assign output pointers: */
1426- *priftsnumpairs=riftsnumpairs;
1427- *priftspairs=riftspairs;
1428- return noerr;
1429-}/*}}}*/
1430-int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){/*{{{*/
1431-
1432- int i;
1433- int noerr=1;
1434-
1435- /*output: */
1436- int riftflag=0;
1437- int numrifts=0;
1438-
1439- int maxmark=1; //default marker for regular segments
1440-
1441- /*Any marker >=2 indicates a certain rift: */
1442- numrifts=0;
1443- for (i=0;i<nsegs;i++){
1444- if (segmentmarkerlist[i]>maxmark){
1445- numrifts++;
1446- maxmark=segmentmarkerlist[i];
1447- }
1448- }
1449- if(numrifts)riftflag=1;
1450-
1451- /*Assign output pointers:*/
1452- *priftflag=riftflag;
1453- *pnumrifts=numrifts;
1454- return noerr;
1455-}/*}}}*/
1456-int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){/*{{{*/
1457-
1458- int noerr=1;
1459- int i,j,k,counter;
1460-
1461- /*intermediary: */
1462- int *riftsegments = NULL;
1463- int *riftpairs = NULL;
1464- int numsegs;
1465-
1466- /*ordering and copy: */
1467- int *order = NULL;
1468- int *riftsegments_copy = NULL;
1469- int *riftpairs_copy = NULL;
1470-
1471- /*node and element manipulation: */
1472- int node1,node2,node3,node4,temp_node,tip1,tip2,node;
1473- int el2;
1474- int already_ordered=0;
1475-
1476- /*output: */
1477- int* riftstips=NULL;
1478-
1479- /*Allocate byproduct of this routine, riftstips: */
1480- riftstips=xNew<int>(numrifts*2);
1481-
1482- /*Go through all rifts: */
1483- for (i=0;i<numrifts;i++){
1484- riftsegments = riftssegments[i];
1485- riftpairs = riftspairs[i];
1486- numsegs = riftsnumsegments[i];
1487-
1488- /*Allocate copy of riftsegments and riftpairs,
1489- *as well as ordering vector: */
1490- riftsegments_copy=xNew<int>(numsegs*3);
1491- riftpairs_copy=xNew<int>(numsegs*2);
1492- order=xNew<int>(numsegs);
1493-
1494- /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
1495- tip1=-1;
1496- tip2=-1;
1497-
1498- for (j=0;j<numsegs;j++){
1499- el2=*(riftpairs+2*j+1);
1500- node1=*(riftsegments+3*j+0);
1501- node2=*(riftsegments+3*j+1);
1502- /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
1503- *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
1504- for (k=0;k<numsegs;k++){
1505- if (*(riftsegments+3*k+2)==el2){
1506- node3=*(riftsegments+3*k+0);
1507- node4=*(riftsegments+3*k+1);
1508- break;
1509- }
1510- }
1511- /* Make sure node3 faces node1 and node4 faces node2: */
1512- _assert_(node1<nods+1 && node4<nods+1);
1513- _assert_(node1>0 && node4>0);
1514- if ((x[node1-1]==x[node4-1]) && (y[node1-1]==y[node4-1])){
1515- /*Swap node3 and node4:*/
1516- temp_node=node3;
1517- node3=node4;
1518- node4=temp_node;
1519- }
1520-
1521- /*Figure out if a tip is on this element: */
1522- if (node3==node1){
1523- /*node1 is a tip*/
1524- if (tip1==-1) {
1525- tip1=node1;
1526- continue;
1527- }
1528- if ((tip2==-1) && (node1!=tip1)){
1529- tip2=node1;
1530- break;
1531- }
1532- }
1533-
1534- if (node4==node2){
1535- /*node2 is a tip*/
1536- if (tip1==-1){
1537- tip1=node2;
1538- continue;
1539- }
1540- if ((tip2==-1) && (node2!=tip1)){
1541- tip2=node2;
1542- break;
1543- }
1544- }
1545- }
1546-
1547- /*Record tips in riftstips: */
1548- *(riftstips+2*i+0)=tip1;
1549- *(riftstips+2*i+1)=tip2;
1550-
1551- /*We have the two tips for this rift. Go from tip1 to tip2, and figure out the order in which segments are sequential.
1552- *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
1553- node=tip1;
1554- for (counter=0;counter<numsegs;counter++){
1555- for (j=0;j<numsegs;j++){
1556- node1=*(riftsegments+3*j+0);
1557- node2=*(riftsegments+3*j+1);
1558-
1559- if ((node1==node) || (node2==node)){
1560- /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
1561- already_ordered=0;
1562- for (k=0;k<counter;k++){
1563- if(order[k]==j){
1564- already_ordered=1;
1565- break;
1566- }
1567- }
1568- if (!already_ordered){
1569- order[counter]=j;
1570- if(node1==node){
1571- node=node2;
1572- }
1573- else if(node2==node){
1574- node=node1;
1575- }
1576- break;
1577- }
1578- }
1579- }
1580- }
1581-
1582- /*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
1583- for (j=0;j<numsegs;j++){
1584- _assert_(order[j]<numsegs);
1585- *(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
1586- *(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
1587- *(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
1588- *(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
1589- *(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
1590- }
1591-
1592- for (j=0;j<numsegs;j++){
1593- *(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
1594- *(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
1595- *(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
1596- *(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
1597- *(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
1598- }
1599-
1600- xDelete<int>(order);
1601- xDelete<int>(riftsegments_copy);
1602- xDelete<int>(riftpairs_copy);
1603-
1604- }
1605-
1606- /*Assign output pointer:*/
1607- *priftstips=riftstips;
1608- return noerr;
1609-}/*}}}*/
1610-int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,/*{{{*/
1611- int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
1612-
1613- int noerr=1;
1614- int i,j,k,k0;
1615-
1616- double el1,el2,node1,node2,node3,node4;
1617- double temp_node;
1618-
1619- /*output: */
1620- double **riftspenaltypairs = NULL;
1621- double *riftpenaltypairs = NULL;
1622- int *riftsnumpenaltypairs = NULL;
1623-
1624- /*intermediary: */
1625- int numsegs;
1626- int* riftsegments=NULL;
1627- int* riftpairs=NULL;
1628- int counter;
1629- double normal[2];
1630- double length;
1631- int k1,k2;
1632-
1633- /*Allocate: */
1634- riftspenaltypairs=xNew<double*>(numrifts);
1635- riftsnumpenaltypairs=xNew<int>(numrifts);
1636-
1637- for(i=0;i<numrifts;i++){
1638- numsegs=riftsnumsegs[i];
1639- riftsegments=riftssegments[i];
1640- riftpairs=riftspairs[i];
1641-
1642- /*allocate riftpenaltypairs, and riftnumpenaltypairs: */
1643- if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
1644-
1645- /*Go through only one flank of the rifts, not counting the tips: */
1646- counter=0;
1647- for(j=0;j<(numsegs/2);j++){
1648- el1=*(riftpairs+2*j+0);
1649- el2=*(riftpairs+2*j+1);
1650- node1=*(riftsegments+3*j+0);
1651- node2=*(riftsegments+3*j+1);
1652- /*Find segment index to recover node3 and node4, facing node1 and node2: */
1653- k0=-1;
1654- for(k=0;k<numsegs;k++){
1655- if(*(riftsegments+3*k+2)==el2){
1656- k0=k;
1657- break;
1658- }
1659- }
1660- node3=*(riftsegments+3*k0+0);
1661- node4=*(riftsegments+3*k0+1);
1662-
1663- /* Make sure node3 faces node1 and node4 faces node2: */
1664- if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
1665- /*Swap node3 and node4:*/
1666- temp_node=node3;
1667- node3=node4;
1668- node4=temp_node;
1669- }
1670- /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
1671- *this segment, and its length: */
1672- normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
1673- normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
1674- length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
1675-
1676- /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
1677- * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4,
1678- * only once. We'll add the normals and the lengths : */
1679-
1680- if(node1!=node3){ //exclude tips from loads
1681- k1=-1;
1682- for(k=0;k<counter;k++){
1683- if( (*(riftpenaltypairs+k*7+0))==node1){
1684- k1=k;
1685- break;
1686- }
1687- }
1688- if(k1==-1){
1689- *(riftpenaltypairs+counter*7+0)=node1;
1690- *(riftpenaltypairs+counter*7+1)=node3;
1691- *(riftpenaltypairs+counter*7+2)=el1;
1692- *(riftpenaltypairs+counter*7+3)=el2;
1693- *(riftpenaltypairs+counter*7+4)=normal[0];
1694- *(riftpenaltypairs+counter*7+5)=normal[1];
1695- *(riftpenaltypairs+counter*7+6)=length/2;
1696- counter++;
1697- }
1698- else{
1699- *(riftpenaltypairs+k1*7+4)+=normal[0];
1700- *(riftpenaltypairs+k1*7+5)+=normal[1];
1701- *(riftpenaltypairs+k1*7+6)+=length/2;
1702- }
1703- }
1704- if(node2!=node4){
1705- k2=-1;
1706- for(k=0;k<counter;k++){
1707- if( (*(riftpenaltypairs+k*7+0))==node2){
1708- k2=k;
1709- break;
1710- }
1711- }
1712- if(k2==-1){
1713- *(riftpenaltypairs+counter*7+0)=node2;
1714- *(riftpenaltypairs+counter*7+1)=node4;
1715- *(riftpenaltypairs+counter*7+2)=el1;
1716- *(riftpenaltypairs+counter*7+3)=el2;
1717- *(riftpenaltypairs+counter*7+4)=normal[0];
1718- *(riftpenaltypairs+counter*7+5)=normal[1];
1719- *(riftpenaltypairs+counter*7+6)=length/2;
1720- counter++;
1721- }
1722- else{
1723- *(riftpenaltypairs+k2*7+4)+=normal[0];
1724- *(riftpenaltypairs+k2*7+5)+=normal[1];
1725- *(riftpenaltypairs+k2*7+6)+=length/2;
1726- }
1727- }
1728- }
1729- /*Renormalize normals: */
1730- for(j=0;j<counter;j++){
1731- double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
1732- *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
1733- *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
1734- }
1735-
1736- riftspenaltypairs[i]=riftpenaltypairs;
1737- riftsnumpenaltypairs[i]=(numsegs/2-1);
1738- }
1739-
1740- /*Assign output pointers: */
1741- *priftspenaltypairs=riftspenaltypairs;
1742- *priftsnumpenaltypairs=riftsnumpenaltypairs;
1743- return noerr;
1744-}/*}}}*/
1745-int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
1746-
1747- int noerr=1;
1748- int i,j,k;
1749- int node1,node2,node3;
1750- int el;
1751- double pair[2];
1752- int pair_count=0;
1753- int triple=0;
1754-
1755- /*Recover input: */
1756- int *index = *pindex;
1757- int nel = *pnel;
1758- double *x = *px;
1759- double *y = *py;
1760- int nods = *pnods;
1761-
1762- for (i=0;i<num_seg;i++){
1763- node1=*(segments+3*i+0);
1764- node2=*(segments+3*i+1);
1765- /*Find all elements connected to [node1 node2]: */
1766- pair_count=0;
1767- for (j=0;j<nel;j++){
1768- if (*(index+3*j+0)==node1){
1769- if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
1770- pair[pair_count]=j;
1771- pair_count++;
1772- }
1773- }
1774- if (*(index+3*j+1)==node1){
1775- if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
1776- pair[pair_count]=j;
1777- pair_count++;
1778- }
1779- }
1780- if (*(index+3*j+2)==node1){
1781- if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
1782- pair[pair_count]=j;
1783- pair_count++;
1784- }
1785- }
1786- }
1787- /*Ok, we have pair_count elements connected to this segment. For each of these elements,
1788- *figure out if the third node also belongs to a segment: */
1789- if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to 2 elements
1790- continue;
1791- }
1792- else{
1793- for (j=0;j<pair_count;j++){
1794- el=(int)pair[j];
1795- triple=0;
1796- /*First find node3: */
1797- if (*(index+3*el+0)==node1){
1798- if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
1799- else node3=*(index+3*el+1);
1800- }
1801- if (*(index+3*el+1)==node1){
1802- if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
1803- else node3=*(index+3*el+0);
1804- }
1805- if (*(index+3*el+2)==node1){
1806- if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
1807- else node3=*(index+3*el+0);
1808- }
1809- /*Ok, we have node3. Does node3 belong to a segment? : */
1810- for (k=0;k<num_seg;k++){
1811- if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
1812- triple=1;
1813- break;
1814- }
1815- }
1816- if(triple==1){
1817- /*el is a corner element: we need to split it in 3 triangles: */
1818- x=xReNew<double>(x,nods,nods+1);
1819- y=xReNew<double>(y,nods,nods+1);
1820- x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
1821- y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
1822- index=xReNew<int>(index,nel*3,(nel+2*3));
1823- /*First, reassign element el: */
1824- *(index+3*el+0)=node1;
1825- *(index+3*el+1)=node2;
1826- *(index+3*el+2)=nods+1;
1827- /*Other two elements: */
1828- *(index+3*nel+0)=node2;
1829- *(index+3*nel+1)=node3;
1830- *(index+3*nel+2)=nods+1;
1831-
1832- *(index+3*(nel+1)+0)=node3;
1833- *(index+3*(nel+1)+1)=node1;
1834- *(index+3*(nel+1)+2)=nods+1;
1835- /*we need to change the segment elements corresponding to el: */
1836- for (k=0;k<num_seg;k++){
1837- if (*(segments+3*k+2)==(el+1)){
1838- if ( ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node2)) || ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node1))) *(segments+3*k+2)=el+1;
1839- if ( ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node3)) || ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node2))) *(segments+3*k+2)=nel+1;
1840- if ( ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node1)) || ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node3))) *(segments+3*k+2)=nel+2;
1841- }
1842- }
1843-
1844- nods=nods+1;
1845- nel=nel+2;
1846- i=0;
1847- break;
1848- }
1849- } //for (j=0;j<pair_count;j++)
1850- }
1851- }// for (i=0;i<num_seg;i++)
1852-
1853- /*Assign output pointers: */
1854- *pindex=index;
1855- *pnel=nel;
1856- *px=x;
1857- *py=y;
1858- *pnods=nods;
1859- return noerr;
1860-}/*}}}*/
1861Index: ../trunk-jpl/src/c/shared/TriMesh/trimesh.h
1862===================================================================
1863--- ../trunk-jpl/src/c/shared/TriMesh/trimesh.h (revision 22497)
1864+++ ../trunk-jpl/src/c/shared/TriMesh/trimesh.h (nonexistent)
1865@@ -1,33 +0,0 @@
1866-/*!\file: trimesh.h
1867- * \brief
1868- */
1869-
1870-#ifndef _SHARED_TRIMESH_H
1871-#define _SHARED_TRIMESH_H
1872-
1873-#include <stdio.h>
1874-#include <math.h>
1875-
1876-//#define REAL double //took it out because it may conflict with stdlib.h defines. put back if necessary
1877-int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel);
1878-int OrderSegments(int** psegments,int nseg, int* index,int nel);
1879-int GridInsideHole(double* px0,double* py0,int n,double* x,double* y);
1880-int FindElement(int A,int B,int* index,int nel);
1881-int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist);
1882-int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
1883-int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
1884-int IsNeighbor(int el1,int el2,int* index);
1885-int IsOnRift(int el,int nriftsegs,int* riftsegments);
1886-void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index, int nsegs,int* segments);
1887-int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift,int segmentnumber, int nriftsegs,int* riftsegments, int node,int* index,int nel);
1888-int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
1889-int FindElement(double A,double B,int* index,int nel);
1890-int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
1891-int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels);
1892-int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels);
1893-int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,
1894- int* riftsnumsegments,int** riftspairs,int* riftstips,double* x,double* y);
1895-int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg);
1896-int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y);
1897-
1898-#endif /* _SHARED_TRIMESH_H */
1899Index: ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp
1900===================================================================
1901--- ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp (revision 22497)
1902+++ ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp (nonexistent)
1903@@ -1,46 +0,0 @@
1904-/*
1905- * OrderSegments.c:
1906- * reorder segments so that their normals point outside the domain outline.
1907- */
1908-#include "./trimesh.h"
1909-
1910-int OrderSegments(int** psegments,int nseg,int* index,int nel){
1911-
1912- /*vertex indices: */
1913- int A,B;
1914-
1915- /*element index*/
1916- int el;
1917-
1918- /*Recover segments: */
1919- int* segments=*psegments;
1920-
1921- for(int i=0;i<nseg;i++){
1922- A=segments[3*i+0];
1923- B=segments[3*i+1];
1924- el=segments[3*i+2]-1; //after AssociateSegmentToElement, el was a matlab index, we need the c index now.
1925-
1926- if (index[3*el+0]==A){
1927- if (index[3*el+2]==B){
1928- segments[3*i+0]=B;
1929- segments[3*i+1]=A;
1930- }
1931- }
1932- else if (index[3*el+1]==A){
1933- if (index[3*el+0]==B){
1934- segments[3*i+0]=B;
1935- segments[3*i+1]=A;
1936- }
1937- }
1938- else{
1939- if (index[3*el+1]==B){
1940- segments[3*i+0]=B;
1941- segments[3*i+1]=A;
1942- }
1943- }
1944- }
1945-
1946- /*Assign output pointers: */
1947- *psegments=segments;
1948- return 1;
1949-}
1950Index: ../trunk-jpl/src/c/shared/TriMesh
1951===================================================================
1952--- ../trunk-jpl/src/c/shared/TriMesh (revision 22497)
1953+++ ../trunk-jpl/src/c/shared/TriMesh (nonexistent)
1954
1955Property changes on: ../trunk-jpl/src/c/shared/TriMesh
1956___________________________________________________________________
1957Deleted: svn:ignore
1958## -1,2 +0,0 ##
1959-.deps
1960-.dirstamp
1961Index: ../trunk-jpl/src/c/shared/shared.h
1962===================================================================
1963--- ../trunk-jpl/src/c/shared/shared.h (revision 22497)
1964+++ ../trunk-jpl/src/c/shared/shared.h (revision 22498)
1965@@ -18,7 +18,7 @@
1966 #include "./Sorting/sorting.h"
1967 #include "./String/sharedstring.h"
1968 #include "./Threads/issm_threads.h"
1969-#include "./TriMesh/trimesh.h"
1970+#include "./Triangle/triangle.h"
1971 #include "./LatLong/latlong.h"
1972
1973 #endif
1974Index: ../trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp
1975===================================================================
1976--- ../trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp (nonexistent)
1977+++ ../trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp (revision 22498)
1978@@ -0,0 +1,46 @@
1979+/*
1980+ * OrderSegments.c:
1981+ * reorder segments so that their normals point outside the domain outline.
1982+ */
1983+#include "./trimesh.h"
1984+
1985+int OrderSegments(int** psegments,int nseg,int* index,int nel){
1986+
1987+ /*vertex indices: */
1988+ int A,B;
1989+
1990+ /*element index*/
1991+ int el;
1992+
1993+ /*Recover segments: */
1994+ int* segments=*psegments;
1995+
1996+ for(int i=0;i<nseg;i++){
1997+ A=segments[3*i+0];
1998+ B=segments[3*i+1];
1999+ el=segments[3*i+2]-1; //after AssociateSegmentToElement, el was a matlab index, we need the c index now.
2000+
2001+ if (index[3*el+0]==A){
2002+ if (index[3*el+2]==B){
2003+ segments[3*i+0]=B;
2004+ segments[3*i+1]=A;
2005+ }
2006+ }
2007+ else if (index[3*el+1]==A){
2008+ if (index[3*el+0]==B){
2009+ segments[3*i+0]=B;
2010+ segments[3*i+1]=A;
2011+ }
2012+ }
2013+ else{
2014+ if (index[3*el+1]==B){
2015+ segments[3*i+0]=B;
2016+ segments[3*i+1]=A;
2017+ }
2018+ }
2019+ }
2020+
2021+ /*Assign output pointers: */
2022+ *psegments=segments;
2023+ return 1;
2024+}
2025Index: ../trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp
2026===================================================================
2027--- ../trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp (nonexistent)
2028+++ ../trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp (revision 22498)
2029@@ -0,0 +1,912 @@
2030+/*
2031+ * TriangleUtils: mesh manipulation routines:
2032+ */
2033+
2034+#include <stdio.h>
2035+
2036+#include "./triangle.h"
2037+#include "../Exceptions/exceptions.h"
2038+#include "../MemOps/MemOps.h"
2039+
2040+#define RIFTPENALTYPAIRSWIDTH 8
2041+int IsGridOnRift(int* riftsegments, int nriftsegs, int node){/*{{{*/
2042+
2043+ /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
2044+ *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).*/
2045+
2046+ int i;
2047+ int j;
2048+ int count;
2049+
2050+ count=0;
2051+ for (i=0;i<nriftsegs;i++){
2052+ for (j=0;j<2;j++){
2053+ if ((*(riftsegments+4*i+2+j))==node) count++;
2054+ }
2055+ }
2056+ if (count==2){
2057+ return 1;
2058+ }
2059+ else{
2060+ return 0;
2061+ }
2062+}/*}}}*/
2063+int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){/*{{{*/
2064+
2065+ /*From a node, recover all the elements that are connected to it: */
2066+ int i,j;
2067+ int noerr=1;
2068+
2069+ int max_number_elements=12;
2070+ int current_size;
2071+ int NumGridElements;
2072+ int* GridElements=NULL;
2073+ int* GridElementsRealloc=NULL;
2074+
2075+ /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own
2076+ * the node. We start by allocating GridElements with that size, and realloc
2077+ * more if needed.*/
2078+
2079+ current_size=max_number_elements;
2080+ NumGridElements=0;
2081+ GridElements=xNew<int>(max_number_elements);
2082+
2083+ for (i=0;i<nel;i++){
2084+ for (j=0;j<3;j++){
2085+ if (index[3*i+j]==node){
2086+ if (NumGridElements<=(current_size-1)){
2087+ GridElements[NumGridElements]=i;
2088+ NumGridElements++;
2089+ break;
2090+ }
2091+ else{
2092+ /*Reallocate another max_number_elements slots in the GridElements: */
2093+ GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
2094+ if (!GridElementsRealloc){
2095+ noerr=0;
2096+ goto cleanup_and_return;
2097+ }
2098+ current_size+=max_number_elements;
2099+ GridElements=GridElementsRealloc;
2100+ GridElements[NumGridElements]=i;
2101+ NumGridElements++;
2102+ break;
2103+ }
2104+ }
2105+ }
2106+ }
2107+ cleanup_and_return:
2108+ if(!noerr){
2109+ xDelete<int>(GridElements);
2110+ }
2111+ /*Allocate return pointers: */
2112+ *pGridElements=GridElements;
2113+ *pNumGridElements=NumGridElements;
2114+ return noerr;
2115+}/*}}}*/
2116+int IsNeighbor(int el1,int el2,int* index){/*{{{*/
2117+ /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
2118+ int i,j;
2119+ int count=0;
2120+ for (i=0;i<3;i++){
2121+ for (j=0;j<3;j++){
2122+ if (index[3*el1+i]==index[3*el2+j])count++;
2123+ }
2124+ }
2125+ if (count==2){
2126+ return 1;
2127+ }
2128+ else{
2129+ return 0;
2130+ }
2131+}/*}}}*/
2132+int IsOnRift(int el,int nriftsegs,int* riftsegments){/*{{{*/
2133+ /*From a list of elements segments, figure out if el belongs to it: */
2134+ int i;
2135+ for (i=0;i<nriftsegs;i++){
2136+ if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
2137+ return 1;
2138+ }
2139+ }
2140+ return 0;
2141+}/*}}}*/
2142+void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){/*{{{*/
2143+
2144+ int i,counter;
2145+ int el,el2;
2146+
2147+ int nriftsegs;
2148+ int* riftsegments=NULL;
2149+ int* riftsegments_uncompressed=NULL;
2150+ int element_nodes[3];
2151+
2152+ /*Allocate segmentflags: */
2153+ riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
2154+
2155+ /*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary
2156+ *or a hole: */
2157+ nriftsegs=0;
2158+ for (i=0;i<nsegs;i++){
2159+ el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
2160+ /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
2161+ *then proceed to find another element that owns the segment. If we don't find it, we know
2162+ *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
2163+ element_nodes[0]=*(index+3*el+0);
2164+ element_nodes[1]=*(index+3*el+1);
2165+ element_nodes[2]=*(index+3*el+2);
2166+
2167+ index[3*el+0]=-1;
2168+ index[3*el+1]=-1;
2169+ index[3*el+2]=-1;
2170+
2171+ el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
2172+
2173+ /*Restore index: */
2174+ index[3*el+0]=element_nodes[0];
2175+ index[3*el+1]=element_nodes[1];
2176+ index[3*el+2]=element_nodes[2];
2177+
2178+ if (el2!=-1){
2179+ /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
2180+ riftsegments_uncompressed[5*i+0]=1;
2181+ riftsegments_uncompressed[5*i+1]=el;
2182+ riftsegments_uncompressed[5*i+2]=el2;
2183+ riftsegments_uncompressed[5*i+3]=segments[3*i+0];
2184+ riftsegments_uncompressed[5*i+4]=segments[3*i+1];
2185+ nriftsegs++;
2186+ }
2187+ }
2188+
2189+ /*Compress riftsegments_uncompressed:*/
2190+ riftsegments=xNew<int>(nriftsegs*4);
2191+ counter=0;
2192+ for (i=0;i<nsegs;i++){
2193+ if (riftsegments_uncompressed[5*i+0]){
2194+ riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
2195+ riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
2196+ riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
2197+ riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
2198+ counter++;
2199+ }
2200+ }
2201+ xDelete<int>(riftsegments_uncompressed);
2202+
2203+ /*Assign output pointers: */
2204+ *priftsegments=riftsegments;
2205+ *pnriftsegs=nriftsegs;
2206+}/*}}}*/
2207+int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){/*{{{*/
2208+
2209+ int noerr=1;
2210+ int k,l,counter;
2211+ int newel;
2212+
2213+ int* GridElements=NULL;
2214+ int NumGridElements;
2215+
2216+ /*Output: */
2217+ int NumGridElementListOnOneSideOfRift;
2218+ int* GridElementListOnOneSideOfRift=NULL;
2219+
2220+ /*Build a list of all the elements connected to this node: */
2221+ GridElementsList(&GridElements,&NumGridElements,node,index,nel);
2222+
2223+ /*Figure out the list of elements that are on the same side of the rift. To do so, we start from one
2224+ * side of the rift and keep rotating in the same direction:*/
2225+ GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
2226+ //bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
2227+ GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there
2228+ for a rotation direction, we 'll take it out when we are
2229+ done rotating*/
2230+ GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
2231+ counter=1;
2232+ for (;;){
2233+ /*Find neighbour of element GridElementListOnOneSideOfRift[counter], not
2234+ * equal to GridElementListOnOneSideOfRift[counter-1]*/
2235+ for (k=0;k<NumGridElements;k++){
2236+ if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
2237+ /*Verify this element is not already in our list of element on the same side of the rift: */
2238+ newel=1;
2239+ for (l=0;l<=counter;l++){
2240+ if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
2241+ newel=0;
2242+ break;
2243+ }
2244+ }
2245+ if (newel){
2246+ counter++;
2247+ GridElementListOnOneSideOfRift[counter]=GridElements[k];
2248+ if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
2249+ break;
2250+ }
2251+ k=-1;
2252+ }
2253+ }
2254+ }
2255+ /*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
2256+ NumGridElementListOnOneSideOfRift=counter;
2257+ for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
2258+ GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
2259+ }
2260+ break;
2261+ }// for (;;)
2262+
2263+ /*Free ressources: */
2264+ xDelete<int>(GridElements);
2265+ /*Assign output pointers: */
2266+ *pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
2267+ *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
2268+ return noerr;
2269+}/*}}}*/
2270+int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){/*{{{*/
2271+
2272+ int noerr=1;
2273+ int i,j,k;
2274+ int el1,el2;
2275+
2276+ int *segments = NULL;
2277+ int *segmentmarkerlist = NULL;
2278+ int nsegs;
2279+
2280+ /*Recover input: */
2281+ segments = *psegments;
2282+ segmentmarkerlist = *psegmentmarkerlist;
2283+ nsegs = *pnsegs;
2284+
2285+ /*Reallocate segments: */
2286+ segments =xReNew<int>(segments, nsegs*3,(nsegs+nriftsegs)*3);
2287+ segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
2288+
2289+ /*First, update the existing segments to the new nodes :*/
2290+ for (i=0;i<nriftsegs;i++){
2291+ el1=riftsegments[4*i+0];
2292+ el2=riftsegments[4*i+1];
2293+ for (j=0;j<nsegs;j++){
2294+ if (segments[3*j+2]==(el1+1)){
2295+ /*segment j is the same as rift segment i.Let's update segments[j][:] using element el1 and the corresponding rift segment.
2296+ *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
2297+ *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
2298+ for (k=0;k<3;k++){
2299+ if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
2300+ *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
2301+ break;
2302+ }
2303+ }
2304+ for (k=0;k<3;k++){
2305+ if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
2306+ *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
2307+ break;
2308+ }
2309+ }
2310+ /*Deal with el2: */
2311+ *(segments+3*(nsegs+i)+2)=el2+1;
2312+ *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
2313+ for (k=0;k<3;k++){
2314+ if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
2315+ *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
2316+ break;
2317+ }
2318+ }
2319+ for (k=0;k<3;k++){
2320+ if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
2321+ *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
2322+ break;
2323+ }
2324+ }
2325+ }
2326+ if (*(segments+3*j+2)==(el2+1)){
2327+ /*segment j is the same as rift segment i.*/
2328+ /*Let's update segments[j][:] using element el2 and the corresponding rift segment: */
2329+ for (k=0;k<3;k++){
2330+ if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+0)-1])){
2331+ *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
2332+ break;
2333+ }
2334+ }
2335+ for (k=0;k<3;k++){
2336+ if ((x[*(index+el2*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el2*3+k)-1]==y[*(segments+3*j+1)-1])){
2337+ *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
2338+ break;
2339+ }
2340+ }
2341+ /*Deal with el1: */
2342+ *(segments+3*(nsegs+i)+2)=el1+1;
2343+ *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
2344+ for (k=0;k<3;k++){
2345+ if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+0)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+0)-1])){
2346+ *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
2347+ break;
2348+ }
2349+ }
2350+ for (k=0;k<3;k++){
2351+ if ((x[*(index+el1*3+k)-1]==x[*(segments+3*j+1)-1]) && (y[*(index+el1*3+k)-1]==y[*(segments+3*j+1)-1])){
2352+ *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
2353+ break;
2354+ }
2355+ }
2356+ }
2357+ }
2358+ }
2359+ nsegs+=nriftsegs;
2360+
2361+ /*Assign output pointers: */
2362+ *psegments=segments;
2363+ *psegmentmarkerlist=segmentmarkerlist;
2364+ *pnsegs=nsegs;
2365+
2366+ return noerr;
2367+}/*}}}*/
2368+int FindElement(int A,int B,int* index,int nel){/*{{{*/
2369+
2370+ int el=-1;
2371+ for (int n=0;n<nel;n++){
2372+ if(((index[3*n+0]==A) || (index[3*n+1]==A) || (index[3*n+2]==A)) && ((index[3*n+0]==B) || (index[3*n+1]==B) || (index[3*n+2]==B))){
2373+ el=n;
2374+ break;
2375+ }
2376+ }
2377+ return el;
2378+}/*}}}*/
2379+int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){/*{{{*/
2380+
2381+ /*Using segment markers, wring out the rift segments from the segments. Rift markers are
2382+ *of the form 2+i where i=0 to number of rifts */
2383+
2384+ int noerr=1;
2385+ int i,j,counter;
2386+
2387+ /*input: */
2388+ int *segments = NULL;
2389+ int *segmentmarkerlist = NULL;
2390+ int numsegs;
2391+
2392+ /*output: */
2393+ int new_numsegs;
2394+ int *riftsnumsegs = NULL;
2395+ int **riftssegments = NULL;
2396+ int *new_segments = NULL;
2397+ int *new_segmentmarkers = NULL;
2398+
2399+ /*intermediary: */
2400+ int* riftsegment=NULL;
2401+
2402+ /*Recover input arguments: */
2403+ segments = *psegments;
2404+ numsegs = *pnumsegs;
2405+ segmentmarkerlist = *psegmentmarkerlist;
2406+
2407+ /*First, figure out how many segments will be left in 'segments': */
2408+ counter=0;
2409+ for (i=0;i<numsegs;i++){
2410+ if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
2411+ }
2412+ /*Allocate new segments: */
2413+ new_numsegs=counter;
2414+ new_segments=xNew<int>(new_numsegs*3);
2415+ new_segmentmarkers=xNew<int>(new_numsegs);
2416+
2417+ /*Copy new segments info : */
2418+ counter=0;
2419+ for (i=0;i<numsegs;i++){
2420+ if (segmentmarkerlist[i]==1){
2421+ new_segments[3*counter+0]=segments[3*i+0];
2422+ new_segments[3*counter+1]=segments[3*i+1];
2423+ new_segments[3*counter+2]=segments[3*i+2];
2424+ new_segmentmarkers[counter]=segmentmarkerlist[i];
2425+ counter++;
2426+ }
2427+ }
2428+
2429+ /*Now deal with rift segments: */
2430+ riftsnumsegs=xNew<int>(numrifts);
2431+ riftssegments=xNew<int*>(numrifts);
2432+ for (i=0;i<numrifts;i++){
2433+ /*Figure out how many segments for rift i: */
2434+ counter=0;
2435+ for (j=0;j<numsegs;j++){
2436+ if (segmentmarkerlist[j]==2+i)counter++;
2437+ }
2438+ riftsnumsegs[i]=counter;
2439+ riftsegment=xNew<int>(counter*3);
2440+ /*Copy new segments info :*/
2441+ counter=0;
2442+ for (j=0;j<numsegs;j++){
2443+ if (segmentmarkerlist[j]==(2+i)){
2444+ riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods+1);
2445+ riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods+1);
2446+ riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel+1);
2447+ counter++;
2448+ }
2449+ }
2450+ *(riftssegments+i)=riftsegment;
2451+ }
2452+
2453+ /*Free ressources: */
2454+ xDelete<int>(segments);
2455+
2456+ /*Assign output pointers: */
2457+ *psegments=new_segments;
2458+ *psegmentmarkerlist=new_segmentmarkers;
2459+ *pnumsegs=new_numsegs;
2460+ *pnumrifts=numrifts;
2461+ *priftssegments=riftssegments;
2462+ *priftsnumsegs=riftsnumsegs;
2463+ return noerr;
2464+}/*}}}*/
2465+int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){/*{{{*/
2466+
2467+ int noerr=1;
2468+ int i,j,k;
2469+
2470+ /*output: */
2471+ int *riftsnumpairs = NULL;
2472+ int **riftspairs = NULL;
2473+
2474+ /*intermediary :*/
2475+ int numsegs;
2476+ int* segments=NULL;
2477+ int* pairs=NULL;
2478+ int node1,node2,node3,node4;
2479+
2480+ riftsnumpairs=xNew<int>(numrifts);
2481+ riftspairs=xNew<int*>(numrifts);
2482+ for (i=0;i<numrifts;i++){
2483+ segments=riftssegments[i];
2484+ numsegs =riftsnumsegments[i];
2485+ riftsnumpairs[i]=numsegs;
2486+ pairs=xNew<int>(2*numsegs);
2487+ for (j=0;j<numsegs;j++){
2488+ pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
2489+ node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
2490+ /*Find element facing on other side of rift: */
2491+ for (k=0;k<numsegs;k++){
2492+ if (k==j)continue;
2493+ node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
2494+ /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
2495+ if ( ((x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]))
2496+ || ((x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1])) ){
2497+ /*We found the corresponding element: */
2498+ pairs[2*j+1]=segments[3*k+2];
2499+ break;
2500+ }
2501+ }
2502+ }
2503+ riftspairs[i]=pairs;
2504+ }
2505+
2506+ /*Assign output pointers: */
2507+ *priftsnumpairs=riftsnumpairs;
2508+ *priftspairs=riftspairs;
2509+ return noerr;
2510+}/*}}}*/
2511+int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){/*{{{*/
2512+
2513+ int i;
2514+ int noerr=1;
2515+
2516+ /*output: */
2517+ int riftflag=0;
2518+ int numrifts=0;
2519+
2520+ int maxmark=1; //default marker for regular segments
2521+
2522+ /*Any marker >=2 indicates a certain rift: */
2523+ numrifts=0;
2524+ for (i=0;i<nsegs;i++){
2525+ if (segmentmarkerlist[i]>maxmark){
2526+ numrifts++;
2527+ maxmark=segmentmarkerlist[i];
2528+ }
2529+ }
2530+ if(numrifts)riftflag=1;
2531+
2532+ /*Assign output pointers:*/
2533+ *priftflag=riftflag;
2534+ *pnumrifts=numrifts;
2535+ return noerr;
2536+}/*}}}*/
2537+int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){/*{{{*/
2538+
2539+ int noerr=1;
2540+ int i,j,k,counter;
2541+
2542+ /*intermediary: */
2543+ int *riftsegments = NULL;
2544+ int *riftpairs = NULL;
2545+ int numsegs;
2546+
2547+ /*ordering and copy: */
2548+ int *order = NULL;
2549+ int *riftsegments_copy = NULL;
2550+ int *riftpairs_copy = NULL;
2551+
2552+ /*node and element manipulation: */
2553+ int node1,node2,node3,node4,temp_node,tip1,tip2,node;
2554+ int el2;
2555+ int already_ordered=0;
2556+
2557+ /*output: */
2558+ int* riftstips=NULL;
2559+
2560+ /*Allocate byproduct of this routine, riftstips: */
2561+ riftstips=xNew<int>(numrifts*2);
2562+
2563+ /*Go through all rifts: */
2564+ for (i=0;i<numrifts;i++){
2565+ riftsegments = riftssegments[i];
2566+ riftpairs = riftspairs[i];
2567+ numsegs = riftsnumsegments[i];
2568+
2569+ /*Allocate copy of riftsegments and riftpairs,
2570+ *as well as ordering vector: */
2571+ riftsegments_copy=xNew<int>(numsegs*3);
2572+ riftpairs_copy=xNew<int>(numsegs*2);
2573+ order=xNew<int>(numsegs);
2574+
2575+ /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
2576+ tip1=-1;
2577+ tip2=-1;
2578+
2579+ for (j=0;j<numsegs;j++){
2580+ el2=*(riftpairs+2*j+1);
2581+ node1=*(riftsegments+3*j+0);
2582+ node2=*(riftsegments+3*j+1);
2583+ /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
2584+ *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
2585+ for (k=0;k<numsegs;k++){
2586+ if (*(riftsegments+3*k+2)==el2){
2587+ node3=*(riftsegments+3*k+0);
2588+ node4=*(riftsegments+3*k+1);
2589+ break;
2590+ }
2591+ }
2592+ /* Make sure node3 faces node1 and node4 faces node2: */
2593+ _assert_(node1<nods+1 && node4<nods+1);
2594+ _assert_(node1>0 && node4>0);
2595+ if ((x[node1-1]==x[node4-1]) && (y[node1-1]==y[node4-1])){
2596+ /*Swap node3 and node4:*/
2597+ temp_node=node3;
2598+ node3=node4;
2599+ node4=temp_node;
2600+ }
2601+
2602+ /*Figure out if a tip is on this element: */
2603+ if (node3==node1){
2604+ /*node1 is a tip*/
2605+ if (tip1==-1) {
2606+ tip1=node1;
2607+ continue;
2608+ }
2609+ if ((tip2==-1) && (node1!=tip1)){
2610+ tip2=node1;
2611+ break;
2612+ }
2613+ }
2614+
2615+ if (node4==node2){
2616+ /*node2 is a tip*/
2617+ if (tip1==-1){
2618+ tip1=node2;
2619+ continue;
2620+ }
2621+ if ((tip2==-1) && (node2!=tip1)){
2622+ tip2=node2;
2623+ break;
2624+ }
2625+ }
2626+ }
2627+
2628+ /*Record tips in riftstips: */
2629+ *(riftstips+2*i+0)=tip1;
2630+ *(riftstips+2*i+1)=tip2;
2631+
2632+ /*We have the two tips for this rift. Go from tip1 to tip2, and figure out the order in which segments are sequential.
2633+ *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
2634+ node=tip1;
2635+ for (counter=0;counter<numsegs;counter++){
2636+ for (j=0;j<numsegs;j++){
2637+ node1=*(riftsegments+3*j+0);
2638+ node2=*(riftsegments+3*j+1);
2639+
2640+ if ((node1==node) || (node2==node)){
2641+ /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
2642+ already_ordered=0;
2643+ for (k=0;k<counter;k++){
2644+ if(order[k]==j){
2645+ already_ordered=1;
2646+ break;
2647+ }
2648+ }
2649+ if (!already_ordered){
2650+ order[counter]=j;
2651+ if(node1==node){
2652+ node=node2;
2653+ }
2654+ else if(node2==node){
2655+ node=node1;
2656+ }
2657+ break;
2658+ }
2659+ }
2660+ }
2661+ }
2662+
2663+ /*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
2664+ for (j=0;j<numsegs;j++){
2665+ _assert_(order[j]<numsegs);
2666+ *(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
2667+ *(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
2668+ *(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
2669+ *(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
2670+ *(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
2671+ }
2672+
2673+ for (j=0;j<numsegs;j++){
2674+ *(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
2675+ *(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
2676+ *(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
2677+ *(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
2678+ *(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
2679+ }
2680+
2681+ xDelete<int>(order);
2682+ xDelete<int>(riftsegments_copy);
2683+ xDelete<int>(riftpairs_copy);
2684+
2685+ }
2686+
2687+ /*Assign output pointer:*/
2688+ *priftstips=riftstips;
2689+ return noerr;
2690+}/*}}}*/
2691+int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,/*{{{*/
2692+ int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
2693+
2694+ int noerr=1;
2695+ int i,j,k,k0;
2696+
2697+ double el1,el2,node1,node2,node3,node4;
2698+ double temp_node;
2699+
2700+ /*output: */
2701+ double **riftspenaltypairs = NULL;
2702+ double *riftpenaltypairs = NULL;
2703+ int *riftsnumpenaltypairs = NULL;
2704+
2705+ /*intermediary: */
2706+ int numsegs;
2707+ int* riftsegments=NULL;
2708+ int* riftpairs=NULL;
2709+ int counter;
2710+ double normal[2];
2711+ double length;
2712+ int k1,k2;
2713+
2714+ /*Allocate: */
2715+ riftspenaltypairs=xNew<double*>(numrifts);
2716+ riftsnumpenaltypairs=xNew<int>(numrifts);
2717+
2718+ for(i=0;i<numrifts;i++){
2719+ numsegs=riftsnumsegs[i];
2720+ riftsegments=riftssegments[i];
2721+ riftpairs=riftspairs[i];
2722+
2723+ /*allocate riftpenaltypairs, and riftnumpenaltypairs: */
2724+ if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
2725+
2726+ /*Go through only one flank of the rifts, not counting the tips: */
2727+ counter=0;
2728+ for(j=0;j<(numsegs/2);j++){
2729+ el1=*(riftpairs+2*j+0);
2730+ el2=*(riftpairs+2*j+1);
2731+ node1=*(riftsegments+3*j+0);
2732+ node2=*(riftsegments+3*j+1);
2733+ /*Find segment index to recover node3 and node4, facing node1 and node2: */
2734+ k0=-1;
2735+ for(k=0;k<numsegs;k++){
2736+ if(*(riftsegments+3*k+2)==el2){
2737+ k0=k;
2738+ break;
2739+ }
2740+ }
2741+ node3=*(riftsegments+3*k0+0);
2742+ node4=*(riftsegments+3*k0+1);
2743+
2744+ /* Make sure node3 faces node1 and node4 faces node2: */
2745+ if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
2746+ /*Swap node3 and node4:*/
2747+ temp_node=node3;
2748+ node3=node4;
2749+ node4=temp_node;
2750+ }
2751+ /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
2752+ *this segment, and its length: */
2753+ normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
2754+ normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
2755+ length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
2756+
2757+ /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
2758+ * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4,
2759+ * only once. We'll add the normals and the lengths : */
2760+
2761+ if(node1!=node3){ //exclude tips from loads
2762+ k1=-1;
2763+ for(k=0;k<counter;k++){
2764+ if( (*(riftpenaltypairs+k*7+0))==node1){
2765+ k1=k;
2766+ break;
2767+ }
2768+ }
2769+ if(k1==-1){
2770+ *(riftpenaltypairs+counter*7+0)=node1;
2771+ *(riftpenaltypairs+counter*7+1)=node3;
2772+ *(riftpenaltypairs+counter*7+2)=el1;
2773+ *(riftpenaltypairs+counter*7+3)=el2;
2774+ *(riftpenaltypairs+counter*7+4)=normal[0];
2775+ *(riftpenaltypairs+counter*7+5)=normal[1];
2776+ *(riftpenaltypairs+counter*7+6)=length/2;
2777+ counter++;
2778+ }
2779+ else{
2780+ *(riftpenaltypairs+k1*7+4)+=normal[0];
2781+ *(riftpenaltypairs+k1*7+5)+=normal[1];
2782+ *(riftpenaltypairs+k1*7+6)+=length/2;
2783+ }
2784+ }
2785+ if(node2!=node4){
2786+ k2=-1;
2787+ for(k=0;k<counter;k++){
2788+ if( (*(riftpenaltypairs+k*7+0))==node2){
2789+ k2=k;
2790+ break;
2791+ }
2792+ }
2793+ if(k2==-1){
2794+ *(riftpenaltypairs+counter*7+0)=node2;
2795+ *(riftpenaltypairs+counter*7+1)=node4;
2796+ *(riftpenaltypairs+counter*7+2)=el1;
2797+ *(riftpenaltypairs+counter*7+3)=el2;
2798+ *(riftpenaltypairs+counter*7+4)=normal[0];
2799+ *(riftpenaltypairs+counter*7+5)=normal[1];
2800+ *(riftpenaltypairs+counter*7+6)=length/2;
2801+ counter++;
2802+ }
2803+ else{
2804+ *(riftpenaltypairs+k2*7+4)+=normal[0];
2805+ *(riftpenaltypairs+k2*7+5)+=normal[1];
2806+ *(riftpenaltypairs+k2*7+6)+=length/2;
2807+ }
2808+ }
2809+ }
2810+ /*Renormalize normals: */
2811+ for(j=0;j<counter;j++){
2812+ double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
2813+ *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
2814+ *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
2815+ }
2816+
2817+ riftspenaltypairs[i]=riftpenaltypairs;
2818+ riftsnumpenaltypairs[i]=(numsegs/2-1);
2819+ }
2820+
2821+ /*Assign output pointers: */
2822+ *priftspenaltypairs=riftspenaltypairs;
2823+ *priftsnumpenaltypairs=riftsnumpenaltypairs;
2824+ return noerr;
2825+}/*}}}*/
2826+int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
2827+
2828+ int noerr=1;
2829+ int i,j,k;
2830+ int node1,node2,node3;
2831+ int el;
2832+ double pair[2];
2833+ int pair_count=0;
2834+ int triple=0;
2835+
2836+ /*Recover input: */
2837+ int *index = *pindex;
2838+ int nel = *pnel;
2839+ double *x = *px;
2840+ double *y = *py;
2841+ int nods = *pnods;
2842+
2843+ for (i=0;i<num_seg;i++){
2844+ node1=*(segments+3*i+0);
2845+ node2=*(segments+3*i+1);
2846+ /*Find all elements connected to [node1 node2]: */
2847+ pair_count=0;
2848+ for (j=0;j<nel;j++){
2849+ if (*(index+3*j+0)==node1){
2850+ if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
2851+ pair[pair_count]=j;
2852+ pair_count++;
2853+ }
2854+ }
2855+ if (*(index+3*j+1)==node1){
2856+ if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
2857+ pair[pair_count]=j;
2858+ pair_count++;
2859+ }
2860+ }
2861+ if (*(index+3*j+2)==node1){
2862+ if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
2863+ pair[pair_count]=j;
2864+ pair_count++;
2865+ }
2866+ }
2867+ }
2868+ /*Ok, we have pair_count elements connected to this segment. For each of these elements,
2869+ *figure out if the third node also belongs to a segment: */
2870+ if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to 2 elements
2871+ continue;
2872+ }
2873+ else{
2874+ for (j=0;j<pair_count;j++){
2875+ el=(int)pair[j];
2876+ triple=0;
2877+ /*First find node3: */
2878+ if (*(index+3*el+0)==node1){
2879+ if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
2880+ else node3=*(index+3*el+1);
2881+ }
2882+ if (*(index+3*el+1)==node1){
2883+ if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
2884+ else node3=*(index+3*el+0);
2885+ }
2886+ if (*(index+3*el+2)==node1){
2887+ if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
2888+ else node3=*(index+3*el+0);
2889+ }
2890+ /*Ok, we have node3. Does node3 belong to a segment? : */
2891+ for (k=0;k<num_seg;k++){
2892+ if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
2893+ triple=1;
2894+ break;
2895+ }
2896+ }
2897+ if(triple==1){
2898+ /*el is a corner element: we need to split it in 3 triangles: */
2899+ x=xReNew<double>(x,nods,nods+1);
2900+ y=xReNew<double>(y,nods,nods+1);
2901+ x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
2902+ y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
2903+ index=xReNew<int>(index,nel*3,(nel+2*3));
2904+ /*First, reassign element el: */
2905+ *(index+3*el+0)=node1;
2906+ *(index+3*el+1)=node2;
2907+ *(index+3*el+2)=nods+1;
2908+ /*Other two elements: */
2909+ *(index+3*nel+0)=node2;
2910+ *(index+3*nel+1)=node3;
2911+ *(index+3*nel+2)=nods+1;
2912+
2913+ *(index+3*(nel+1)+0)=node3;
2914+ *(index+3*(nel+1)+1)=node1;
2915+ *(index+3*(nel+1)+2)=nods+1;
2916+ /*we need to change the segment elements corresponding to el: */
2917+ for (k=0;k<num_seg;k++){
2918+ if (*(segments+3*k+2)==(el+1)){
2919+ if ( ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node2)) || ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node1))) *(segments+3*k+2)=el+1;
2920+ if ( ((*(segments+3*k+0)==node2) && (*(segments+3*k+1)==node3)) || ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node2))) *(segments+3*k+2)=nel+1;
2921+ if ( ((*(segments+3*k+0)==node3) && (*(segments+3*k+1)==node1)) || ((*(segments+3*k+0)==node1) && (*(segments+3*k+1)==node3))) *(segments+3*k+2)=nel+2;
2922+ }
2923+ }
2924+
2925+ nods=nods+1;
2926+ nel=nel+2;
2927+ i=0;
2928+ break;
2929+ }
2930+ } //for (j=0;j<pair_count;j++)
2931+ }
2932+ }// for (i=0;i<num_seg;i++)
2933+
2934+ /*Assign output pointers: */
2935+ *pindex=index;
2936+ *pnel=nel;
2937+ *px=x;
2938+ *py=y;
2939+ *pnods=nods;
2940+ return noerr;
2941+}/*}}}*/
2942Index: ../trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp
2943===================================================================
2944--- ../trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp (nonexistent)
2945+++ ../trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp (revision 22498)
2946@@ -0,0 +1,96 @@
2947+/*
2948+ * SplitMeshForRifts.c:
2949+ */
2950+#include "./trimesh.h"
2951+#include "../MemOps/MemOps.h"
2952+
2953+int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist){
2954+
2955+ /*Some notes on dimensions:
2956+ index of size nelx3
2957+ x and y of size nodsx1
2958+ segments of size nsegsx3*/
2959+
2960+ int i,j,k,l;
2961+ int node;
2962+ int el;
2963+ int nriftsegs;
2964+ int* riftsegments=NULL;
2965+ int* flags=NULL;
2966+ int NumGridElementListOnOneSideOfRift;
2967+ int* GridElementListOnOneSideOfRift=NULL;
2968+
2969+ /*Recover input: */
2970+ int nel = *pnel;
2971+ int *index = *pindex;
2972+ int nods = *pnods;
2973+ double *x = *px;
2974+ double *y = *py;
2975+ int nsegs = *pnsegs;
2976+ int *segments = *psegments;
2977+ int *segmentmarkerlist = *psegmentmarkerlist;
2978+
2979+ /*Establish list of segments that belong to a rift: */
2980+ /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
2981+ RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);
2982+
2983+ /*Go through all nodes of the rift segments, and start splitting the mesh: */
2984+ flags=xNewZeroInit<int>(nods); //to make sure we don't split the same nodes twice!
2985+ for (i=0;i<nriftsegs;i++){
2986+ for (j=0;j<2;j++){
2987+
2988+ node=riftsegments[4*i+j+2];
2989+ if(flags[node-1]){
2990+ /*This node was already split, skip:*/
2991+ continue;
2992+ }
2993+ else{
2994+ flags[node-1]=1;
2995+ }
2996+
2997+ if(IsGridOnRift(riftsegments,nriftsegs,node)){
2998+
2999+ DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
3000+
3001+ /*Summary: we have for node, a list of elements
3002+ * (GridElementListOnOneSideOfRift, of size
3003+ * NumGridElementListOnOneSideOfRift) that all contain node
3004+ *and that are on the same side of the rift. For all these
3005+ elements, we clone node into another node, and we swap all
3006+ instances of node in the triangulation *for those elements, to the
3007+ new node.*/
3008+
3009+ //create new node
3010+ x=xReNew<double>(x,nods,nods+1);
3011+ y=xReNew<double>(y,nods,nods+1);
3012+ x[nods]=x[node-1]; //matlab indexing
3013+ y[nods]=y[node-1]; //matlab indexing
3014+
3015+ //augment number of nodes
3016+ nods++;
3017+
3018+ //change elements owning this node
3019+ for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
3020+ el=GridElementListOnOneSideOfRift[k];
3021+ for (l=0;l<3;l++){
3022+ if (index[3*el+l]==node) index[3*el+l]=nods; //again, matlab indexing.
3023+ }
3024+ }
3025+ }// if(IsGridOnRift(riftsegments,nriftsegs,node))
3026+ } //for(j=0;j<2;j++)
3027+ } //for (i=0;i<nriftsegs;i++)
3028+
3029+ /*update segments: they got modified completely by adding new nodes.*/
3030+ UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
3031+
3032+ /*Assign output pointers: */
3033+ *pnel=nel;
3034+ *pindex=index;
3035+ *pnods=nods;
3036+ *px=x;
3037+ *py=y;
3038+ *pnsegs=nsegs;
3039+ *psegments=segments;
3040+ *psegmentmarkerlist=segmentmarkerlist;
3041+ return 1;
3042+}
3043Index: ../trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp
3044===================================================================
3045--- ../trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp (nonexistent)
3046+++ ../trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp (revision 22498)
3047@@ -0,0 +1,51 @@
3048+/*
3049+ * GridInsideHole.c:
3050+ * from a convex set of points, figure out a point that for sure lies inside the profile.
3051+ */
3052+
3053+#include <math.h>
3054+
3055+#include "./trimesh.h"
3056+#include "../Exp/exp.h"
3057+
3058+#undef M_PI
3059+#define M_PI 3.141592653589793238462643
3060+
3061+int GridInsideHole(double* px0,double* py0,int n,double* x,double* y){
3062+
3063+ double flag=0.0;
3064+ double xA,xB,xC,xD,xE;
3065+ double yA,yB,yC,yD,yE;
3066+
3067+ /*Take first and last vertices: */
3068+ xA=x[0];
3069+ yA=y[0];
3070+ xB=x[n-1];
3071+ yB=y[n-1];
3072+
3073+ /*Figure out middle of segment [A B]: */
3074+ xC=(xA+xB)/2;
3075+ yC=(yA+yB)/2;
3076+
3077+ /*D and E are on each side of segment [A B], on the median line between segment [A B],
3078+ *at an angle of 10 degree (less than the minimum 30 enforced by the quality of the mesh: */
3079+ xD=xC+tan(10./180.*M_PI)*(yC-yA);
3080+ yD=yC+tan(10./180.*M_PI)*(xA-xC);
3081+ xE=xC-tan(10./180.*M_PI)*(yC-yA);
3082+ yE=yC-tan(10./180.*M_PI)*(xA-xC);
3083+
3084+ /*Either E or D is inside profile (x,y): */
3085+ IsInPolySerial(&flag,&xD,&yD,1,x,y,n,2);
3086+ /*FIXME: used to be 'flag' and not '!flag', check*/
3087+ if(!flag){
3088+ /*D is inside the poly: */
3089+ *px0=xD;
3090+ *py0=yD;
3091+ }
3092+ else{
3093+ /*E is inside the poly: */
3094+ *px0=xE;
3095+ *py0=yE;
3096+ }
3097+ return 1;
3098+}
3099Index: ../trunk-jpl/src/c/shared/Triangle/triangle.h
3100===================================================================
3101--- ../trunk-jpl/src/c/shared/Triangle/triangle.h (nonexistent)
3102+++ ../trunk-jpl/src/c/shared/Triangle/triangle.h (revision 22498)
3103@@ -0,0 +1,33 @@
3104+/*!\file: triangle.h
3105+ * \brief
3106+ */
3107+
3108+#ifndef _SHARED_TRIANGLE_H
3109+#define _SHARED_TRIANGLE_H
3110+
3111+#include <stdio.h>
3112+#include <math.h>
3113+
3114+//#define REAL double //took it out because it may conflict with stdlib.h defines. put back if necessary
3115+int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel);
3116+int OrderSegments(int** psegments,int nseg, int* index,int nel);
3117+int GridInsideHole(double* px0,double* py0,int n,double* x,double* y);
3118+int FindElement(int A,int B,int* index,int nel);
3119+int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist);
3120+int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
3121+int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
3122+int IsNeighbor(int el1,int el2,int* index);
3123+int IsOnRift(int el,int nriftsegs,int* riftsegments);
3124+void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index, int nsegs,int* segments);
3125+int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift,int segmentnumber, int nriftsegs,int* riftsegments, int node,int* index,int nel);
3126+int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
3127+int FindElement(double A,double B,int* index,int nel);
3128+int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
3129+int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels);
3130+int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels);
3131+int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,
3132+ int* riftsnumsegments,int** riftspairs,int* riftstips,double* x,double* y);
3133+int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg);
3134+int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y);
3135+
3136+#endif /* _SHARED_TRIANGLE_H */
3137Index: ../trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp
3138===================================================================
3139--- ../trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp (nonexistent)
3140+++ ../trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp (revision 22498)
3141@@ -0,0 +1,24 @@
3142+/*!\file: AssociateSegmentToElement.cpp
3143+ * \brief for each segment, look for the corresponding element.
3144+ */
3145+
3146+#include "./trimesh.h"
3147+
3148+int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel){
3149+
3150+ /*node indices: */
3151+ int A,B;
3152+
3153+ /*Recover segments: */
3154+ int* segments=*psegments;
3155+
3156+ for(int i=0;i<nseg;i++){
3157+ A=segments[3*i+0];
3158+ B=segments[3*i+1];
3159+ segments[3*i+2]=FindElement(A,B,index,nel)+1; //matlab indexing.
3160+ }
3161+
3162+ /*Assign output pointers: */
3163+ *psegments=segments;
3164+ return 1;
3165+}
3166Index: ../trunk-jpl/src/c/shared/Triangle
3167===================================================================
3168--- ../trunk-jpl/src/c/shared/Triangle (nonexistent)
3169+++ ../trunk-jpl/src/c/shared/Triangle (revision 22498)
3170
3171Property changes on: ../trunk-jpl/src/c/shared/Triangle
3172___________________________________________________________________
3173Added: svn:ignore
3174## -0,0 +1,2 ##
3175+.deps
3176+.dirstamp
3177Index: ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.h
3178===================================================================
3179--- ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.h (revision 22497)
3180+++ ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.h (nonexistent)
3181@@ -1,65 +0,0 @@
3182-/*
3183- * TriMeshProcessRifts.h
3184- */
3185-
3186-#ifndef _TRIMESH_PROCESSRIFTS_H_
3187-#define _TRIMESH_PROCESSRIFTS_H_
3188-
3189-#ifdef HAVE_CONFIG_H
3190-#include <config.h>
3191-#else
3192-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
3193-#endif
3194-
3195-/*For python modules: needs to come before header files inclusion*/
3196-#ifdef _HAVE_PYTHON_
3197-#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
3198-#endif
3199-
3200-#include "../bindings.h"
3201-#include "../../c/main/globals.h"
3202-#include "../../c/modules/modules.h"
3203-#include "../../c/shared/shared.h"
3204-
3205-#undef __FUNCT__
3206-#define __FUNCT__ "TriMeshProcessRifts"
3207-
3208-#ifdef _HAVE_MATLAB_MODULES_
3209-/* serial input macros: */
3210-#define INDEXIN prhs[0]
3211-#define XIN prhs[1]
3212-#define YIN prhs[2]
3213-#define SEGMENTSIN prhs[3]
3214-#define SEGMENTMARKERSIN prhs[4]
3215-/* serial output macros: */
3216-#define INDEXOUT (mxArray**)&plhs[0]
3217-#define XOUT (mxArray**)&plhs[1]
3218-#define YOUT (mxArray**)&plhs[2]
3219-#define SEGMENTSOUT (mxArray**)&plhs[3]
3220-#define SEGMENTMARKERSOUT (mxArray**)&plhs[4]
3221-#define RIFTSTRUCT (mxArray**)&plhs[5]
3222-#endif
3223-
3224-#ifdef _HAVE_PYTHON_MODULES_
3225-/* serial input macros: */
3226-#define INDEXIN PyTuple_GetItem(args,0)
3227-#define XIN PyTuple_GetItem(args,1)
3228-#define YIN PyTuple_GetItem(args,2)
3229-#define SEGMENTSIN PyTuple_GetItem(args,3)
3230-#define SEGMENTMARKERSIN PyTuple_GetItem(args,4)
3231-/* serial output macros: */
3232-#define INDEXOUT output,0
3233-#define XOUT output,1
3234-#define YOUT output,2
3235-#define SEGMENTSOUT output,3
3236-#define SEGMENTMARKERSOUT output,4
3237-#define RIFTSTRUCT output,5
3238-#endif
3239-
3240-/* serial arg counts: */
3241-#undef NLHS
3242-#define NLHS 6
3243-#undef NRHS
3244-#define NRHS 5
3245-
3246-#endif
3247Index: ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp
3248===================================================================
3249--- ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 22497)
3250+++ ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp (nonexistent)
3251@@ -1,59 +0,0 @@
3252-/*!\file: TriMeshProcessRifts.cpp
3253- * \brief split a mesh where a rift (or fault) is present
3254- */
3255-
3256-#include "./TriMeshProcessRifts.h"
3257-
3258-void TriMeshProcessRiftsUsage(void){/*{{{*/
3259- _printf_("\n");
3260- _printf_(" usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(index1,x1,y1,segments1,segmentmarkers1) \n");
3261- _printf_(" where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
3262- _printf_(" index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
3263-}/*}}}*/
3264-WRAPPER(TriMeshProcessRifts_python){
3265-
3266- /* returned quantities: */
3267- RiftStruct *riftstruct = NULL;
3268-
3269- /* input: */
3270- int nel,nods;
3271- int *index = NULL;
3272- double *x = NULL;
3273- double *y = NULL;
3274- int *segments = NULL;
3275- int *segmentmarkers = NULL;
3276- int num_seg;
3277-
3278- /*Boot module*/
3279- MODULEBOOT();
3280-
3281- /*checks on arguments on the matlab side: */
3282- CHECKARGUMENTS(NLHS,NRHS,&TriMeshProcessRiftsUsage);
3283-
3284- /*Fetch data */
3285- FetchData(&index,&nel,NULL,INDEXIN);
3286- FetchData(&x,&nods,XIN);
3287- FetchData(&y,NULL,YIN);
3288- FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
3289- FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
3290-
3291- /*call x layer*/
3292- TriMeshProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,&riftstruct);
3293-
3294- /*Output : */
3295- WriteData(INDEXOUT,index,nel,3);
3296- WriteData(XOUT,x,nods,1);
3297- WriteData(YOUT,y,nods,1);
3298- WriteData(SEGMENTSOUT,segments,num_seg,3);
3299- WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
3300- WriteData(RIFTSTRUCT,riftstruct);
3301-
3302- /*end module: */
3303- delete riftstruct;
3304- xDelete<int>(index);
3305- xDelete<double>(x);
3306- xDelete<double>(y);
3307- xDelete<int>(segments);
3308- xDelete<int>(segmentmarkers );
3309- MODULEEND();
3310-}
3311Index: ../trunk-jpl/src/wrappers/TriMeshProcessRifts
3312===================================================================
3313--- ../trunk-jpl/src/wrappers/TriMeshProcessRifts (revision 22497)
3314+++ ../trunk-jpl/src/wrappers/TriMeshProcessRifts (nonexistent)
3315
3316Property changes on: ../trunk-jpl/src/wrappers/TriMeshProcessRifts
3317___________________________________________________________________
3318Deleted: svn:ignore
3319## -1,2 +0,0 ##
3320-.deps
3321-.dirstamp
3322Index: ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp
3323===================================================================
3324--- ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp (revision 22497)
3325+++ ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp (nonexistent)
3326@@ -1,63 +0,0 @@
3327-/*
3328- * TriMesh: mesh a domain using an .exp file
3329- */
3330-
3331-#include "./TriMesh.h"
3332-
3333-void TriMeshUsage(void){/*{{{*/
3334- _printf_("\n");
3335- _printf_(" usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,area) \n");
3336- _printf_(" where: index,x,y defines a triangulation, segments is an array made \n");
3337- _printf_(" of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment, \n");
3338- _printf_(" outlinefilename an Argus domain outline file, \n");
3339- _printf_(" area is the maximum area desired for any element of the resulting mesh, \n");
3340- _printf_("\n");
3341-}/*}}}*/
3342-WRAPPER(TriMesh_python){
3343-
3344- /*intermediary: */
3345- double area;
3346- Contours *domain = NULL;
3347- Contours *rifts = NULL;
3348-
3349- /* output: */
3350- int *index = NULL;
3351- double *x = NULL;
3352- double *y = NULL;
3353- int *segments = NULL;
3354- int *segmentmarkerlist = NULL;
3355- int nel,nods,nsegs;
3356-
3357- /*Boot module: */
3358- MODULEBOOT();
3359-
3360- /*checks on arguments: */
3361- CHECKARGUMENTS(NLHS,NRHS,&TriMeshUsage);
3362-
3363- /*Fetch data needed for meshing: */
3364- FetchData(&domain,DOMAINOUTLINE);
3365- FetchData(&rifts,RIFTSOUTLINE);
3366- FetchData(&area,AREA);
3367-
3368- /*call x core: */
3369- TriMeshx(&index,&x,&y,&segments,&segmentmarkerlist,&nel,&nods,&nsegs,domain,rifts,area);
3370-
3371- /*write outputs: */
3372- WriteData(INDEX,index,nel,3);
3373- WriteData(X,x,nods);
3374- WriteData(Y,y,nods);
3375- WriteData(SEGMENTS,segments,nsegs,3);
3376- WriteData(SEGMENTMARKERLIST,segmentmarkerlist,nsegs);
3377-
3378- /*free ressources: */
3379- delete domain;
3380- delete rifts;
3381- xDelete<int>(index);
3382- xDelete<double>(x);
3383- xDelete<double>(y);
3384- xDelete<int>(segments);
3385- xDelete<int>(segmentmarkerlist);
3386-
3387- /*end module: */
3388- MODULEEND();
3389-}
3390Index: ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h
3391===================================================================
3392--- ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h (revision 22497)
3393+++ ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h (nonexistent)
3394@@ -1,83 +0,0 @@
3395-/*
3396- TriMesh.h
3397-*/
3398-
3399-#ifndef _TRIMESH_H
3400-#define _TRIMESH_H
3401-
3402-#ifdef HAVE_CONFIG_H
3403- #include <config.h>
3404-#else
3405- #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
3406-#endif
3407-
3408-/*For python modules: needs to come before header files inclusion*/
3409-#ifdef _HAVE_PYTHON_
3410-#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
3411-#endif
3412-
3413-#ifdef _HAVE_JAVASCRIPT_MODULES_
3414-#undef _DO_NOT_LOAD_GLOBALS_ /*only module where this needs to be undefined, so as to
3415- not include IssmComm several times in the javascript Modle construct.*/
3416-#endif
3417-
3418-/*Header files: */
3419-#include "../bindings.h"
3420-#include "../../c/main/globals.h"
3421-#include "../../c/toolkits/toolkits.h"
3422-#include "../../c/modules/modules.h"
3423-#include "../../c/shared/shared.h"
3424-#include "../../c/shared/io/io.h"
3425-
3426-#undef __FUNCT__
3427-#define __FUNCT__ "TriMesh"
3428-
3429-#ifdef _HAVE_MATLAB_MODULES_
3430-/* serial input macros: */
3431-#define DOMAINOUTLINE prhs[0]
3432-#define RIFTSOUTLINE prhs[1]
3433-#define AREA prhs[2]
3434-/* serial output macros: */
3435-#define INDEX (mxArray**)&plhs[0]
3436-#define X (mxArray**)&plhs[1]
3437-#define Y (mxArray**)&plhs[2]
3438-#define SEGMENTS (mxArray**)&plhs[3]
3439-#define SEGMENTMARKERLIST (mxArray**)&plhs[4]
3440-#endif
3441-
3442-#ifdef _HAVE_PYTHON_MODULES_
3443-/* serial input macros: */
3444-#define DOMAINOUTLINE PyTuple_GetItem(args,0)
3445-#define RIFTSOUTLINE PyTuple_GetItem(args,1)
3446-#define AREA PyTuple_GetItem(args,2)
3447-/* serial output macros: */
3448-#define INDEX output,0
3449-#define X output,1
3450-#define Y output,2
3451-#define SEGMENTS output,3
3452-#define SEGMENTMARKERLIST output,4
3453-#endif
3454-
3455-#ifdef _HAVE_JAVASCRIPT_MODULES_
3456-/* serial input macros: */
3457-#define DOMAINOUTLINE domainx,domainy,domainnods
3458-#define RIFTSOUTLINE NULL,NULL,0
3459-#define AREA areain
3460-/* serial output macros: */
3461-#define INDEX pindex,pnel
3462-#define X px,pnods
3463-#define Y py,pnods
3464-#define SEGMENTS psegments,pnsegs
3465-#define SEGMENTMARKERLIST psegmentmarkers,pnsegs
3466-#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)
3467-#define _DO_NOT_LOAD_GLOBALS_//we only load globals for TriMeshModule.js, not other modules!
3468-#endif
3469-
3470-
3471-/* serial arg counts: */
3472-#undef NLHS
3473-#define NLHS 5
3474-#undef NRHS
3475-#define NRHS 3
3476-
3477-#endif /* _TRIMESH_H */
3478Index: ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js
3479===================================================================
3480--- ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js (revision 22497)
3481+++ ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js (nonexistent)
3482@@ -1,82 +0,0 @@
3483-function TriMesh(md,domain,rifts, area){
3484-/*TriMesh
3485- usage: var array = TriMesh(domain,rifts,area);
3486- where: array is made of [index,x,y,segments,segmentmarkers]
3487- and index,x,y defines a triangulation, segments is an array made
3488- of exterior segments to the mesh domain outline, segmentmarkers is an array
3489- flagging each segment, domain a js array defining the domain outline (sames for
3490- rifts) and area is the maximum area desired for any element of the resulting mesh.
3491-
3492- Ok, for now, we are not dealing with rifts. Also, the domain is made of only one
3493- profile, this to avoid passing a double** pointer to js.
3494-*/
3495-
3496- //Dynamic allocations: {{{
3497- //Retrieve domain arrays, and allocate on Module heap:
3498-
3499- //input
3500- var dx=new Float64Array(domain['x']); var nx=dx.length * dx.BYTES_PER_ELEMENT;
3501- var dxPtr= Module._malloc(nx); var domainxHeap = new Uint8Array(Module.HEAPU8.buffer,dxPtr,nx);
3502- domainxHeap.set(new Uint8Array(dx.buffer)); var domainx=domainxHeap.byteOffset;
3503-
3504- var dy=new Float64Array(domain['y']); var ny=dy.length * dy.BYTES_PER_ELEMENT;
3505- var dyPtr = Module._malloc(ny); var domainyHeap = new Uint8Array(Module.HEAPU8.buffer,dyPtr,ny);
3506- domainyHeap.set(new Uint8Array(dy.buffer)); var domainy=domainyHeap.byteOffset;
3507-
3508- //output
3509- var nel,indexlinear,index,nods,x,y;
3510- var pnel= Module._malloc(4);
3511- var pindex= Module._malloc(4);
3512- var pnods= Module._malloc(4);
3513- var px= Module._malloc(4);
3514- var py= Module._malloc(4);
3515- var psegments= Module._malloc(4);
3516- var psegmentmarkers= Module._malloc(4);
3517- var pnsegs= Module._malloc(4);
3518- //}}}
3519-
3520- //Declare TriMesh module:
3521- TriMeshModule = Module.cwrap('TriMeshModule','number',['number','number','number','number','number','number','number','number','number','number','number','number']);
3522-
3523- //Call TriMesh module:
3524- TriMeshModule(pindex,px,py,pnel,pnods,psegments,psegmentmarkers,pnsegs, domainx,domainy,dx.length,area);
3525-
3526- /*Dynamic copying from heap: {{{*/
3527- //recover mesh:
3528- nel = Module.getValue(pnel, 'i32');
3529- var indexptr = Module.getValue(pindex,'i32');
3530- indexlinear = Module.HEAPF64.slice(indexptr /8, indexptr/8 + nel*3);
3531- index = ListToMatrix(indexlinear,3);
3532-
3533- nods = Module.getValue(pnods, 'i32');
3534- var xptr = Module.getValue(px,'i32');
3535- var yptr = Module.getValue(py,'i32');
3536- x = Module.HEAPF64.slice(xptr /8, xptr/8 + nods);
3537- y = Module.HEAPF64.slice(yptr /8, yptr/8 + nods);
3538-
3539- nsegs = Module.getValue(pnsegs, 'i32');
3540- var segmentsptr = Module.getValue(psegments,'i32');
3541- segmentslinear = Module.HEAPF64.slice(segmentsptr /8, segmentsptr/8 + nsegs*3);
3542- segments = ListToMatrix(segmentslinear,3);
3543-
3544- var segmentmarkersptr = Module.getValue(psegmentmarkers,'i32');
3545- segmentmarkers = Module.HEAPF64.slice(segmentmarkersptr /8, segmentmarkersptr/8 + nsegs);
3546- /*}}}*/
3547-
3548- var return_array=[index,x,y,segments,segmentmarkers];
3549-
3550- /*Free ressources: */
3551- Module._free(pindex);
3552- Module._free(indexlinear);
3553- Module._free(px);
3554- Module._free(x);
3555- Module._free(py);
3556- Module._free(y);
3557- Module._free(pnel);
3558- Module._free(pnods);
3559- Module._free(psegments);
3560- Module._free(psegmentmarkers);
3561- Module._free(pnsegs);
3562-
3563- return return_array;
3564-}
3565Index: ../trunk-jpl/src/wrappers/TriMesh
3566===================================================================
3567--- ../trunk-jpl/src/wrappers/TriMesh (revision 22497)
3568+++ ../trunk-jpl/src/wrappers/TriMesh (nonexistent)
3569
3570Property changes on: ../trunk-jpl/src/wrappers/TriMesh
3571___________________________________________________________________
3572Deleted: svn:ignore
3573## -1,2 +0,0 ##
3574-.deps
3575-.dirstamp
3576Index: ../trunk-jpl/src/wrappers/Triangle/Triangle.js
3577===================================================================
3578--- ../trunk-jpl/src/wrappers/Triangle/Triangle.js (nonexistent)
3579+++ ../trunk-jpl/src/wrappers/Triangle/Triangle.js (revision 22498)
3580@@ -0,0 +1,82 @@
3581+function Triangle(md,domain,rifts, area){
3582+/*Triangle
3583+ usage: var array = Triangle(domain,rifts,area);
3584+ where: array is made of [index,x,y,segments,segmentmarkers]
3585+ and index,x,y defines a triangulation, segments is an array made
3586+ of exterior segments to the mesh domain outline, segmentmarkers is an array
3587+ flagging each segment, domain a js array defining the domain outline (sames for
3588+ rifts) and area is the maximum area desired for any element of the resulting mesh.
3589+
3590+ Ok, for now, we are not dealing with rifts. Also, the domain is made of only one
3591+ profile, this to avoid passing a double** pointer to js.
3592+*/
3593+
3594+ //Dynamic allocations: {{{
3595+ //Retrieve domain arrays, and allocate on Module heap:
3596+
3597+ //input
3598+ var dx=new Float64Array(domain['x']); var nx=dx.length * dx.BYTES_PER_ELEMENT;
3599+ var dxPtr= Module._malloc(nx); var domainxHeap = new Uint8Array(Module.HEAPU8.buffer,dxPtr,nx);
3600+ domainxHeap.set(new Uint8Array(dx.buffer)); var domainx=domainxHeap.byteOffset;
3601+
3602+ var dy=new Float64Array(domain['y']); var ny=dy.length * dy.BYTES_PER_ELEMENT;
3603+ var dyPtr = Module._malloc(ny); var domainyHeap = new Uint8Array(Module.HEAPU8.buffer,dyPtr,ny);
3604+ domainyHeap.set(new Uint8Array(dy.buffer)); var domainy=domainyHeap.byteOffset;
3605+
3606+ //output
3607+ var nel,indexlinear,index,nods,x,y;
3608+ var pnel= Module._malloc(4);
3609+ var pindex= Module._malloc(4);
3610+ var pnods= Module._malloc(4);
3611+ var px= Module._malloc(4);
3612+ var py= Module._malloc(4);
3613+ var psegments= Module._malloc(4);
3614+ var psegmentmarkers= Module._malloc(4);
3615+ var pnsegs= Module._malloc(4);
3616+ //}}}
3617+
3618+ //Declare Triangle module:
3619+ TriangleModule = Module.cwrap('TriangleModule','number',['number','number','number','number','number','number','number','number','number','number','number','number']);
3620+
3621+ //Call Triangle module:
3622+ TriangleModule(pindex,px,py,pnel,pnods,psegments,psegmentmarkers,pnsegs, domainx,domainy,dx.length,area);
3623+
3624+ /*Dynamic copying from heap: {{{*/
3625+ //recover mesh:
3626+ nel = Module.getValue(pnel, 'i32');
3627+ var indexptr = Module.getValue(pindex,'i32');
3628+ indexlinear = Module.HEAPF64.slice(indexptr /8, indexptr/8 + nel*3);
3629+ index = ListToMatrix(indexlinear,3);
3630+
3631+ nods = Module.getValue(pnods, 'i32');
3632+ var xptr = Module.getValue(px,'i32');
3633+ var yptr = Module.getValue(py,'i32');
3634+ x = Module.HEAPF64.slice(xptr /8, xptr/8 + nods);
3635+ y = Module.HEAPF64.slice(yptr /8, yptr/8 + nods);
3636+
3637+ nsegs = Module.getValue(pnsegs, 'i32');
3638+ var segmentsptr = Module.getValue(psegments,'i32');
3639+ segmentslinear = Module.HEAPF64.slice(segmentsptr /8, segmentsptr/8 + nsegs*3);
3640+ segments = ListToMatrix(segmentslinear,3);
3641+
3642+ var segmentmarkersptr = Module.getValue(psegmentmarkers,'i32');
3643+ segmentmarkers = Module.HEAPF64.slice(segmentmarkersptr /8, segmentmarkersptr/8 + nsegs);
3644+ /*}}}*/
3645+
3646+ var return_array=[index,x,y,segments,segmentmarkers];
3647+
3648+ /*Free ressources: */
3649+ Module._free(pindex);
3650+ Module._free(indexlinear);
3651+ Module._free(px);
3652+ Module._free(x);
3653+ Module._free(py);
3654+ Module._free(y);
3655+ Module._free(pnel);
3656+ Module._free(pnods);
3657+ Module._free(psegments);
3658+ Module._free(psegmentmarkers);
3659+ Module._free(pnsegs);
3660+
3661+ return return_array;
3662+}
3663Index: ../trunk-jpl/src/wrappers/Triangle/Triangle.cpp
3664===================================================================
3665--- ../trunk-jpl/src/wrappers/Triangle/Triangle.cpp (nonexistent)
3666+++ ../trunk-jpl/src/wrappers/Triangle/Triangle.cpp (revision 22498)
3667@@ -0,0 +1,63 @@
3668+/*
3669+ * Triangle: mesh a domain using an .exp file
3670+ */
3671+
3672+#include "./Triangle.h"
3673+
3674+void TriangleUsage(void){/*{{{*/
3675+ _printf_("\n");
3676+ _printf_(" usage: [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,area) \n");
3677+ _printf_(" where: index,x,y defines a triangulation, segments is an array made \n");
3678+ _printf_(" of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment, \n");
3679+ _printf_(" outlinefilename an Argus domain outline file, \n");
3680+ _printf_(" area is the maximum area desired for any element of the resulting mesh, \n");
3681+ _printf_("\n");
3682+}/*}}}*/
3683+WRAPPER(Triangle_python){
3684+
3685+ /*intermediary: */
3686+ double area;
3687+ Contours *domain = NULL;
3688+ Contours *rifts = NULL;
3689+
3690+ /* output: */
3691+ int *index = NULL;
3692+ double *x = NULL;
3693+ double *y = NULL;
3694+ int *segments = NULL;
3695+ int *segmentmarkerlist = NULL;
3696+ int nel,nods,nsegs;
3697+
3698+ /*Boot module: */
3699+ MODULEBOOT();
3700+
3701+ /*checks on arguments: */
3702+ CHECKARGUMENTS(NLHS,NRHS,&TriangleUsage);
3703+
3704+ /*Fetch data needed for meshing: */
3705+ FetchData(&domain,DOMAINOUTLINE);
3706+ FetchData(&rifts,RIFTSOUTLINE);
3707+ FetchData(&area,AREA);
3708+
3709+ /*call x core: */
3710+ Trianglex(&index,&x,&y,&segments,&segmentmarkerlist,&nel,&nods,&nsegs,domain,rifts,area);
3711+
3712+ /*write outputs: */
3713+ WriteData(INDEX,index,nel,3);
3714+ WriteData(X,x,nods);
3715+ WriteData(Y,y,nods);
3716+ WriteData(SEGMENTS,segments,nsegs,3);
3717+ WriteData(SEGMENTMARKERLIST,segmentmarkerlist,nsegs);
3718+
3719+ /*free ressources: */
3720+ delete domain;
3721+ delete rifts;
3722+ xDelete<int>(index);
3723+ xDelete<double>(x);
3724+ xDelete<double>(y);
3725+ xDelete<int>(segments);
3726+ xDelete<int>(segmentmarkerlist);
3727+
3728+ /*end module: */
3729+ MODULEEND();
3730+}
3731Index: ../trunk-jpl/src/wrappers/Triangle/Triangle.h
3732===================================================================
3733--- ../trunk-jpl/src/wrappers/Triangle/Triangle.h (nonexistent)
3734+++ ../trunk-jpl/src/wrappers/Triangle/Triangle.h (revision 22498)
3735@@ -0,0 +1,83 @@
3736+/*
3737+ Triangle.h
3738+*/
3739+
3740+#ifndef _TRIANGLE_H
3741+#define _TRIANGLE_H
3742+
3743+#ifdef HAVE_CONFIG_H
3744+ #include <config.h>
3745+#else
3746+ #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
3747+#endif
3748+
3749+/*For python modules: needs to come before header files inclusion*/
3750+#ifdef _HAVE_PYTHON_
3751+#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
3752+#endif
3753+
3754+#ifdef _HAVE_JAVASCRIPT_MODULES_
3755+#undef _DO_NOT_LOAD_GLOBALS_ /*only module where this needs to be undefined, so as to
3756+ not include IssmComm several times in the javascript Modle construct.*/
3757+#endif
3758+
3759+/*Header files: */
3760+#include "../bindings.h"
3761+#include "../../c/main/globals.h"
3762+#include "../../c/toolkits/toolkits.h"
3763+#include "../../c/modules/modules.h"
3764+#include "../../c/shared/shared.h"
3765+#include "../../c/shared/io/io.h"
3766+
3767+#undef __FUNCT__
3768+#define __FUNCT__ "Triangle"
3769+
3770+#ifdef _HAVE_MATLAB_MODULES_
3771+/* serial input macros: */
3772+#define DOMAINOUTLINE prhs[0]
3773+#define RIFTSOUTLINE prhs[1]
3774+#define AREA prhs[2]
3775+/* serial output macros: */
3776+#define INDEX (mxArray**)&plhs[0]
3777+#define X (mxArray**)&plhs[1]
3778+#define Y (mxArray**)&plhs[2]
3779+#define SEGMENTS (mxArray**)&plhs[3]
3780+#define SEGMENTMARKERLIST (mxArray**)&plhs[4]
3781+#endif
3782+
3783+#ifdef _HAVE_PYTHON_MODULES_
3784+/* serial input macros: */
3785+#define DOMAINOUTLINE PyTuple_GetItem(args,0)
3786+#define RIFTSOUTLINE PyTuple_GetItem(args,1)
3787+#define AREA PyTuple_GetItem(args,2)
3788+/* serial output macros: */
3789+#define INDEX output,0
3790+#define X output,1
3791+#define Y output,2
3792+#define SEGMENTS output,3
3793+#define SEGMENTMARKERLIST output,4
3794+#endif
3795+
3796+#ifdef _HAVE_JAVASCRIPT_MODULES_
3797+/* serial input macros: */
3798+#define DOMAINOUTLINE domainx,domainy,domainnods
3799+#define RIFTSOUTLINE NULL,NULL,0
3800+#define AREA areain
3801+/* serial output macros: */
3802+#define INDEX pindex,pnel
3803+#define X px,pnods
3804+#define Y py,pnods
3805+#define SEGMENTS psegments,pnsegs
3806+#define SEGMENTMARKERLIST psegmentmarkers,pnsegs
3807+#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)
3808+#define _DO_NOT_LOAD_GLOBALS_//we only load globals for TriangleModule.js, not other modules!
3809+#endif
3810+
3811+
3812+/* serial arg counts: */
3813+#undef NLHS
3814+#define NLHS 5
3815+#undef NRHS
3816+#define NRHS 3
3817+
3818+#endif /* _TRIANGLE_H */
3819Index: ../trunk-jpl/src/wrappers/Triangle
3820===================================================================
3821--- ../trunk-jpl/src/wrappers/Triangle (nonexistent)
3822+++ ../trunk-jpl/src/wrappers/Triangle (revision 22498)
3823
3824Property changes on: ../trunk-jpl/src/wrappers/Triangle
3825___________________________________________________________________
3826Added: svn:ignore
3827## -0,0 +1,2 ##
3828+.deps
3829+.dirstamp
3830Index: ../trunk-jpl/src/wrappers/javascript/Makefile.am
3831===================================================================
3832--- ../trunk-jpl/src/wrappers/javascript/Makefile.am (revision 22497)
3833+++ ../trunk-jpl/src/wrappers/javascript/Makefile.am (revision 22498)
3834@@ -6,7 +6,7 @@
3835 #define prefix (from http://www.gnu.org/software/autoconf/manual/autoconf-2.67/html_node/Defining-Directories.html)
3836 AM_CPPFLAGS+= -DISSM_PREFIX='"$(prefix)"'
3837
3838-js_scripts = ${ISSM_DIR}/src/wrappers/TriMesh/TriMesh.js \
3839+js_scripts = ${ISSM_DIR}/src/wrappers/Triangle/Triangle.js \
3840 ${ISSM_DIR}/src/wrappers/NodeConnectivity/NodeConnectivity.js\
3841 ${ISSM_DIR}/src/wrappers/ContourToMesh/ContourToMesh.js\
3842 ${ISSM_DIR}/src/wrappers/ElementConnectivity/ElementConnectivity.js\
3843@@ -80,7 +80,7 @@
3844 libISSMApi_la_LDFLAGS = -static
3845 endif
3846
3847-IssmModule_SOURCES = ../TriMesh/TriMesh.cpp \
3848+IssmModule_SOURCES = ../Triangle/Triangle.cpp \
3849 ../NodeConnectivity/NodeConnectivity.cpp\
3850 ../ContourToMesh/ContourToMesh.cpp\
3851 ../ElementConnectivity/ElementConnectivity.cpp\
3852@@ -88,6 +88,6 @@
3853 ../IssmConfig/IssmConfig.cpp\
3854 ../Issm/issm.cpp
3855
3856-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
3857+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
3858 IssmModule_LDADD = ${deps} $(TRIANGLELIB) $(GSLLIB)
3859 #}}}
3860Index: ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp
3861===================================================================
3862--- ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp (nonexistent)
3863+++ ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp (revision 22498)
3864@@ -0,0 +1,59 @@
3865+/*!\file: ProcessRifts.cpp
3866+ * \brief split a mesh where a rift (or fault) is present
3867+ */
3868+
3869+#include "./ProcessRifts.h"
3870+
3871+void ProcessRiftsUsage(void){/*{{{*/
3872+ _printf_("\n");
3873+ _printf_(" usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1) \n");
3874+ _printf_(" where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
3875+ _printf_(" index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
3876+}/*}}}*/
3877+WRAPPER(ProcessRifts_python){
3878+
3879+ /* returned quantities: */
3880+ RiftStruct *riftstruct = NULL;
3881+
3882+ /* input: */
3883+ int nel,nods;
3884+ int *index = NULL;
3885+ double *x = NULL;
3886+ double *y = NULL;
3887+ int *segments = NULL;
3888+ int *segmentmarkers = NULL;
3889+ int num_seg;
3890+
3891+ /*Boot module*/
3892+ MODULEBOOT();
3893+
3894+ /*checks on arguments on the matlab side: */
3895+ CHECKARGUMENTS(NLHS,NRHS,&ProcessRiftsUsage);
3896+
3897+ /*Fetch data */
3898+ FetchData(&index,&nel,NULL,INDEXIN);
3899+ FetchData(&x,&nods,XIN);
3900+ FetchData(&y,NULL,YIN);
3901+ FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
3902+ FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
3903+
3904+ /*call x layer*/
3905+ ProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,&riftstruct);
3906+
3907+ /*Output : */
3908+ WriteData(INDEXOUT,index,nel,3);
3909+ WriteData(XOUT,x,nods,1);
3910+ WriteData(YOUT,y,nods,1);
3911+ WriteData(SEGMENTSOUT,segments,num_seg,3);
3912+ WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
3913+ WriteData(RIFTSTRUCT,riftstruct);
3914+
3915+ /*end module: */
3916+ delete riftstruct;
3917+ xDelete<int>(index);
3918+ xDelete<double>(x);
3919+ xDelete<double>(y);
3920+ xDelete<int>(segments);
3921+ xDelete<int>(segmentmarkers );
3922+ MODULEEND();
3923+}
3924Index: ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h
3925===================================================================
3926--- ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h (nonexistent)
3927+++ ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h (revision 22498)
3928@@ -0,0 +1,65 @@
3929+/*
3930+ * ProcessRifts.h
3931+ */
3932+
3933+#ifndef _PROCESSRIFTS_H_
3934+#define _PROCESSRIFTS_H_
3935+
3936+#ifdef HAVE_CONFIG_H
3937+#include <config.h>
3938+#else
3939+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
3940+#endif
3941+
3942+/*For python modules: needs to come before header files inclusion*/
3943+#ifdef _HAVE_PYTHON_
3944+#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
3945+#endif
3946+
3947+#include "../bindings.h"
3948+#include "../../c/main/globals.h"
3949+#include "../../c/modules/modules.h"
3950+#include "../../c/shared/shared.h"
3951+
3952+#undef __FUNCT__
3953+#define __FUNCT__ "ProcessRifts"
3954+
3955+#ifdef _HAVE_MATLAB_MODULES_
3956+/* serial input macros: */
3957+#define INDEXIN prhs[0]
3958+#define XIN prhs[1]
3959+#define YIN prhs[2]
3960+#define SEGMENTSIN prhs[3]
3961+#define SEGMENTMARKERSIN prhs[4]
3962+/* serial output macros: */
3963+#define INDEXOUT (mxArray**)&plhs[0]
3964+#define XOUT (mxArray**)&plhs[1]
3965+#define YOUT (mxArray**)&plhs[2]
3966+#define SEGMENTSOUT (mxArray**)&plhs[3]
3967+#define SEGMENTMARKERSOUT (mxArray**)&plhs[4]
3968+#define RIFTSTRUCT (mxArray**)&plhs[5]
3969+#endif
3970+
3971+#ifdef _HAVE_PYTHON_MODULES_
3972+/* serial input macros: */
3973+#define INDEXIN PyTuple_GetItem(args,0)
3974+#define XIN PyTuple_GetItem(args,1)
3975+#define YIN PyTuple_GetItem(args,2)
3976+#define SEGMENTSIN PyTuple_GetItem(args,3)
3977+#define SEGMENTMARKERSIN PyTuple_GetItem(args,4)
3978+/* serial output macros: */
3979+#define INDEXOUT output,0
3980+#define XOUT output,1
3981+#define YOUT output,2
3982+#define SEGMENTSOUT output,3
3983+#define SEGMENTMARKERSOUT output,4
3984+#define RIFTSTRUCT output,5
3985+#endif
3986+
3987+/* serial arg counts: */
3988+#undef NLHS
3989+#define NLHS 6
3990+#undef NRHS
3991+#define NRHS 5
3992+
3993+#endif
3994Index: ../trunk-jpl/src/wrappers/ProcessRifts
3995===================================================================
3996--- ../trunk-jpl/src/wrappers/ProcessRifts (nonexistent)
3997+++ ../trunk-jpl/src/wrappers/ProcessRifts (revision 22498)
3998
3999Property changes on: ../trunk-jpl/src/wrappers/ProcessRifts
4000___________________________________________________________________
4001Added: svn:ignore
4002## -0,0 +1,2 ##
4003+.deps
4004+.dirstamp
4005Index: ../trunk-jpl/src/wrappers/matlab/Makefile.am
4006===================================================================
4007--- ../trunk-jpl/src/wrappers/matlab/Makefile.am (revision 22497)
4008+++ ../trunk-jpl/src/wrappers/matlab/Makefile.am (revision 22498)
4009@@ -57,8 +57,8 @@
4010 MeshProfileIntersection_matlab.la\
4011 PointCloudFindNeighbors_matlab.la\
4012 PropagateFlagsFromConnectivity_matlab.la\
4013- TriMesh_matlab.la\
4014- TriMeshProcessRifts_matlab.la\
4015+ Triangle_matlab.la\
4016+ ProcessRifts_matlab.la\
4017 Scotch_matlab.la
4018
4019 if CHACO
4020@@ -232,11 +232,11 @@
4021 ShpRead_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
4022 ShpRead_matlab_la_LIBADD = ${deps} $(SHAPELIBLIB) $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
4023
4024-TriMesh_matlab_la_SOURCES = ../TriMesh/TriMesh.cpp
4025-TriMesh_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
4026-TriMesh_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(TRIANGLELIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
4027+Triangle_matlab_la_SOURCES = ../Triangle/Triangle.cpp
4028+Triangle_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
4029+Triangle_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(TRIANGLELIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
4030
4031-TriMeshProcessRifts_matlab_la_SOURCES = ../TriMeshProcessRifts/TriMeshProcessRifts.cpp
4032-TriMeshProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
4033-TriMeshProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
4034+ProcessRifts_matlab_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp
4035+ProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
4036+ProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
4037 #}}}
4038Index: ../trunk-jpl/src/wrappers/python/Makefile.am
4039===================================================================
4040--- ../trunk-jpl/src/wrappers/python/Makefile.am (revision 22497)
4041+++ ../trunk-jpl/src/wrappers/python/Makefile.am (revision 22498)
4042@@ -40,8 +40,8 @@
4043 IssmConfig_python.la\
4044 MeshProfileIntersection_python.la\
4045 NodeConnectivity_python.la\
4046- TriMesh_python.la\
4047- TriMeshProcessRifts_python.la
4048+ Triangle_python.la\
4049+ ProcessRifts_python.la
4050 endif
4051 #}}}
4052 #Flags and libraries {{{
4053@@ -140,11 +140,11 @@
4054 NodeConnectivity_python_la_CXXFLAGS = ${AM_CXXFLAGS}
4055 NodeConnectivity_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
4056
4057-TriMesh_python_la_SOURCES = ../TriMesh/TriMesh.cpp
4058-TriMesh_python_la_CXXFLAGS = ${AM_CXXFLAGS}
4059-TriMesh_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)
4060+Triangle_python_la_SOURCES = ../Triangle/Triangle.cpp
4061+Triangle_python_la_CXXFLAGS = ${AM_CXXFLAGS}
4062+Triangle_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)
4063
4064-TriMeshProcessRifts_python_la_SOURCES = ../TriMeshProcessRifts/TriMeshProcessRifts.cpp
4065-TriMeshProcessRifts_python_la_CXXFLAGS = ${AM_CXXFLAGS}
4066-TriMeshProcessRifts_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
4067+ProcessRifts_python_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp
4068+ProcessRifts_python_la_CXXFLAGS = ${AM_CXXFLAGS}
4069+ProcessRifts_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
4070 #}}}
4071Index: ../trunk-jpl/src/m/mesh/triangle.js
4072===================================================================
4073--- ../trunk-jpl/src/m/mesh/triangle.js (revision 22497)
4074+++ ../trunk-jpl/src/m/mesh/triangle.js (revision 22498)
4075@@ -1,7 +1,7 @@
4076 function triangle(md){
4077 //TRIANGLE - create model mesh using the triangle package
4078 //
4079-// This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
4080+// This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
4081 // where md is a @model object, domainname is the name of an Argus domain outline file,
4082 // and resolution is a characteristic length for the mesh (same unit as the domain outline
4083 // unit). Riftname is an optional argument (Argus domain outline) describing rifts.
4084@@ -35,7 +35,7 @@
4085 var area=Math.pow(resolution,2);
4086
4087 //Call mesher:
4088- var return_array=TriMesh(md, domain, rifts, area);
4089+ var return_array=Triangle(md, domain, rifts, area);
4090
4091 //Plug into md:
4092 md.mesh.elements=return_array[0];
4093Index: ../trunk-jpl/src/m/mesh/triangle2dvertical.m
4094===================================================================
4095--- ../trunk-jpl/src/m/mesh/triangle2dvertical.m (revision 22497)
4096+++ ../trunk-jpl/src/m/mesh/triangle2dvertical.m (revision 22498)
4097@@ -1,7 +1,7 @@
4098 function md=triangle(md,domainname,resolution)
4099 %TRIANGLE - create model mesh using the triangle package
4100 %
4101-% This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
4102+% This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
4103 % where md is a @model object, domainname is the name of an Argus domain outline file,
4104 % and resolution is a characteristic length for the mesh (same unit as the domain outline
4105 % unit). Riftname is an optional argument (Argus domain outline) describing rifts.
4106@@ -23,8 +23,8 @@
4107
4108 area=resolution^2;
4109
4110-%Mesh using TriMesh
4111-[elements,x,z,segments,segmentmarkers]=TriMesh(domainname,'',area);
4112+%Mesh using Triangle
4113+[elements,x,z,segments,segmentmarkers]=Triangle(domainname,'',area);
4114
4115 %check that all the created nodes belong to at least one element
4116 orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
4117Index: ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py
4118===================================================================
4119--- ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 22497)
4120+++ ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 22498)
4121@@ -1,5 +1,5 @@
4122 import numpy as np
4123-from TriMeshProcessRifts import TriMeshProcessRifts
4124+from ProcessRifts import ProcessRifts
4125 from ContourToMesh import ContourToMesh
4126 from meshprocessoutsiderifts import meshprocessoutsiderifts
4127 from GetAreas import GetAreas
4128@@ -21,7 +21,7 @@
4129 """
4130
4131 #Call MEX file
4132- 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)
4133+ 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)
4134 md.mesh.elements=md.mesh.elements.astype(int)
4135 md.mesh.x=md.mesh.x.reshape(-1)
4136 md.mesh.y=md.mesh.y.reshape(-1)
4137@@ -28,7 +28,7 @@
4138 md.mesh.segments=md.mesh.segments.astype(int)
4139 md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
4140 if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
4141- raise RuntimeError("TriMeshProcessRifts did not find any rift")
4142+ raise RuntimeError("ProcessRifts did not find any rift")
4143
4144 #Fill in rest of fields:
4145 numrifts=len(md.rifts.riftstruct)
4146Index: ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m
4147===================================================================
4148--- ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m (revision 22497)
4149+++ ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m (revision 22498)
4150@@ -24,9 +24,9 @@
4151 end
4152
4153 %Call MEX file
4154-[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);
4155+[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);
4156 if ~isstruct(md.rifts.riftstruct),
4157- error('TriMeshProcessRifts did not find any rift');
4158+ error('ProcessRifts did not find any rift');
4159 end
4160
4161 %Fill in rest of fields:
4162Index: ../trunk-jpl/src/m/mesh/argusmesh.m
4163===================================================================
4164--- ../trunk-jpl/src/m/mesh/argusmesh.m (revision 22497)
4165+++ ../trunk-jpl/src/m/mesh/argusmesh.m (revision 22498)
4166@@ -8,7 +8,7 @@
4167 % md=argusmesh(md,infile)
4168 %
4169 % Example:
4170-% md=argusmesh(md,'TriMesh.exp')
4171+% md=argusmesh(md,'Domain.exp')
4172
4173 %some argument check:
4174 if nargin~=2 | nargout~=1,
4175Index: ../trunk-jpl/src/m/mesh/triangle.py
4176===================================================================
4177--- ../trunk-jpl/src/m/mesh/triangle.py (revision 22497)
4178+++ ../trunk-jpl/src/m/mesh/triangle.py (revision 22498)
4179@@ -1,7 +1,7 @@
4180 import os.path
4181 import numpy as np
4182 from mesh2d import mesh2d
4183-from TriMesh import TriMesh
4184+from Triangle import Triangle
4185 from NodeConnectivity import NodeConnectivity
4186 from ElementConnectivity import ElementConnectivity
4187 import MatlabFuncs as m
4188@@ -10,7 +10,7 @@
4189 """
4190 TRIANGLE - create model mesh using the triangle package
4191
4192- This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
4193+ This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
4194 where md is a @model object, domainname is the name of an Argus domain outline file,
4195 and resolution is a characteristic length for the mesh (same unit as the domain outline
4196 unit). Riftname is an optional argument (Argus domain outline) describing rifts.
4197@@ -47,9 +47,9 @@
4198 if not os.path.exists(domainname):
4199 raise IOError("file '%s' not found" % domainname)
4200
4201- #Mesh using TriMesh
4202+ #Mesh using Triangle
4203 md.mesh=mesh2d()
4204- md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=TriMesh(domainname,riftname,area)
4205+ md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=Triangle(domainname,riftname,area)
4206 md.mesh.elements=md.mesh.elements.astype(int)
4207 md.mesh.segments=md.mesh.segments.astype(int)
4208 md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
4209Index: ../trunk-jpl/src/m/mesh/triangle.m
4210===================================================================
4211--- ../trunk-jpl/src/m/mesh/triangle.m (revision 22497)
4212+++ ../trunk-jpl/src/m/mesh/triangle.m (revision 22498)
4213@@ -1,7 +1,7 @@
4214 function md=triangle(md,domainname,varargin)
4215 %TRIANGLE - create model mesh using the triangle package
4216 %
4217-% This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
4218+% This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
4219 % where md is a @model object, domainname is the name of an Argus domain outline file,
4220 % and resolution is a characteristic length for the mesh (same unit as the domain outline
4221 % unit). Riftname is an optional argument (Argus domain outline) describing rifts.
4222@@ -41,8 +41,8 @@
4223 error(['file "' domainname '" not found']);
4224 end
4225
4226-%Mesh using TriMesh
4227-[elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area);
4228+%Mesh using Triangle
4229+[elements,x,y,segments,segmentmarkers]=Triangle(domainname,riftname,area);
4230
4231 %check that all the created nodes belong to at least one element
4232 removeorphans=1;
4233Index: ../trunk-jpl/src/m/modules/TriMesh.m
4234===================================================================
4235--- ../trunk-jpl/src/m/modules/TriMesh.m (revision 22497)
4236+++ ../trunk-jpl/src/m/modules/TriMesh.m (nonexistent)
4237@@ -1,21 +0,0 @@
4238-function [index,x,y,segments,segmentmarkers] = TriMesh(domainoutlinefilename,rifts,mesh_area);
4239-%TRIMESH - Mesh a domain using an .exp file
4240-%
4241-% Usage:
4242-% [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,mesh_area);
4243-%
4244-% index,x,y: Defines a triangulation
4245-% segments: Array made of exterior segments to the mesh domain outline
4246-% segmentmarkers: Array flagging each segment
4247-%
4248-% domainoutlinefilename: Argus domain outline file
4249-% mesh_area: Maximum area desired for any element of the resulting mesh
4250-
4251-% Check usage
4252-if nargin~=3 && nargout~=5
4253- help TriMesh
4254- error('Wrong usage (see above)');
4255-end
4256-
4257-% Call mex module
4258-[index,x,y,segments,segmentmarkers]=TriMesh_matlab(domainoutlinefilename,rifts,mesh_area);
4259Index: ../trunk-jpl/src/m/modules/TriMesh.py
4260===================================================================
4261--- ../trunk-jpl/src/m/modules/TriMesh.py (revision 22497)
4262+++ ../trunk-jpl/src/m/modules/TriMesh.py (nonexistent)
4263@@ -1,21 +0,0 @@
4264-from TriMesh_python import TriMesh_python
4265-
4266-def TriMesh(domainoutlinefilename,rifts,mesh_area):
4267- """
4268- TRIMESH - Mesh a domain using an .exp file
4269-
4270- Usage:
4271- [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,mesh_area);
4272-
4273- index,x,y: defines a triangulation
4274- segments: An array made of exterior segments to the mesh domain outline
4275- segmentmarkers: An array flagging each segment
4276-
4277- domainoutlinefilename: an Argus domain outline file
4278- mesh_area: The maximum area desired for any element of the resulting mesh
4279- """
4280- # Call mex module
4281- index,x,y,segments,segmentmarkers=TriMesh_python(domainoutlinefilename,rifts,mesh_area)
4282- # Return
4283- return index,x,y,segments,segmentmarkers
4284-
4285Index: ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py
4286===================================================================
4287--- ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py (revision 22497)
4288+++ ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py (nonexistent)
4289@@ -1,16 +0,0 @@
4290-from TriMeshProcessRifts_python import TriMeshProcessRifts_python
4291-
4292-def TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1):
4293- """
4294- TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
4295-
4296- Usage:
4297- [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
4298-
4299- (index1,x1,y1,segments1,segmentmarkers1): An initial triangulation.
4300- [index2,x2,y2,segments2,segmentmarkers2,rifts2]: The resulting triangulation where rifts have been processed.
4301- """
4302- # Call mex module
4303- index2,x2,y2,segments2,segmentmarkers2,rifts2 = TriMeshProcessRifts_python(index1,x1,y1,segments1,segmentmarkers1)
4304- # Return
4305- return index2,x2,y2,segments2,segmentmarkers2,rifts2
4306Index: ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m
4307===================================================================
4308--- ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m (revision 22497)
4309+++ ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m (nonexistent)
4310@@ -1,17 +0,0 @@
4311-function [index2,x2,y2,segments2,segmentmarkers2,rifts2] = TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
4312-%TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
4313-%
4314-% Usage:
4315-% [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
4316-%
4317-% (index1,x1,y1,segments1,segmentmarkers1): An initial triangulation.
4318-% [index2,x2,y2,segments2,segmentmarkers2,rifts2]: The resulting triangulation where rifts have been processed.
4319-
4320-% Check usage
4321-if nargin~=5 && nargout~=6
4322- help TriMeshProcessRifts
4323- error('Wrong usage (see above)');
4324-end
4325-
4326-% Call mex module
4327-[index2,x2,y2,segments2,segmentmarkers2,rifts2] = TriMeshProcessRifts_matlab(index1,x1,y1,segments1,segmentmarkers1);
4328Index: ../trunk-jpl/src/m/modules/ProcessRifts.m
4329===================================================================
4330--- ../trunk-jpl/src/m/modules/ProcessRifts.m (nonexistent)
4331+++ ../trunk-jpl/src/m/modules/ProcessRifts.m (revision 22498)
4332@@ -0,0 +1,17 @@
4333+function [index2,x2,y2,segments2,segmentmarkers2,rifts2] = ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
4334+%TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
4335+%
4336+% Usage:
4337+% [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
4338+%
4339+% (index1,x1,y1,segments1,segmentmarkers1): An initial triangulation.
4340+% [index2,x2,y2,segments2,segmentmarkers2,rifts2]: The resulting triangulation where rifts have been processed.
4341+
4342+% Check usage
4343+if nargin~=5 && nargout~=6
4344+ help ProcessRifts
4345+ error('Wrong usage (see above)');
4346+end
4347+
4348+% Call mex module
4349+[index2,x2,y2,segments2,segmentmarkers2,rifts2] = ProcessRifts_matlab(index1,x1,y1,segments1,segmentmarkers1);
4350Index: ../trunk-jpl/src/m/modules/Triangle.py
4351===================================================================
4352--- ../trunk-jpl/src/m/modules/Triangle.py (nonexistent)
4353+++ ../trunk-jpl/src/m/modules/Triangle.py (revision 22498)
4354@@ -0,0 +1,21 @@
4355+from Triangle_python import Triangle_python
4356+
4357+def Triangle(domainoutlinefilename,rifts,mesh_area):
4358+ """
4359+ TRIMESH - Mesh a domain using an .exp file
4360+
4361+ Usage:
4362+ [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,mesh_area);
4363+
4364+ index,x,y: defines a triangulation
4365+ segments: An array made of exterior segments to the mesh domain outline
4366+ segmentmarkers: An array flagging each segment
4367+
4368+ domainoutlinefilename: an Argus domain outline file
4369+ mesh_area: The maximum area desired for any element of the resulting mesh
4370+ """
4371+ # Call mex module
4372+ index,x,y,segments,segmentmarkers=Triangle_python(domainoutlinefilename,rifts,mesh_area)
4373+ # Return
4374+ return index,x,y,segments,segmentmarkers
4375+
4376Index: ../trunk-jpl/src/m/modules/Triangle.m
4377===================================================================
4378--- ../trunk-jpl/src/m/modules/Triangle.m (nonexistent)
4379+++ ../trunk-jpl/src/m/modules/Triangle.m (revision 22498)
4380@@ -0,0 +1,21 @@
4381+function [index,x,y,segments,segmentmarkers] = Triangle(domainoutlinefilename,rifts,mesh_area);
4382+%TRIMESH - Mesh a domain using an .exp file
4383+%
4384+% Usage:
4385+% [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,mesh_area);
4386+%
4387+% index,x,y: Defines a triangulation
4388+% segments: Array made of exterior segments to the mesh domain outline
4389+% segmentmarkers: Array flagging each segment
4390+%
4391+% domainoutlinefilename: Argus domain outline file
4392+% mesh_area: Maximum area desired for any element of the resulting mesh
4393+
4394+% Check usage
4395+if nargin~=3 && nargout~=5
4396+ help Triangle
4397+ error('Wrong usage (see above)');
4398+end
4399+
4400+% Call mex module
4401+[index,x,y,segments,segmentmarkers]=Triangle_matlab(domainoutlinefilename,rifts,mesh_area);
4402Index: ../trunk-jpl/src/m/modules/ProcessRifts.py
4403===================================================================
4404--- ../trunk-jpl/src/m/modules/ProcessRifts.py (nonexistent)
4405+++ ../trunk-jpl/src/m/modules/ProcessRifts.py (revision 22498)
4406@@ -0,0 +1,16 @@
4407+from ProcessRifts_python import ProcessRifts_python
4408+
4409+def ProcessRifts(index1,x1,y1,segments1,segmentmarkers1):
4410+ """
4411+ TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
4412+
4413+ Usage:
4414+ [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
4415+
4416+ (index1,x1,y1,segments1,segmentmarkers1): An initial triangulation.
4417+ [index2,x2,y2,segments2,segmentmarkers2,rifts2]: The resulting triangulation where rifts have been processed.
4418+ """
4419+ # Call mex module
4420+ index2,x2,y2,segments2,segmentmarkers2,rifts2 = ProcessRifts_python(index1,x1,y1,segments1,segmentmarkers1)
4421+ # Return
4422+ return index2,x2,y2,segments2,segmentmarkers2,rifts2
Note: See TracBrowser for help on using the repository browser.