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

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

CHG: added 21724-22754
File size: 150.8 KB

  • ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.cpp

     
    1 /*!\file TriMeshx
    2  * \brief: x code for TriMesh mesher
    3  */
    4 
    5 /*Header files: {{{*/
    6 #include "./TriMeshx.h"
    7 #include "../../shared/shared.h"
    8 #include "../../toolkits/toolkits.h"
    9 /*ANSI_DECLARATORS needed to call triangle library: */
    10 #if defined(_HAVE_TRIANGLE_)
    11         #ifndef ANSI_DECLARATORS
    12         #define ANSI_DECLARATORS
    13         #endif
    14         #include "triangle.h"
    15         #undef ANSI_DECLARATORS
    16 #endif
    17 /*}}}*/
    18 
    19 void TriMeshx(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnsegs,Contours* domain,Contours* rifts,double area){
    20 
    21 #if !defined(_HAVE_TRIANGLE_)
    22         _error_("triangle has not been installed");
    23 #else
    24         /*indexing: */
    25         int i,j;
    26 
    27         /*output: */
    28         int    *index             = NULL;
    29         double *x                 = NULL;
    30         double *y                 = NULL;
    31         int    *segments          = NULL;
    32         int    *segmentmarkerlist = NULL;
    33 
    34         /*intermediary: */
    35         int counter,counter2,backcounter;
    36         Contour<IssmPDouble> *contour = NULL;
    37 
    38         /* Triangle structures needed to call Triangle library routines: */
    39         struct triangulateio in,out;
    40         char   options[256];
    41 
    42         /*Create initial triangulation to call triangulate(). First number of points:*/
    43         in.numberofpoints=0;
    44         for (i=0;i<domain->Size();i++){
    45                 contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    46                 in.numberofpoints+=contour->nods-1;
    47         }
    48         for (i=0;i<rifts->Size();i++){
    49                 contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    50                 in.numberofpoints+=contour->nods;
    51         }
    52 
    53         /*number of point attributes: */
    54         in.numberofpointattributes=1;
    55 
    56         /*fill in the point list: */
    57         in.pointlist = xNew<REAL>(in.numberofpoints*2);
    58 
    59         counter=0;
    60         for (i=0;i<domain->Size();i++){
    61                 contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    62                 for (j=0;j<contour->nods-1;j++){
    63                         in.pointlist[2*counter+0]=contour->x[j];
    64                         in.pointlist[2*counter+1]=contour->y[j];
    65                         counter++;
    66                 }
    67         }
    68         for (i=0;i<rifts->Size();i++){
    69                 contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    70                 for (j=0;j<contour->nods;j++){
    71                         in.pointlist[2*counter+0]=contour->x[j];
    72                         in.pointlist[2*counter+1]=contour->y[j];
    73                         counter++;
    74                 }
    75         }
    76 
    77         /*fill in the point attribute list: */
    78         in.pointattributelist = xNew<REAL>(in.numberofpoints*in.numberofpointattributes);
    79         for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0;
    80 
    81         /*fill in the point marker list: */
    82         in.pointmarkerlist = xNew<int>(in.numberofpoints);
    83         for(i=0;i<in.numberofpoints;i++) in.pointmarkerlist[i] = 0;
    84 
    85         /*Build segments. First figure out number of segments: holes and closed outlines have as many segments as vertices: */
    86         in.numberofsegments=0;
    87         for (i=0;i<domain->Size();i++){
    88                 contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    89                 in.numberofsegments+=contour->nods-1;
    90         }
    91         for(i=0;i<rifts->Size();i++){
    92                 contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    93                 /*for rifts, we have one less segment as we have vertices*/
    94                 in.numberofsegments+=contour->nods-1;
    95         }
    96 
    97         in.segmentlist = xNew<int>(in.numberofsegments*2);
    98         in.segmentmarkerlist = xNewZeroInit<int>(in.numberofsegments);
    99         counter=0;
    100         backcounter=0;
    101         for (i=0;i<domain->Size();i++){
    102                 contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    103                 for (j=0;j<contour->nods-2;j++){
    104                         in.segmentlist[2*counter+0]=counter;
    105                         in.segmentlist[2*counter+1]=counter+1;
    106                         in.segmentmarkerlist[counter]=0;
    107                         counter++;
    108                 }
    109                 /*Close this profile: */
    110                  in.segmentlist[2*counter+0]=counter;
    111                  in.segmentlist[2*counter+1]=backcounter;
    112                  in.segmentmarkerlist[counter]=0;
    113                  counter++;
    114                  backcounter=counter;
    115         }
    116         counter2=counter;
    117         for (i=0;i<rifts->Size();i++){
    118                 contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    119                 for (j=0;j<(contour->nods-1);j++){
    120                         in.segmentlist[2*counter2+0]=counter;
    121                         in.segmentlist[2*counter2+1]=counter+1;
    122                         in.segmentmarkerlist[counter2]=2+i;
    123                         counter2++;
    124                         counter++;
    125                 }
    126                 counter++;
    127         }
    128 
    129         /*Build regions: */
    130         in.numberofregions = 0;
    131 
    132         /*Build holes: */
    133         in.numberofholes = domain->Size()-1; /*everything is a hole, but for the first profile.*/
    134         if(in.numberofholes){
    135                 in.holelist = xNew<REAL>(in.numberofholes*2);
    136                 for (i=0;i<domain->Size()-1;i++){
    137                         contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
    138                         GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
    139                 }
    140         }
    141 
    142         /* Make necessary initializations so that Triangle can return a triangulation in `out': */
    143         out.pointlist             = (REAL*)NULL;
    144         out.pointattributelist    = (REAL*)NULL;
    145         out.pointmarkerlist       = (int *)NULL;
    146         out.trianglelist          = (int *)NULL;
    147         out.triangleattributelist = (REAL*)NULL;
    148         out.neighborlist          = (int *)NULL;
    149         out.segmentlist           = (int *)NULL;
    150         out.segmentmarkerlist     = (int *)NULL;
    151         out.edgelist              = (int *)NULL;
    152         out.edgemarkerlist        = (int *)NULL;
    153 
    154         /* Triangulate the points:.  Switches are chosen to read and write a  */
    155         /*   PSLG (p), preserve the convex hull (c), number everything from  */
    156         /*   zero (z), assign a regional attribute to each element (A), and  */
    157         /*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
    158         /*   neighbor list (n).                                              */
    159         sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
    160         triangulate(options, &in, &out, NULL);
    161         /*report(&out, 0, 1, 1, 1, 1, 0);*/
    162 
    163         /*Allocate index, x and y: */
    164         index=xNew<int>(3*out.numberoftriangles);
    165         x=xNew<double>(out.numberofpoints);
    166         y=xNew<double>(out.numberofpoints);
    167         segments=xNew<int>(3*out.numberofsegments);
    168         segmentmarkerlist=xNew<int>(out.numberofsegments);
    169 
    170         for (i = 0; i< out.numberoftriangles; i++) {
    171                 for (j = 0; j < out.numberofcorners; j++) {
    172                         index[3*i+j]=(int)out.trianglelist[i*out.numberofcorners+j]+1;
    173                 }
    174         }
    175         for (i = 0; i< out.numberofpoints; i++){
    176                 x[i]=(double)out.pointlist[i*2+0];
    177                 y[i]=(double)out.pointlist[i*2+1];
    178         }
    179         for (i = 0; i<out.numberofsegments;i++){
    180                 segments[3*i+0]=(int)out.segmentlist[i*2+0]+1;
    181                 segments[3*i+1]=(int)out.segmentlist[i*2+1]+1;
    182                 segmentmarkerlist[i]=(int)out.segmentmarkerlist[i];
    183         }
    184 
    185         /*Associate elements with segments: */
    186         AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
    187 
    188         /*Order segments so that their normals point outside the domain: */
    189         OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
    190 
    191         /*Output : */
    192         *pindex=index;
    193         *px=x;
    194         *py=y;
    195         *psegments=segments;
    196         *psegmentmarkerlist=segmentmarkerlist;
    197         *pnels=out.numberoftriangles;
    198         *pnods=out.numberofpoints;
    199         *pnsegs=out.numberofsegments;
    200 #endif
    201 }

  • ../trunk-jpl/src/c/modules/TriMeshx/TriMeshx.h

     
    1 /*!\file:  TriMeshx.h
    2  * \brief header file for TriMeshx module
    3  */
    4 
    5 #ifndef _TRIMESHX_H_
    6 #define _TRIMESHX_H_
    7 
    8 #include <string.h>
    9 #include "../../classes/classes.h"
    10 
    11 /* local prototypes: */
    12 void TriMeshx(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnseg,Contours* domain,Contours* rifts,double area);
    13 #endif  /* _TRIMESHX_H */

  • ../trunk-jpl/src/c/modules/TriMeshx

  • ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.h 

    Property changes on: ../trunk-jpl/src/c/modules/TriMeshx
___________________________________________________________________
Deleted: svn:ignore
## -1,2 +0,0 ##
-.deps
-.dirstamp
     
    1 /*!\file:  TriMeshProcessRiftsx.h
    2  * \brief header file for TriMeshProcessRifts module
    3  */
    4 
    5 #ifndef _TRIMESHPROCESSRIFTX_H
    6 #define _TRIMESHPROCESSRIFTX_H
    7 
    8 class RiftStruct;
    9 
    10 void TriMeshProcessRiftsx(int** pindex,int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct);
    11 
    12 #endif  /* _TRIMESHPROCESSRIFTX_H*/

  • ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp

     
    1 /*!\file:  TriMeshProcessRifts.cpp
    2  * \brief split a mesh where a rift (or fault) is present
    3  */
    4 
    5 #include "./TriMeshProcessRiftsx.h"
    6 #include "../../classes/RiftStruct.h"
    7 #include "../../shared/shared.h"
    8 #include "../../toolkits/toolkits.h"
    9 
    10 void TriMeshProcessRiftsx(int** pindex, int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct){
    11 
    12         /*Output*/
    13         int      numrifts,numrifts0;
    14         int     *riftsnumsegments     = NULL;
    15         int    **riftssegments        = NULL;
    16         int     *riftsnumpairs        = NULL;
    17         int    **riftspairs           = NULL;
    18         int     *riftstips            = NULL;
    19         double **riftspenaltypairs    = NULL;
    20         int     *riftsnumpenaltypairs = NULL;
    21 
    22         /*Recover initial mesh*/
    23         int     nel            = *pnel;
    24         int    *index          = *pindex;
    25         double *x              = *px;
    26         double *y              = *py;
    27         int     nods           = *pnods;
    28         int    *segments       = *psegments;
    29         int    *segmentmarkers = *psegmentmarkers;
    30         int     num_seg        = *pnum_seg;
    31 
    32         /*Intermediary*/
    33         int     riftflag;
    34 
    35         /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
    36          *all the nodes of this element belong to the segments (tends to happen when there are corners: */
    37         RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
    38 
    39         /*Figure out if we have rifts, and how many: */
    40         IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
    41 
    42         if(!riftflag) _error_("No rift present in mesh");
    43 
    44         /*Split mesh*/
    45         SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
    46 
    47         /*Order segments so that their normals point outside the domain: */
    48         OrderSegments(&segments,num_seg, index,nel);
    49 
    50         /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
    51          *segmentmarkerlist:*/
    52         SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
    53 
    54         /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
    55         PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
    56 
    57         /*Order rifts so that they start from one tip, go to the other tip, and back: */
    58         OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
    59 
    60         /*Create penalty pairs, used by Imp: */
    61         PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
    62 
    63         /*Create Riftstruct*/
    64         RiftStruct* riftstruct = new RiftStruct(numrifts,riftsnumsegments,riftssegments,riftsnumpairs,riftspairs,riftsnumpenaltypairs,riftspenaltypairs,riftstips);
    65 
    66         /*Assign output pointers for mesh*/
    67         *pnel            = nel;
    68         *pindex          = index;
    69         *px              = x;
    70         *py              = y;
    71         *pnods           = nods;
    72         *psegments       = segments;
    73         *psegmentmarkers = segmentmarkers;
    74         *pnum_seg        = num_seg;
    75 
    76         /*Assign output pointers for rifts*/
    77         *priftstruct = riftstruct;
    78 }

  • ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx

  • ../trunk-jpl/src/c/modules/modules.h 

    Property changes on: ../trunk-jpl/src/c/modules/TriMeshProcessRiftsx
___________________________________________________________________
Deleted: svn:ignore
## -1,2 +0,0 ##
-.deps
-.dirstamp
     
    9090#include "./SystemMatricesx/SystemMatricesx.h"
    9191#include "./SpcNodesx/SpcNodesx.h"
    9292#include "./SurfaceAreax/SurfaceAreax.h"
    93 #include "./TriMeshx/TriMeshx.h"
    94 #include "./TriMeshProcessRiftsx/TriMeshProcessRiftsx.h"
     93#include "./Trianglex/Trianglex.h"
     94#include "./ProcessRiftsx/ProcessRiftsx.h"
    9595#include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
    9696#include "./ThicknessAlongGradientx/ThicknessAlongGradientx.h"
    9797#include "./ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"

  • ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.h

     
     1/*!\file:  ProcessRiftsx.h
     2 * \brief header file for ProcessRifts module
     3 */
     4
     5#ifndef _PROCESSRIFTX_H
     6#define _PROCESSRIFTX_H
     7
     8class RiftStruct;
     9
     10void ProcessRiftsx(int** pindex,int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct);
     11
     12#endif  /* _PROCESSRIFTX_H*/

  • ../trunk-jpl/src/c/modules/ProcessRiftsx/ProcessRiftsx.cpp

     
     1/*!\file:  ProcessRifts.cpp
     2 * \brief split a mesh where a rift (or fault) is present
     3 */
     4
     5#include "./ProcessRiftsx.h"
     6#include "../../classes/RiftStruct.h"
     7#include "../../shared/shared.h"
     8#include "../../toolkits/toolkits.h"
     9
     10void ProcessRiftsx(int** pindex, int* pnel,double** px,double** py,int* pnods,int** psegments,int** psegmentmarkers,int *pnum_seg,RiftStruct **priftstruct){
     11
     12        /*Output*/
     13        int      numrifts,numrifts0;
     14        int     *riftsnumsegments     = NULL;
     15        int    **riftssegments        = NULL;
     16        int     *riftsnumpairs        = NULL;
     17        int    **riftspairs           = NULL;
     18        int     *riftstips            = NULL;
     19        double **riftspenaltypairs    = NULL;
     20        int     *riftsnumpenaltypairs = NULL;
     21
     22        /*Recover initial mesh*/
     23        int     nel            = *pnel;
     24        int    *index          = *pindex;
     25        double *x              = *px;
     26        double *y              = *py;
     27        int     nods           = *pnods;
     28        int    *segments       = *psegments;
     29        int    *segmentmarkers = *psegmentmarkers;
     30        int     num_seg        = *pnum_seg;
     31
     32        /*Intermediary*/
     33        int     riftflag;
     34
     35        /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
     36         *all the nodes of this element belong to the segments (tends to happen when there are corners: */
     37        RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
     38
     39        /*Figure out if we have rifts, and how many: */
     40        IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
     41
     42        if(!riftflag) _error_("No rift present in mesh");
     43
     44        /*Split mesh*/
     45        SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
     46
     47        /*Order segments so that their normals point outside the domain: */
     48        OrderSegments(&segments,num_seg, index,nel);
     49
     50        /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
     51         *segmentmarkerlist:*/
     52        SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
     53
     54        /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
     55        PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
     56
     57        /*Order rifts so that they start from one tip, go to the other tip, and back: */
     58        OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
     59
     60        /*Create penalty pairs, used by Imp: */
     61        PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
     62
     63        /*Create Riftstruct*/
     64        RiftStruct* riftstruct = new RiftStruct(numrifts,riftsnumsegments,riftssegments,riftsnumpairs,riftspairs,riftsnumpenaltypairs,riftspenaltypairs,riftstips);
     65
     66        /*Assign output pointers for mesh*/
     67        *pnel            = nel;
     68        *pindex          = index;
     69        *px              = x;
     70        *py              = y;
     71        *pnods           = nods;
     72        *psegments       = segments;
     73        *psegmentmarkers = segmentmarkers;
     74        *pnum_seg        = num_seg;
     75
     76        /*Assign output pointers for rifts*/
     77        *priftstruct = riftstruct;
     78}

  • ../trunk-jpl/src/c/modules/ProcessRiftsx

  • ../trunk-jpl/src/c/modules/Trianglex/Trianglex.h 

    Property changes on: ../trunk-jpl/src/c/modules/ProcessRiftsx
___________________________________________________________________
Added: svn:ignore
## -0,0 +1,2 ##
+.deps
+.dirstamp
     
     1/*!\file:  Trianglex.h
     2 * \brief header file for Trianglex module
     3 */
     4
     5#ifndef _TRIANGLEX_H_
     6#define _TRIANGLEX_H_
     7
     8#include <string.h>
     9#include "../../classes/classes.h"
     10
     11/* local prototypes: */
     12void Trianglex(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnseg,Contours* domain,Contours* rifts,double area);
     13#endif  /* _TRIANGLEX_H */

  • ../trunk-jpl/src/c/modules/Trianglex/Trianglex.cpp

     
     1/*!\file Trianglex
     2 * \brief: x code for Triangle mesher
     3 */
     4
     5/*Header files: {{{*/
     6#include "./Trianglex.h"
     7#include "../../shared/shared.h"
     8#include "../../toolkits/toolkits.h"
     9/*ANSI_DECLARATORS needed to call triangle library: */
     10#if defined(_HAVE_TRIANGLE_)
     11        #ifndef ANSI_DECLARATORS
     12        #define ANSI_DECLARATORS
     13        #endif
     14        #include "triangle.h"
     15        #undef ANSI_DECLARATORS
     16#endif
     17/*}}}*/
     18
     19void Trianglex(int** pindex,IssmPDouble** px,IssmPDouble** py,int** psegments,int** psegmentmarkerlist,int* pnels,int* pnods, int* pnsegs,Contours* domain,Contours* rifts,double area){
     20
     21#if !defined(_HAVE_TRIANGLE_)
     22        _error_("triangle has not been installed");
     23#else
     24        /*indexing: */
     25        int i,j;
     26
     27        /*output: */
     28        int    *index             = NULL;
     29        double *x                 = NULL;
     30        double *y                 = NULL;
     31        int    *segments          = NULL;
     32        int    *segmentmarkerlist = NULL;
     33
     34        /*intermediary: */
     35        int counter,counter2,backcounter;
     36        Contour<IssmPDouble> *contour = NULL;
     37
     38        /* Triangle structures needed to call Triangle library routines: */
     39        struct triangulateio in,out;
     40        char   options[256];
     41
     42        /*Create initial triangulation to call triangulate(). First number of points:*/
     43        in.numberofpoints=0;
     44        for (i=0;i<domain->Size();i++){
     45                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
     46                in.numberofpoints+=contour->nods-1;
     47        }
     48        for (i=0;i<rifts->Size();i++){
     49                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
     50                in.numberofpoints+=contour->nods;
     51        }
     52
     53        /*number of point attributes: */
     54        in.numberofpointattributes=1;
     55
     56        /*fill in the point list: */
     57        in.pointlist = xNew<REAL>(in.numberofpoints*2);
     58
     59        counter=0;
     60        for (i=0;i<domain->Size();i++){
     61                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
     62                for (j=0;j<contour->nods-1;j++){
     63                        in.pointlist[2*counter+0]=contour->x[j];
     64                        in.pointlist[2*counter+1]=contour->y[j];
     65                        counter++;
     66                }
     67        }
     68        for (i=0;i<rifts->Size();i++){
     69                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
     70                for (j=0;j<contour->nods;j++){
     71                        in.pointlist[2*counter+0]=contour->x[j];
     72                        in.pointlist[2*counter+1]=contour->y[j];
     73                        counter++;
     74                }
     75        }
     76
     77        /*fill in the point attribute list: */
     78        in.pointattributelist = xNew<REAL>(in.numberofpoints*in.numberofpointattributes);
     79        for (i=0;i<in.numberofpoints;i++) in.pointattributelist[i] = 0.0;
     80
     81        /*fill in the point marker list: */
     82        in.pointmarkerlist = xNew<int>(in.numberofpoints);
     83        for(i=0;i<in.numberofpoints;i++) in.pointmarkerlist[i] = 0;
     84
     85        /*Build segments. First figure out number of segments: holes and closed outlines have as many segments as vertices: */
     86        in.numberofsegments=0;
     87        for (i=0;i<domain->Size();i++){
     88                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
     89                in.numberofsegments+=contour->nods-1;
     90        }
     91        for(i=0;i<rifts->Size();i++){
     92                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
     93                /*for rifts, we have one less segment as we have vertices*/
     94                in.numberofsegments+=contour->nods-1;
     95        }
     96
     97        in.segmentlist = xNew<int>(in.numberofsegments*2);
     98        in.segmentmarkerlist = xNewZeroInit<int>(in.numberofsegments);
     99        counter=0;
     100        backcounter=0;
     101        for (i=0;i<domain->Size();i++){
     102                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
     103                for (j=0;j<contour->nods-2;j++){
     104                        in.segmentlist[2*counter+0]=counter;
     105                        in.segmentlist[2*counter+1]=counter+1;
     106                        in.segmentmarkerlist[counter]=0;
     107                        counter++;
     108                }
     109                /*Close this profile: */
     110                 in.segmentlist[2*counter+0]=counter;
     111                 in.segmentlist[2*counter+1]=backcounter;
     112                 in.segmentmarkerlist[counter]=0;
     113                 counter++;
     114                 backcounter=counter;
     115        }
     116        counter2=counter;
     117        for (i=0;i<rifts->Size();i++){
     118                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
     119                for (j=0;j<(contour->nods-1);j++){
     120                        in.segmentlist[2*counter2+0]=counter;
     121                        in.segmentlist[2*counter2+1]=counter+1;
     122                        in.segmentmarkerlist[counter2]=2+i;
     123                        counter2++;
     124                        counter++;
     125                }
     126                counter++;
     127        }
     128
     129        /*Build regions: */
     130        in.numberofregions = 0;
     131
     132        /*Build holes: */
     133        in.numberofholes = domain->Size()-1; /*everything is a hole, but for the first profile.*/
     134        if(in.numberofholes){
     135                in.holelist = xNew<REAL>(in.numberofholes*2);
     136                for (i=0;i<domain->Size()-1;i++){
     137                        contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
     138                        GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
     139                }
     140        }
     141
     142        /* Make necessary initializations so that Triangle can return a triangulation in `out': */
     143        out.pointlist             = (REAL*)NULL;
     144        out.pointattributelist    = (REAL*)NULL;
     145        out.pointmarkerlist       = (int *)NULL;
     146        out.trianglelist          = (int *)NULL;
     147        out.triangleattributelist = (REAL*)NULL;
     148        out.neighborlist          = (int *)NULL;
     149        out.segmentlist           = (int *)NULL;
     150        out.segmentmarkerlist     = (int *)NULL;
     151        out.edgelist              = (int *)NULL;
     152        out.edgemarkerlist        = (int *)NULL;
     153
     154        /* Triangulate the points:.  Switches are chosen to read and write a  */
     155        /*   PSLG (p), preserve the convex hull (c), number everything from  */
     156        /*   zero (z), assign a regional attribute to each element (A), and  */
     157        /*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
     158        /*   neighbor list (n).                                              */
     159        sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
     160        triangulate(options, &in, &out, NULL);
     161        /*report(&out, 0, 1, 1, 1, 1, 0);*/
     162
     163        /*Allocate index, x and y: */
     164        index=xNew<int>(3*out.numberoftriangles);
     165        x=xNew<double>(out.numberofpoints);
     166        y=xNew<double>(out.numberofpoints);
     167        segments=xNew<int>(3*out.numberofsegments);
     168        segmentmarkerlist=xNew<int>(out.numberofsegments);
     169
     170        for (i = 0; i< out.numberoftriangles; i++) {
     171                for (j = 0; j < out.numberofcorners; j++) {
     172                        index[3*i+j]=(int)out.trianglelist[i*out.numberofcorners+j]+1;
     173                }
     174        }
     175        for (i = 0; i< out.numberofpoints; i++){
     176                x[i]=(double)out.pointlist[i*2+0];
     177                y[i]=(double)out.pointlist[i*2+1];
     178        }
     179        for (i = 0; i<out.numberofsegments;i++){
     180                segments[3*i+0]=(int)out.segmentlist[i*2+0]+1;
     181                segments[3*i+1]=(int)out.segmentlist[i*2+1]+1;
     182                segmentmarkerlist[i]=(int)out.segmentmarkerlist[i];
     183        }
     184
     185        /*Associate elements with segments: */
     186        AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
     187
     188        /*Order segments so that their normals point outside the domain: */
     189        OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
     190
     191        /*Output : */
     192        *pindex=index;
     193        *px=x;
     194        *py=y;
     195        *psegments=segments;
     196        *psegmentmarkerlist=segmentmarkerlist;
     197        *pnels=out.numberoftriangles;
     198        *pnods=out.numberofpoints;
     199        *pnsegs=out.numberofsegments;
     200#endif
     201}

  • ../trunk-jpl/src/c/modules/Trianglex

  • ../trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp 

    Property changes on: ../trunk-jpl/src/c/modules/Trianglex
___________________________________________________________________
Added: svn:ignore
## -0,0 +1,2 ##
+.deps
+.dirstamp
     
    1 /*!\file TriMeshx
    2  * \brief: x code for TriMesh mesher
     1/*!\file CoordinateSystemTransformx
     2 * \brief: x code for CoordinateSystemTransformx
    33 */
    44
    55/*Header files*/

  • ../trunk-jpl/src/c/Makefile.am

     
    560560modules_sources= ./shared/Threads/LaunchThread.cpp\
    561561                        ./shared/Threads/PartitionRange.cpp\
    562562                        ./shared/Exp/exp.cpp\
    563                         ./shared/TriMesh/AssociateSegmentToElement.cpp\
    564                         ./shared/TriMesh/GridInsideHole.cpp\
    565                         ./shared/TriMesh/OrderSegments.cpp\
    566                         ./shared/TriMesh/SplitMeshForRifts.cpp\
    567                         ./shared/TriMesh/TriMeshUtils.cpp\
    568                         ./modules/TriMeshx/TriMeshx.cpp\
    569                         ./modules/TriMeshProcessRiftsx/TriMeshProcessRiftsx.cpp\
     563                        ./shared/Triangle/AssociateSegmentToElement.cpp\
     564                        ./shared/Triangle/GridInsideHole.cpp\
     565                        ./shared/Triangle/OrderSegments.cpp\
     566                        ./shared/Triangle/SplitMeshForRifts.cpp\
     567                        ./shared/Triangle/TriangleUtils.cpp\
     568                        ./modules/Trianglex/Trianglex.cpp\
     569                        ./modules/ProcessRiftsx/ProcessRiftsx.cpp\
    570570                        ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp\
    571571                        ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp\
    572572                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\

  • ../trunk-jpl/src/c/main/issm.js

     
    1010        var dbinaryPtr= Module._malloc(nb); var binHeap = new Uint8Array(Module.HEAPU8.buffer,dbinaryPtr,nb);
    1111        binHeap.set(new Uint8Array(dbinary.buffer)); var binary=binHeap.byteOffset;
    1212
    13         //Declare TriMesh module:
     13        //Declare module:
    1414        issm= Module.cwrap('main','number',['number','number']);
    1515       
    1616        //Call issm:

  • ../trunk-jpl/src/c/shared/TriMesh/SplitMeshForRifts.cpp

     
    1 /*
    2  * SplitMeshForRifts.c:
    3  */
    4 #include "./trimesh.h"
    5 #include "../MemOps/MemOps.h"
    6 
    7 int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist){
    8 
    9         /*Some notes on dimensions:
    10         index  of size nelx3
    11         x and y of size nodsx1
    12         segments of size nsegsx3*/
    13 
    14         int i,j,k,l;
    15         int node;
    16         int el;
    17         int  nriftsegs;
    18         int* riftsegments=NULL;
    19         int* flags=NULL;
    20         int  NumGridElementListOnOneSideOfRift;
    21         int* GridElementListOnOneSideOfRift=NULL;
    22 
    23         /*Recover input: */
    24         int     nel               = *pnel;
    25         int    *index             = *pindex;
    26         int     nods              = *pnods;
    27         double *x                 = *px;
    28         double *y                 = *py;
    29         int     nsegs             = *pnsegs;
    30         int    *segments          = *psegments;
    31         int    *segmentmarkerlist = *psegmentmarkerlist;
    32 
    33         /*Establish list of segments that belong to a rift: */
    34         /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
    35         RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);
    36 
    37         /*Go through all nodes of the rift segments, and start splitting the mesh: */
    38         flags=xNewZeroInit<int>(nods); //to make sure we don't split the same nodes twice!
    39         for (i=0;i<nriftsegs;i++){
    40                 for (j=0;j<2;j++){
    41 
    42                         node=riftsegments[4*i+j+2];
    43                         if(flags[node-1]){
    44                                 /*This node was already split, skip:*/
    45                                 continue;
    46                         }
    47                         else{
    48                                 flags[node-1]=1;
    49                         }
    50 
    51                         if(IsGridOnRift(riftsegments,nriftsegs,node)){
    52 
    53                                 DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
    54 
    55                                 /*Summary: we have for node, a list of elements
    56                                  * (GridElementListOnOneSideOfRift, of size
    57                                  * NumGridElementListOnOneSideOfRift) that all contain node
    58                                  *and that are on the same side of the rift. For all these
    59                                  elements, we clone node into another node, and we swap all
    60                                  instances of node in the triangulation *for those elements, to the
    61                                  new node.*/
    62 
    63                                 //create new node
    64                                 x=xReNew<double>(x,nods,nods+1);
    65                                 y=xReNew<double>(y,nods,nods+1);
    66                                 x[nods]=x[node-1]; //matlab indexing
    67                                 y[nods]=y[node-1]; //matlab indexing
    68 
    69                                 //augment number of nodes
    70                                 nods++;
    71 
    72                                 //change elements owning this node
    73                                 for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
    74                                         el=GridElementListOnOneSideOfRift[k];
    75                                         for (l=0;l<3;l++){
    76                                                 if (index[3*el+l]==node) index[3*el+l]=nods; //again, matlab indexing.
    77                                         }
    78                                 }
    79                         }// if(IsGridOnRift(riftsegments,nriftsegs,node))
    80                 } //for(j=0;j<2;j++)
    81         } //for (i=0;i<nriftsegs;i++)
    82 
    83         /*update segments: they got modified completely by adding new nodes.*/
    84         UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
    85 
    86         /*Assign output pointers: */
    87         *pnel=nel;
    88         *pindex=index;
    89         *pnods=nods;
    90         *px=x;
    91         *py=y;
    92         *pnsegs=nsegs;
    93         *psegments=segments;
    94         *psegmentmarkerlist=segmentmarkerlist;
    95         return 1;
    96 }

  • ../trunk-jpl/src/c/shared/TriMesh/GridInsideHole.cpp

     
    1 /*
    2  * GridInsideHole.c:
    3  * from a convex set of points, figure out a point that for sure lies inside the profile.
    4  */
    5 
    6 #include <math.h>
    7 
    8 #include "./trimesh.h"
    9 #include "../Exp/exp.h"
    10 
    11 #undef M_PI
    12 #define M_PI 3.141592653589793238462643
    13 
    14 int GridInsideHole(double* px0,double* py0,int n,double* x,double* y){
    15 
    16         double flag=0.0;
    17         double xA,xB,xC,xD,xE;
    18         double yA,yB,yC,yD,yE;
    19 
    20         /*Take first and last vertices: */
    21         xA=x[0];
    22         yA=y[0];
    23         xB=x[n-1];
    24         yB=y[n-1];
    25 
    26         /*Figure out middle of segment [A B]: */
    27         xC=(xA+xB)/2;
    28         yC=(yA+yB)/2;
    29 
    30         /*D and E are on each side of segment [A B], on the median line between segment [A  B],
    31          *at an angle of 10 degree (less than the minimum 30 enforced by the quality of the mesh: */
    32         xD=xC+tan(10./180.*M_PI)*(yC-yA);
    33         yD=yC+tan(10./180.*M_PI)*(xA-xC);
    34         xE=xC-tan(10./180.*M_PI)*(yC-yA);
    35         yE=yC-tan(10./180.*M_PI)*(xA-xC);
    36 
    37         /*Either E or D is inside profile (x,y): */
    38         IsInPolySerial(&flag,&xD,&yD,1,x,y,n,2);
    39         /*FIXME: used to be 'flag' and not '!flag', check*/
    40         if(!flag){
    41                 /*D is inside the poly: */
    42                 *px0=xD;
    43                 *py0=yD;
    44         }
    45         else{
    46                 /*E is inside the poly: */
    47                 *px0=xE;
    48                 *py0=yE;
    49         }
    50         return 1;
    51 }

  • ../trunk-jpl/src/c/shared/TriMesh/AssociateSegmentToElement.cpp

     
    1 /*!\file:  AssociateSegmentToElement.cpp
    2  * \brief for each segment, look for the corresponding element.
    3  */
    4 
    5 #include "./trimesh.h"
    6 
    7 int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel){
    8 
    9         /*node indices: */
    10         int A,B;
    11 
    12         /*Recover segments: */
    13         int* segments=*psegments;
    14 
    15         for(int i=0;i<nseg;i++){
    16                 A=segments[3*i+0];
    17                 B=segments[3*i+1];
    18                 segments[3*i+2]=FindElement(A,B,index,nel)+1; //matlab indexing.
    19         }
    20 
    21         /*Assign output pointers: */
    22         *psegments=segments;
    23         return 1;
    24 }

  • ../trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp

     
    1 /*
    2  * TriMeshUtils: mesh manipulation routines:
    3  */
    4 
    5 #include <stdio.h>
    6 
    7 #include "./trimesh.h"
    8 #include "../Exceptions/exceptions.h"
    9 #include "../MemOps/MemOps.h"
    10 
    11 #define RIFTPENALTYPAIRSWIDTH 8
    12 int IsGridOnRift(int* riftsegments, int nriftsegs, int node){/*{{{*/
    13 
    14         /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
    15          *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).*/
    16 
    17         int i;
    18         int j;
    19         int count;
    20 
    21         count=0;
    22         for (i=0;i<nriftsegs;i++){
    23                 for (j=0;j<2;j++){
    24                         if ((*(riftsegments+4*i+2+j))==node) count++;
    25                 }
    26         }
    27         if (count==2){
    28                 return 1;
    29         }
    30         else{
    31                 return 0;
    32         }
    33 }/*}}}*/
    34 int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){/*{{{*/
    35 
    36         /*From a node, recover all the elements that are connected to it: */
    37         int i,j;
    38         int noerr=1;
    39 
    40         int max_number_elements=12;
    41         int current_size;
    42         int NumGridElements;
    43         int* GridElements=NULL;
    44         int* GridElementsRealloc=NULL;
    45 
    46         /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own
    47          * the node. We start by allocating GridElements with that size, and realloc
    48          * more if needed.*/
    49 
    50         current_size=max_number_elements;
    51         NumGridElements=0;
    52         GridElements=xNew<int>(max_number_elements);
    53 
    54         for (i=0;i<nel;i++){
    55                 for (j=0;j<3;j++){
    56                         if (index[3*i+j]==node){
    57                                 if (NumGridElements<=(current_size-1)){
    58                                         GridElements[NumGridElements]=i;
    59                                         NumGridElements++;
    60                                         break;
    61                                 }
    62                                 else{
    63                                         /*Reallocate another max_number_elements slots in the GridElements: */
    64                                         GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
    65                                         if (!GridElementsRealloc){
    66                                                 noerr=0;
    67                                                 goto cleanup_and_return;
    68                                         }
    69                                         current_size+=max_number_elements;
    70                                         GridElements=GridElementsRealloc;
    71                                         GridElements[NumGridElements]=i;
    72                                         NumGridElements++;
    73                                         break;
    74                                 }
    75                         }
    76                 }
    77         }
    78         cleanup_and_return:
    79         if(!noerr){
    80                 xDelete<int>(GridElements);
    81         }
    82         /*Allocate return pointers: */
    83         *pGridElements=GridElements;
    84         *pNumGridElements=NumGridElements;
    85         return noerr;
    86 }/*}}}*/
    87 int IsNeighbor(int el1,int el2,int* index){/*{{{*/
    88         /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
    89         int i,j;
    90         int count=0;
    91         for (i=0;i<3;i++){
    92                 for (j=0;j<3;j++){
    93                         if (index[3*el1+i]==index[3*el2+j])count++;
    94                 }
    95         }
    96         if (count==2){
    97                 return 1;
    98         }
    99         else{
    100                 return 0;
    101         }
    102 }/*}}}*/
    103 int IsOnRift(int el,int nriftsegs,int* riftsegments){/*{{{*/
    104         /*From a list of elements segments, figure out if el belongs to it: */
    105         int i;
    106         for (i=0;i<nriftsegs;i++){
    107                 if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
    108                         return 1;
    109                 }
    110         }
    111         return 0;
    112 }/*}}}*/
    113 void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){/*{{{*/
    114 
    115         int i,counter;
    116         int el,el2;
    117 
    118         int  nriftsegs;
    119         int* riftsegments=NULL;
    120         int* riftsegments_uncompressed=NULL;
    121         int  element_nodes[3];
    122 
    123         /*Allocate segmentflags: */
    124         riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
    125 
    126         /*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary
    127          *or a hole: */
    128         nriftsegs=0;
    129         for (i=0;i<nsegs;i++){
    130                 el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
    131                 /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
    132                  *then  proceed to find another element that owns the segment. If we don't find it, we know
    133                  *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
    134                 element_nodes[0]=*(index+3*el+0);
    135                 element_nodes[1]=*(index+3*el+1);
    136                 element_nodes[2]=*(index+3*el+2);
    137 
    138                 index[3*el+0]=-1;
    139                 index[3*el+1]=-1;
    140                 index[3*el+2]=-1;
    141 
    142                 el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
    143 
    144                 /*Restore index: */
    145                 index[3*el+0]=element_nodes[0];
    146                 index[3*el+1]=element_nodes[1];
    147                 index[3*el+2]=element_nodes[2];
    148 
    149                 if (el2!=-1){
    150                         /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
    151                     riftsegments_uncompressed[5*i+0]=1;
    152                     riftsegments_uncompressed[5*i+1]=el;
    153                     riftsegments_uncompressed[5*i+2]=el2;
    154                     riftsegments_uncompressed[5*i+3]=segments[3*i+0];
    155                          riftsegments_uncompressed[5*i+4]=segments[3*i+1];
    156                          nriftsegs++;
    157                 }
    158         }
    159 
    160         /*Compress riftsegments_uncompressed:*/
    161         riftsegments=xNew<int>(nriftsegs*4);
    162         counter=0;
    163         for (i=0;i<nsegs;i++){
    164                 if (riftsegments_uncompressed[5*i+0]){
    165                         riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
    166                         riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
    167                         riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
    168                         riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
    169                         counter++;
    170                 }
    171         }
    172         xDelete<int>(riftsegments_uncompressed);
    173 
    174         /*Assign output pointers: */
    175         *priftsegments=riftsegments;
    176         *pnriftsegs=nriftsegs;
    177 }/*}}}*/
    178 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){/*{{{*/
    179 
    180         int noerr=1;
    181         int k,l,counter;
    182         int newel;
    183 
    184         int* GridElements=NULL;
    185         int  NumGridElements;
    186 
    187         /*Output: */
    188         int NumGridElementListOnOneSideOfRift;
    189         int* GridElementListOnOneSideOfRift=NULL;
    190 
    191         /*Build a list of all the elements connected to this node: */
    192         GridElementsList(&GridElements,&NumGridElements,node,index,nel);
    193 
    194         /*Figure out the list of elements  that are on the same side of the rift. To do so, we start from one
    195          * side of the rift and keep rotating in the same direction:*/
    196         GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
    197         //bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
    198         GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there
    199                                                                                                                            for a rotation direction, we 'll take it out when we are
    200                                                                                                                            done rotating*/
    201         GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
    202         counter=1;
    203         for (;;){
    204                 /*Find neighbour of element GridElementListOnOneSideOfRift[counter], not
    205                  * equal to GridElementListOnOneSideOfRift[counter-1]*/
    206                 for (k=0;k<NumGridElements;k++){
    207                         if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
    208                                 /*Verify this element is not already in our list of element on the same side of the rift: */
    209                                 newel=1;
    210                                 for (l=0;l<=counter;l++){
    211                                         if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
    212                                                 newel=0;
    213                                                 break;
    214                                         }
    215                                 }
    216                                 if (newel){
    217                                         counter++;
    218                                         GridElementListOnOneSideOfRift[counter]=GridElements[k];
    219                                         if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
    220                                                 break;
    221                                         }
    222                                         k=-1;
    223                                 }
    224                         }
    225                 }
    226                 /*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
    227                 NumGridElementListOnOneSideOfRift=counter;
    228                 for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
    229                         GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
    230                 }
    231                 break;
    232         }// for (;;)
    233 
    234         /*Free ressources: */
    235         xDelete<int>(GridElements);
    236         /*Assign output pointers: */
    237         *pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
    238         *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
    239         return noerr;
    240 }/*}}}*/
    241 int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){/*{{{*/
    242 
    243         int noerr=1;
    244         int i,j,k;
    245         int el1,el2;
    246 
    247         int *segments          = NULL;
    248         int *segmentmarkerlist = NULL;
    249         int  nsegs;
    250 
    251         /*Recover input: */
    252         segments          = *psegments;
    253         segmentmarkerlist = *psegmentmarkerlist;
    254         nsegs             = *pnsegs;
    255 
    256         /*Reallocate segments: */
    257         segments         =xReNew<int>(segments,         nsegs*3,(nsegs+nriftsegs)*3);
    258         segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
    259 
    260         /*First, update the existing segments to the new nodes :*/
    261         for (i=0;i<nriftsegs;i++){
    262                 el1=riftsegments[4*i+0];
    263                 el2=riftsegments[4*i+1];
    264                 for (j=0;j<nsegs;j++){
    265                         if (segments[3*j+2]==(el1+1)){
    266                                 /*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
    267                                  *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
    268                                  *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
    269                                 for (k=0;k<3;k++){
    270                                         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])){
    271                                                 *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
    272                                                 break;
    273                                         }
    274                                 }
    275                                 for (k=0;k<3;k++){
    276                                         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])){
    277                                                 *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
    278                                                 break;
    279                                         }
    280                                 }
    281                                 /*Deal with el2: */
    282                                 *(segments+3*(nsegs+i)+2)=el2+1;
    283                                 *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
    284                                 for (k=0;k<3;k++){
    285                                         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])){
    286                                                 *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
    287                                                 break;
    288                                         }
    289                                 }
    290                                 for (k=0;k<3;k++){
    291                                         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])){
    292                                                 *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
    293                                                 break;
    294                                         }
    295                                 }
    296                         }
    297                         if (*(segments+3*j+2)==(el2+1)){
    298                                 /*segment j is the same as rift segment i.*/
    299                                 /*Let's update segments[j][:] using  element el2 and the corresponding rift segment: */
    300                                 for (k=0;k<3;k++){
    301                                         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])){
    302                                                 *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
    303                                                 break;
    304                                         }
    305                                 }
    306                                 for (k=0;k<3;k++){
    307                                         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])){
    308                                                 *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
    309                                                 break;
    310                                         }
    311                                 }
    312                                 /*Deal with el1: */
    313                                 *(segments+3*(nsegs+i)+2)=el1+1;
    314                                 *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
    315                                 for (k=0;k<3;k++){
    316                                         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])){
    317                                                 *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
    318                                                 break;
    319                                         }
    320                                 }
    321                                 for (k=0;k<3;k++){
    322                                         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])){
    323                                                 *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
    324                                                 break;
    325                                         }
    326                                 }
    327                         }
    328                 }
    329         }
    330         nsegs+=nriftsegs;
    331 
    332         /*Assign output pointers: */
    333         *psegments=segments;
    334         *psegmentmarkerlist=segmentmarkerlist;
    335         *pnsegs=nsegs;
    336 
    337         return noerr;
    338 }/*}}}*/
    339 int FindElement(int A,int B,int* index,int nel){/*{{{*/
    340 
    341         int el=-1;
    342         for (int n=0;n<nel;n++){
    343                 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))){
    344                         el=n;
    345                         break;
    346                 }
    347         }
    348         return el;
    349 }/*}}}*/
    350 int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){/*{{{*/
    351 
    352         /*Using segment markers, wring out the rift segments from the segments. Rift markers are
    353          *of the form 2+i where i=0 to number of rifts */
    354 
    355         int noerr=1;
    356         int i,j,counter;
    357 
    358         /*input: */
    359         int *segments          = NULL;
    360         int *segmentmarkerlist = NULL;
    361         int numsegs;
    362 
    363         /*output: */
    364         int   new_numsegs;
    365         int  *riftsnumsegs       = NULL;
    366         int **riftssegments      = NULL;
    367         int  *new_segments       = NULL;
    368         int  *new_segmentmarkers = NULL;
    369 
    370         /*intermediary: */
    371         int* riftsegment=NULL;
    372 
    373         /*Recover input arguments: */
    374         segments          = *psegments;
    375         numsegs           = *pnumsegs;
    376         segmentmarkerlist = *psegmentmarkerlist;
    377 
    378         /*First, figure out  how many segments will be left in 'segments': */
    379         counter=0;
    380         for (i=0;i<numsegs;i++){
    381                 if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
    382         }
    383         /*Allocate new segments: */
    384         new_numsegs=counter;
    385         new_segments=xNew<int>(new_numsegs*3);
    386         new_segmentmarkers=xNew<int>(new_numsegs);
    387 
    388         /*Copy new segments info : */
    389         counter=0;
    390         for (i=0;i<numsegs;i++){
    391                 if (segmentmarkerlist[i]==1){
    392                         new_segments[3*counter+0]=segments[3*i+0];
    393                         new_segments[3*counter+1]=segments[3*i+1];
    394                         new_segments[3*counter+2]=segments[3*i+2];
    395                         new_segmentmarkers[counter]=segmentmarkerlist[i];
    396                         counter++;
    397                 }
    398         }
    399 
    400         /*Now deal with rift segments: */
    401         riftsnumsegs=xNew<int>(numrifts);
    402         riftssegments=xNew<int*>(numrifts);
    403         for (i=0;i<numrifts;i++){
    404                 /*Figure out how many segments for rift i: */
    405                 counter=0;
    406                 for (j=0;j<numsegs;j++){
    407                         if (segmentmarkerlist[j]==2+i)counter++;
    408                 }
    409                 riftsnumsegs[i]=counter;
    410                 riftsegment=xNew<int>(counter*3);
    411                 /*Copy new segments info :*/
    412                 counter=0;
    413                 for (j=0;j<numsegs;j++){
    414                         if (segmentmarkerlist[j]==(2+i)){
    415                                 riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods+1);
    416                                 riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods+1);
    417                                 riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel+1);
    418                                 counter++;
    419                         }
    420                 }
    421                 *(riftssegments+i)=riftsegment;
    422         }
    423 
    424         /*Free ressources: */
    425         xDelete<int>(segments);
    426 
    427         /*Assign output pointers: */
    428         *psegments=new_segments;
    429         *psegmentmarkerlist=new_segmentmarkers;
    430         *pnumsegs=new_numsegs;
    431         *pnumrifts=numrifts;
    432         *priftssegments=riftssegments;
    433         *priftsnumsegs=riftsnumsegs;
    434         return noerr;
    435 }/*}}}*/
    436 int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){/*{{{*/
    437 
    438         int noerr=1;
    439         int i,j,k;
    440 
    441         /*output: */
    442         int  *riftsnumpairs = NULL;
    443         int **riftspairs    = NULL;
    444 
    445         /*intermediary :*/
    446         int  numsegs;
    447         int* segments=NULL;
    448         int* pairs=NULL;
    449         int  node1,node2,node3,node4;
    450 
    451         riftsnumpairs=xNew<int>(numrifts);
    452         riftspairs=xNew<int*>(numrifts);
    453         for (i=0;i<numrifts;i++){
    454                 segments=riftssegments[i];
    455                 numsegs =riftsnumsegments[i];
    456                 riftsnumpairs[i]=numsegs;
    457                 pairs=xNew<int>(2*numsegs);
    458                 for (j=0;j<numsegs;j++){
    459                         pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
    460                         node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
    461                         /*Find element facing on other side of rift: */
    462                         for (k=0;k<numsegs;k++){
    463                                 if (k==j)continue;
    464                                 node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
    465                                 /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
    466                                 if (   ((x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]))
    467                                     || ((x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1]))  ){
    468                                         /*We found the corresponding element: */
    469                                         pairs[2*j+1]=segments[3*k+2];
    470                                         break;
    471                                 }
    472                         }
    473                 }
    474                 riftspairs[i]=pairs;
    475         }
    476 
    477         /*Assign output pointers: */
    478         *priftsnumpairs=riftsnumpairs;
    479         *priftspairs=riftspairs;
    480         return noerr;
    481 }/*}}}*/
    482 int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){/*{{{*/
    483 
    484         int i;
    485         int noerr=1;
    486 
    487         /*output: */
    488         int riftflag=0;
    489         int numrifts=0;
    490 
    491         int maxmark=1; //default marker for regular segments
    492 
    493         /*Any marker >=2 indicates a certain rift: */
    494         numrifts=0;
    495         for (i=0;i<nsegs;i++){
    496                 if (segmentmarkerlist[i]>maxmark){
    497                         numrifts++;
    498                         maxmark=segmentmarkerlist[i];
    499                 }
    500         }
    501         if(numrifts)riftflag=1;
    502 
    503         /*Assign output pointers:*/
    504         *priftflag=riftflag;
    505         *pnumrifts=numrifts;
    506         return noerr;
    507 }/*}}}*/
    508 int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){/*{{{*/
    509 
    510         int noerr=1;
    511         int i,j,k,counter;
    512 
    513         /*intermediary: */
    514         int *riftsegments = NULL;
    515         int *riftpairs    = NULL;
    516         int numsegs;
    517 
    518         /*ordering and copy: */
    519         int *order             = NULL;
    520         int *riftsegments_copy = NULL;
    521         int *riftpairs_copy    = NULL;
    522 
    523         /*node and element manipulation: */
    524         int node1,node2,node3,node4,temp_node,tip1,tip2,node;
    525         int el2;
    526         int already_ordered=0;
    527 
    528         /*output: */
    529         int* riftstips=NULL;
    530 
    531         /*Allocate byproduct of this routine, riftstips: */
    532         riftstips=xNew<int>(numrifts*2);
    533 
    534         /*Go through all rifts: */
    535         for (i=0;i<numrifts;i++){
    536                 riftsegments = riftssegments[i];
    537                 riftpairs    = riftspairs[i];
    538                 numsegs      = riftsnumsegments[i];
    539 
    540                 /*Allocate copy of riftsegments and riftpairs,
    541                  *as well as ordering vector: */
    542                 riftsegments_copy=xNew<int>(numsegs*3);
    543                 riftpairs_copy=xNew<int>(numsegs*2);
    544                 order=xNew<int>(numsegs);
    545 
    546                 /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
    547                 tip1=-1;
    548                 tip2=-1;
    549 
    550                 for (j=0;j<numsegs;j++){
    551                         el2=*(riftpairs+2*j+1);
    552                         node1=*(riftsegments+3*j+0);
    553                         node2=*(riftsegments+3*j+1);
    554                         /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
    555                          *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
    556                         for (k=0;k<numsegs;k++){
    557                                 if (*(riftsegments+3*k+2)==el2){
    558                                         node3=*(riftsegments+3*k+0);
    559                                         node4=*(riftsegments+3*k+1);
    560                                         break;
    561                                 }
    562                         }
    563                         /* Make sure node3 faces node1 and node4 faces node2: */
    564                         _assert_(node1<nods+1 && node4<nods+1);
    565                         _assert_(node1>0 && node4>0);
    566                         if ((x[node1-1]==x[node4-1]) && (y[node1-1]==y[node4-1])){
    567                                 /*Swap node3 and node4:*/
    568                                 temp_node=node3;
    569                                 node3=node4;
    570                                 node4=temp_node;
    571                         }
    572 
    573                         /*Figure out if a tip is on this element: */
    574                         if (node3==node1){
    575                                 /*node1 is a tip*/
    576                                 if (tip1==-1) {
    577                                         tip1=node1;
    578                                         continue;
    579                                 }
    580                                 if ((tip2==-1) && (node1!=tip1)){
    581                                         tip2=node1;
    582                                         break;
    583                                 }
    584                         }
    585 
    586                         if (node4==node2){
    587                                 /*node2 is a tip*/
    588                                 if (tip1==-1){
    589                                         tip1=node2;
    590                                         continue;
    591                                 }
    592                                 if ((tip2==-1) && (node2!=tip1)){
    593                                         tip2=node2;
    594                                         break;
    595                                 }
    596                         }
    597                 }
    598 
    599                 /*Record tips in riftstips: */
    600                 *(riftstips+2*i+0)=tip1;
    601                 *(riftstips+2*i+1)=tip2;
    602 
    603                 /*We have the two tips for this rift.  Go from tip1 to tip2, and figure out the order in which segments are sequential.
    604                  *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
    605                 node=tip1;
    606                 for (counter=0;counter<numsegs;counter++){
    607                         for (j=0;j<numsegs;j++){
    608                                 node1=*(riftsegments+3*j+0);
    609                                 node2=*(riftsegments+3*j+1);
    610 
    611                                 if ((node1==node) || (node2==node)){
    612                                         /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
    613                                         already_ordered=0;
    614                                         for (k=0;k<counter;k++){
    615                                                 if(order[k]==j){
    616                                                         already_ordered=1;
    617                                                         break;
    618                                                 }
    619                                         }
    620                                         if (!already_ordered){
    621                                                 order[counter]=j;
    622                                                 if(node1==node){
    623                                                         node=node2;
    624                                                 }
    625                                                 else if(node2==node){
    626                                                         node=node1;
    627                                                 }
    628                                                 break;
    629                                         }
    630                                 }
    631                         }
    632                 }
    633 
    634                 /*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
    635                 for (j=0;j<numsegs;j++){
    636                         _assert_(order[j]<numsegs);
    637                         *(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
    638                         *(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
    639                         *(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
    640                         *(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
    641                         *(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
    642                 }
    643 
    644                 for (j=0;j<numsegs;j++){
    645                         *(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
    646                         *(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
    647                         *(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
    648                         *(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
    649                         *(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
    650                 }
    651 
    652                 xDelete<int>(order);
    653                 xDelete<int>(riftsegments_copy);
    654                 xDelete<int>(riftpairs_copy);
    655 
    656         }
    657 
    658         /*Assign output pointer:*/
    659         *priftstips=riftstips;
    660         return noerr;
    661 }/*}}}*/
    662 int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,/*{{{*/
    663                 int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
    664 
    665         int noerr=1;
    666         int i,j,k,k0;
    667 
    668         double el1,el2,node1,node2,node3,node4;
    669         double temp_node;
    670 
    671         /*output: */
    672         double **riftspenaltypairs    = NULL;
    673         double  *riftpenaltypairs     = NULL;
    674         int     *riftsnumpenaltypairs = NULL;
    675 
    676         /*intermediary: */
    677         int numsegs;
    678         int* riftsegments=NULL;
    679         int* riftpairs=NULL;
    680         int counter;
    681         double normal[2];
    682         double length;
    683         int    k1,k2;
    684 
    685         /*Allocate: */
    686         riftspenaltypairs=xNew<double*>(numrifts);
    687         riftsnumpenaltypairs=xNew<int>(numrifts);
    688 
    689         for(i=0;i<numrifts;i++){
    690                 numsegs=riftsnumsegs[i];
    691                 riftsegments=riftssegments[i];
    692                 riftpairs=riftspairs[i];
    693 
    694                 /*allocate riftpenaltypairs, and riftnumpenaltypairs: */
    695                 if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
    696 
    697                 /*Go through only one flank of the rifts, not counting the tips: */
    698                 counter=0;
    699                 for(j=0;j<(numsegs/2);j++){
    700                         el1=*(riftpairs+2*j+0);
    701                         el2=*(riftpairs+2*j+1);
    702                         node1=*(riftsegments+3*j+0);
    703                         node2=*(riftsegments+3*j+1);
    704                         /*Find segment index to recover node3 and node4, facing node1 and node2: */
    705                         k0=-1;
    706                         for(k=0;k<numsegs;k++){
    707                                 if(*(riftsegments+3*k+2)==el2){
    708                                         k0=k;
    709                                         break;
    710                                 }
    711                         }
    712                         node3=*(riftsegments+3*k0+0);
    713                         node4=*(riftsegments+3*k0+1);
    714 
    715                         /* Make sure node3 faces node1 and node4 faces node2: */
    716                         if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
    717                                 /*Swap node3 and node4:*/
    718                                 temp_node=node3;
    719                                 node3=node4;
    720                                 node4=temp_node;
    721                         }       
    722                         /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
    723                          *this segment, and its length: */
    724                         normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
    725                         normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
    726                         length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
    727 
    728                         /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
    729                          * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4,
    730                          * only once. We'll add the normals and the lengths : */
    731 
    732                         if(node1!=node3){ //exclude tips from loads
    733                                 k1=-1;
    734                                 for(k=0;k<counter;k++){
    735                                         if( (*(riftpenaltypairs+k*7+0))==node1){
    736                                                 k1=k;
    737                                                 break;
    738                                         }
    739                                 }
    740                                 if(k1==-1){
    741                                         *(riftpenaltypairs+counter*7+0)=node1;
    742                                         *(riftpenaltypairs+counter*7+1)=node3;
    743                                         *(riftpenaltypairs+counter*7+2)=el1;
    744                                         *(riftpenaltypairs+counter*7+3)=el2;
    745                                         *(riftpenaltypairs+counter*7+4)=normal[0];
    746                                         *(riftpenaltypairs+counter*7+5)=normal[1];
    747                                         *(riftpenaltypairs+counter*7+6)=length/2;
    748                                         counter++;
    749                                 }
    750                                 else{
    751                                         *(riftpenaltypairs+k1*7+4)+=normal[0];
    752                                         *(riftpenaltypairs+k1*7+5)+=normal[1];
    753                                         *(riftpenaltypairs+k1*7+6)+=length/2;
    754                                 }
    755                         }
    756                         if(node2!=node4){
    757                                 k2=-1;
    758                                 for(k=0;k<counter;k++){
    759                                         if( (*(riftpenaltypairs+k*7+0))==node2){
    760                                                 k2=k;
    761                                                 break;
    762                                         }
    763                                 }
    764                                 if(k2==-1){
    765                                         *(riftpenaltypairs+counter*7+0)=node2;
    766                                         *(riftpenaltypairs+counter*7+1)=node4;
    767                                         *(riftpenaltypairs+counter*7+2)=el1;
    768                                         *(riftpenaltypairs+counter*7+3)=el2;
    769                                         *(riftpenaltypairs+counter*7+4)=normal[0];
    770                                         *(riftpenaltypairs+counter*7+5)=normal[1];
    771                                         *(riftpenaltypairs+counter*7+6)=length/2;
    772                                         counter++;
    773                                 }
    774                                 else{
    775                                         *(riftpenaltypairs+k2*7+4)+=normal[0];
    776                                         *(riftpenaltypairs+k2*7+5)+=normal[1];
    777                                         *(riftpenaltypairs+k2*7+6)+=length/2;
    778                                 }
    779                         }
    780                 }
    781                 /*Renormalize normals: */
    782                 for(j=0;j<counter;j++){
    783                         double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
    784                         *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
    785                         *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
    786                 }
    787 
    788                 riftspenaltypairs[i]=riftpenaltypairs;
    789                 riftsnumpenaltypairs[i]=(numsegs/2-1);
    790         }
    791 
    792         /*Assign output pointers: */
    793         *priftspenaltypairs=riftspenaltypairs;
    794         *priftsnumpenaltypairs=riftsnumpenaltypairs;
    795         return noerr;
    796 }/*}}}*/
    797 int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
    798 
    799         int noerr=1;
    800         int i,j,k;
    801         int node1,node2,node3;
    802         int el;
    803         double  pair[2];
    804         int     pair_count=0;
    805         int     triple=0;
    806 
    807         /*Recover input: */
    808         int    *index = *pindex;
    809         int     nel   = *pnel;
    810         double *x     = *px;
    811         double *y     = *py;
    812         int     nods  = *pnods;
    813 
    814         for (i=0;i<num_seg;i++){
    815                 node1=*(segments+3*i+0);
    816                 node2=*(segments+3*i+1);
    817                 /*Find all elements connected to [node1 node2]: */
    818                 pair_count=0;
    819                 for (j=0;j<nel;j++){
    820                         if (*(index+3*j+0)==node1){
    821                                 if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
    822                                         pair[pair_count]=j;
    823                                         pair_count++;
    824                                 }
    825                         }
    826                         if (*(index+3*j+1)==node1){
    827                                 if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
    828                                         pair[pair_count]=j;
    829                                         pair_count++;
    830                                 }
    831                         }
    832                         if (*(index+3*j+2)==node1){
    833                                 if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
    834                                         pair[pair_count]=j;
    835                                         pair_count++;
    836                                 }
    837                         }
    838                 }
    839                 /*Ok, we have pair_count elements connected to this segment. For each of these elements,
    840                  *figure out if the third node also belongs to a segment: */
    841                 if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to  2 elements
    842                         continue;
    843                 }
    844                 else{
    845                         for (j=0;j<pair_count;j++){
    846                                 el=(int)pair[j];
    847                                 triple=0;
    848                                 /*First find node3: */
    849                                 if (*(index+3*el+0)==node1){
    850                                         if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
    851                                         else node3=*(index+3*el+1);
    852                                 }
    853                                 if (*(index+3*el+1)==node1){
    854                                         if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
    855                                         else node3=*(index+3*el+0);
    856                                 }
    857                                 if (*(index+3*el+2)==node1){
    858                                         if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
    859                                         else node3=*(index+3*el+0);
    860                                 }
    861                                 /*Ok, we have node3. Does node3 belong to a segment? : */
    862                                 for (k=0;k<num_seg;k++){
    863                                         if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
    864                                                 triple=1;
    865                                                 break;
    866                                         }
    867                                 }
    868                                 if(triple==1){
    869                                         /*el is a corner element: we need to split it in 3 triangles: */
    870                                         x=xReNew<double>(x,nods,nods+1);
    871                                         y=xReNew<double>(y,nods,nods+1);
    872                                         x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
    873                                         y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
    874                                         index=xReNew<int>(index,nel*3,(nel+2*3));
    875                                         /*First, reassign element el: */
    876                                         *(index+3*el+0)=node1;
    877                                         *(index+3*el+1)=node2;
    878                                         *(index+3*el+2)=nods+1;
    879                                         /*Other two elements: */
    880                                         *(index+3*nel+0)=node2;
    881                                         *(index+3*nel+1)=node3;
    882                                         *(index+3*nel+2)=nods+1;
    883 
    884                                         *(index+3*(nel+1)+0)=node3;
    885                                         *(index+3*(nel+1)+1)=node1;
    886                                         *(index+3*(nel+1)+2)=nods+1;
    887                                         /*we need  to change the segment elements corresponding to el: */
    888                                         for (k=0;k<num_seg;k++){
    889                                                 if (*(segments+3*k+2)==(el+1)){
    890                                                         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;
    891                                                         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;
    892                                                         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;
    893                                                 }
    894                                         }
    895 
    896                                         nods=nods+1;
    897                                         nel=nel+2;
    898                                         i=0;
    899                                         break;
    900                                 }
    901                         } //for (j=0;j<pair_count;j++)
    902                 }
    903         }// for (i=0;i<num_seg;i++)
    904 
    905         /*Assign output pointers: */
    906         *pindex=index;
    907         *pnel=nel;
    908         *px=x;
    909         *py=y;
    910         *pnods=nods;
    911         return noerr;
    912 }/*}}}*/

  • ../trunk-jpl/src/c/shared/TriMesh/trimesh.h

     
    1 /*!\file:  trimesh.h
    2  * \brief
    3  */
    4 
    5 #ifndef _SHARED_TRIMESH_H
    6 #define _SHARED_TRIMESH_H
    7 
    8 #include <stdio.h>
    9 #include <math.h>
    10 
    11 //#define REAL double //took  it out because it may conflict with stdlib.h defines. put back if necessary
    12 int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel);
    13 int OrderSegments(int** psegments,int nseg, int* index,int nel);
    14 int GridInsideHole(double* px0,double* py0,int n,double* x,double* y);
    15 int FindElement(int A,int B,int* index,int nel);
    16 int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist);
    17 int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
    18 int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
    19 int IsNeighbor(int el1,int el2,int* index);
    20 int IsOnRift(int el,int nriftsegs,int* riftsegments);
    21 void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index, int nsegs,int* segments);
    22 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift,int segmentnumber, int nriftsegs,int* riftsegments, int node,int* index,int nel);
    23 int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
    24 int FindElement(double A,double B,int* index,int nel);
    25 int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
    26 int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels);
    27 int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels);
    28 int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int**  riftssegments,
    29                 int* riftsnumsegments,int** riftspairs,int* riftstips,double* x,double* y);
    30 int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg);
    31 int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y);
    32 
    33 #endif  /* _SHARED_TRIMESH_H */

  • ../trunk-jpl/src/c/shared/TriMesh/OrderSegments.cpp

     
    1 /*
    2  * OrderSegments.c:
    3  * reorder segments so that their normals point outside the domain outline.
    4  */
    5 #include "./trimesh.h"
    6 
    7 int OrderSegments(int** psegments,int nseg,int* index,int nel){
    8 
    9         /*vertex indices: */
    10         int A,B;
    11 
    12         /*element index*/
    13         int el;
    14 
    15         /*Recover segments: */
    16         int* segments=*psegments;
    17 
    18         for(int i=0;i<nseg;i++){
    19                 A=segments[3*i+0];
    20                 B=segments[3*i+1];
    21                 el=segments[3*i+2]-1; //after AssociateSegmentToElement, el was a matlab index, we need the c index now.
    22 
    23                 if (index[3*el+0]==A){
    24                         if (index[3*el+2]==B){
    25                                 segments[3*i+0]=B;
    26                                 segments[3*i+1]=A;
    27                         }
    28                 }
    29                 else if (index[3*el+1]==A){
    30                         if (index[3*el+0]==B){
    31                                 segments[3*i+0]=B;
    32                                 segments[3*i+1]=A;
    33                         }
    34                 }
    35                 else{
    36                         if (index[3*el+1]==B){
    37                                 segments[3*i+0]=B;
    38                                 segments[3*i+1]=A;
    39                         }
    40                 }
    41         }
    42 
    43         /*Assign output pointers: */
    44         *psegments=segments;
    45         return 1;
    46 }

  • ../trunk-jpl/src/c/shared/TriMesh

  • ../trunk-jpl/src/c/shared/shared.h 

    Property changes on: ../trunk-jpl/src/c/shared/TriMesh
___________________________________________________________________
Deleted: svn:ignore
## -1,2 +0,0 ##
-.deps
-.dirstamp
     
    1818#include "./Sorting/sorting.h"
    1919#include "./String/sharedstring.h"
    2020#include "./Threads/issm_threads.h"
    21 #include "./TriMesh/trimesh.h"
     21#include "./Triangle/triangle.h"
    2222#include "./LatLong/latlong.h"
    2323
    2424#endif

  • ../trunk-jpl/src/c/shared/Triangle/OrderSegments.cpp

     
     1/*
     2 * OrderSegments.c:
     3 * reorder segments so that their normals point outside the domain outline.
     4 */
     5#include "./trimesh.h"
     6
     7int OrderSegments(int** psegments,int nseg,int* index,int nel){
     8
     9        /*vertex indices: */
     10        int A,B;
     11
     12        /*element index*/
     13        int el;
     14
     15        /*Recover segments: */
     16        int* segments=*psegments;
     17
     18        for(int i=0;i<nseg;i++){
     19                A=segments[3*i+0];
     20                B=segments[3*i+1];
     21                el=segments[3*i+2]-1; //after AssociateSegmentToElement, el was a matlab index, we need the c index now.
     22
     23                if (index[3*el+0]==A){
     24                        if (index[3*el+2]==B){
     25                                segments[3*i+0]=B;
     26                                segments[3*i+1]=A;
     27                        }
     28                }
     29                else if (index[3*el+1]==A){
     30                        if (index[3*el+0]==B){
     31                                segments[3*i+0]=B;
     32                                segments[3*i+1]=A;
     33                        }
     34                }
     35                else{
     36                        if (index[3*el+1]==B){
     37                                segments[3*i+0]=B;
     38                                segments[3*i+1]=A;
     39                        }
     40                }
     41        }
     42
     43        /*Assign output pointers: */
     44        *psegments=segments;
     45        return 1;
     46}

  • ../trunk-jpl/src/c/shared/Triangle/TriangleUtils.cpp

     
     1/*
     2 * TriangleUtils: mesh manipulation routines:
     3 */
     4
     5#include <stdio.h>
     6
     7#include "./triangle.h"
     8#include "../Exceptions/exceptions.h"
     9#include "../MemOps/MemOps.h"
     10
     11#define RIFTPENALTYPAIRSWIDTH 8
     12int IsGridOnRift(int* riftsegments, int nriftsegs, int node){/*{{{*/
     13
     14        /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
     15         *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).*/
     16
     17        int i;
     18        int j;
     19        int count;
     20
     21        count=0;
     22        for (i=0;i<nriftsegs;i++){
     23                for (j=0;j<2;j++){
     24                        if ((*(riftsegments+4*i+2+j))==node) count++;
     25                }
     26        }
     27        if (count==2){
     28                return 1;
     29        }
     30        else{
     31                return 0;
     32        }
     33}/*}}}*/
     34int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){/*{{{*/
     35
     36        /*From a node, recover all the elements that are connected to it: */
     37        int i,j;
     38        int noerr=1;
     39
     40        int max_number_elements=12;
     41        int current_size;
     42        int NumGridElements;
     43        int* GridElements=NULL;
     44        int* GridElementsRealloc=NULL;
     45
     46        /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own
     47         * the node. We start by allocating GridElements with that size, and realloc
     48         * more if needed.*/
     49
     50        current_size=max_number_elements;
     51        NumGridElements=0;
     52        GridElements=xNew<int>(max_number_elements);
     53
     54        for (i=0;i<nel;i++){
     55                for (j=0;j<3;j++){
     56                        if (index[3*i+j]==node){
     57                                if (NumGridElements<=(current_size-1)){
     58                                        GridElements[NumGridElements]=i;
     59                                        NumGridElements++;
     60                                        break;
     61                                }
     62                                else{
     63                                        /*Reallocate another max_number_elements slots in the GridElements: */
     64                                        GridElementsRealloc=xReNew<int>(GridElements,current_size,(current_size+max_number_elements));
     65                                        if (!GridElementsRealloc){
     66                                                noerr=0;
     67                                                goto cleanup_and_return;
     68                                        }
     69                                        current_size+=max_number_elements;
     70                                        GridElements=GridElementsRealloc;
     71                                        GridElements[NumGridElements]=i;
     72                                        NumGridElements++;
     73                                        break;
     74                                }
     75                        }
     76                }
     77        }
     78        cleanup_and_return:
     79        if(!noerr){
     80                xDelete<int>(GridElements);
     81        }
     82        /*Allocate return pointers: */
     83        *pGridElements=GridElements;
     84        *pNumGridElements=NumGridElements;
     85        return noerr;
     86}/*}}}*/
     87int IsNeighbor(int el1,int el2,int* index){/*{{{*/
     88        /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
     89        int i,j;
     90        int count=0;
     91        for (i=0;i<3;i++){
     92                for (j=0;j<3;j++){
     93                        if (index[3*el1+i]==index[3*el2+j])count++;
     94                }
     95        }
     96        if (count==2){
     97                return 1;
     98        }
     99        else{
     100                return 0;
     101        }
     102}/*}}}*/
     103int IsOnRift(int el,int nriftsegs,int* riftsegments){/*{{{*/
     104        /*From a list of elements segments, figure out if el belongs to it: */
     105        int i;
     106        for (i=0;i<nriftsegs;i++){
     107                if ((*(riftsegments+4*i+0)==el) || (*(riftsegments+4*i+1)==el)){
     108                        return 1;
     109                }
     110        }
     111        return 0;
     112}/*}}}*/
     113void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){/*{{{*/
     114
     115        int i,counter;
     116        int el,el2;
     117
     118        int  nriftsegs;
     119        int* riftsegments=NULL;
     120        int* riftsegments_uncompressed=NULL;
     121        int  element_nodes[3];
     122
     123        /*Allocate segmentflags: */
     124        riftsegments_uncompressed=xNewZeroInit<int>(nsegs*5);
     125
     126        /*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary
     127         *or a hole: */
     128        nriftsegs=0;
     129        for (i=0;i<nsegs;i++){
     130                el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
     131                /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
     132                 *then  proceed to find another element that owns the segment. If we don't find it, we know
     133                 *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
     134                element_nodes[0]=*(index+3*el+0);
     135                element_nodes[1]=*(index+3*el+1);
     136                element_nodes[2]=*(index+3*el+2);
     137
     138                index[3*el+0]=-1;
     139                index[3*el+1]=-1;
     140                index[3*el+2]=-1;
     141
     142                el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
     143
     144                /*Restore index: */
     145                index[3*el+0]=element_nodes[0];
     146                index[3*el+1]=element_nodes[1];
     147                index[3*el+2]=element_nodes[2];
     148
     149                if (el2!=-1){
     150                        /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
     151                    riftsegments_uncompressed[5*i+0]=1;
     152                    riftsegments_uncompressed[5*i+1]=el;
     153                    riftsegments_uncompressed[5*i+2]=el2;
     154                    riftsegments_uncompressed[5*i+3]=segments[3*i+0];
     155                         riftsegments_uncompressed[5*i+4]=segments[3*i+1];
     156                         nriftsegs++;
     157                }
     158        }
     159
     160        /*Compress riftsegments_uncompressed:*/
     161        riftsegments=xNew<int>(nriftsegs*4);
     162        counter=0;
     163        for (i=0;i<nsegs;i++){
     164                if (riftsegments_uncompressed[5*i+0]){
     165                        riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
     166                        riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
     167                        riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
     168                        riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
     169                        counter++;
     170                }
     171        }
     172        xDelete<int>(riftsegments_uncompressed);
     173
     174        /*Assign output pointers: */
     175        *priftsegments=riftsegments;
     176        *pnriftsegs=nriftsegs;
     177}/*}}}*/
     178int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){/*{{{*/
     179
     180        int noerr=1;
     181        int k,l,counter;
     182        int newel;
     183
     184        int* GridElements=NULL;
     185        int  NumGridElements;
     186
     187        /*Output: */
     188        int NumGridElementListOnOneSideOfRift;
     189        int* GridElementListOnOneSideOfRift=NULL;
     190
     191        /*Build a list of all the elements connected to this node: */
     192        GridElementsList(&GridElements,&NumGridElements,node,index,nel);
     193
     194        /*Figure out the list of elements  that are on the same side of the rift. To do so, we start from one
     195         * side of the rift and keep rotating in the same direction:*/
     196        GridElementListOnOneSideOfRift=xNew<int>(NumGridElements);
     197        //bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
     198        GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there
     199                                                                                                                           for a rotation direction, we 'll take it out when we are
     200                                                                                                                           done rotating*/
     201        GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
     202        counter=1;
     203        for (;;){
     204                /*Find neighbour of element GridElementListOnOneSideOfRift[counter], not
     205                 * equal to GridElementListOnOneSideOfRift[counter-1]*/
     206                for (k=0;k<NumGridElements;k++){
     207                        if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
     208                                /*Verify this element is not already in our list of element on the same side of the rift: */
     209                                newel=1;
     210                                for (l=0;l<=counter;l++){
     211                                        if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
     212                                                newel=0;
     213                                                break;
     214                                        }
     215                                }
     216                                if (newel){
     217                                        counter++;
     218                                        GridElementListOnOneSideOfRift[counter]=GridElements[k];
     219                                        if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
     220                                                break;
     221                                        }
     222                                        k=-1;
     223                                }
     224                        }
     225                }
     226                /*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
     227                NumGridElementListOnOneSideOfRift=counter;
     228                for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
     229                        GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
     230                }
     231                break;
     232        }// for (;;)
     233
     234        /*Free ressources: */
     235        xDelete<int>(GridElements);
     236        /*Assign output pointers: */
     237        *pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
     238        *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
     239        return noerr;
     240}/*}}}*/
     241int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){/*{{{*/
     242
     243        int noerr=1;
     244        int i,j,k;
     245        int el1,el2;
     246
     247        int *segments          = NULL;
     248        int *segmentmarkerlist = NULL;
     249        int  nsegs;
     250
     251        /*Recover input: */
     252        segments          = *psegments;
     253        segmentmarkerlist = *psegmentmarkerlist;
     254        nsegs             = *pnsegs;
     255
     256        /*Reallocate segments: */
     257        segments         =xReNew<int>(segments,         nsegs*3,(nsegs+nriftsegs)*3);
     258        segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
     259
     260        /*First, update the existing segments to the new nodes :*/
     261        for (i=0;i<nriftsegs;i++){
     262                el1=riftsegments[4*i+0];
     263                el2=riftsegments[4*i+1];
     264                for (j=0;j<nsegs;j++){
     265                        if (segments[3*j+2]==(el1+1)){
     266                                /*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
     267                                 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
     268                                 *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
     269                                for (k=0;k<3;k++){
     270                                        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])){
     271                                                *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
     272                                                break;
     273                                        }
     274                                }
     275                                for (k=0;k<3;k++){
     276                                        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])){
     277                                                *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
     278                                                break;
     279                                        }
     280                                }
     281                                /*Deal with el2: */
     282                                *(segments+3*(nsegs+i)+2)=el2+1;
     283                                *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
     284                                for (k=0;k<3;k++){
     285                                        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])){
     286                                                *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
     287                                                break;
     288                                        }
     289                                }
     290                                for (k=0;k<3;k++){
     291                                        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])){
     292                                                *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
     293                                                break;
     294                                        }
     295                                }
     296                        }
     297                        if (*(segments+3*j+2)==(el2+1)){
     298                                /*segment j is the same as rift segment i.*/
     299                                /*Let's update segments[j][:] using  element el2 and the corresponding rift segment: */
     300                                for (k=0;k<3;k++){
     301                                        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])){
     302                                                *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
     303                                                break;
     304                                        }
     305                                }
     306                                for (k=0;k<3;k++){
     307                                        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])){
     308                                                *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
     309                                                break;
     310                                        }
     311                                }
     312                                /*Deal with el1: */
     313                                *(segments+3*(nsegs+i)+2)=el1+1;
     314                                *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
     315                                for (k=0;k<3;k++){
     316                                        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])){
     317                                                *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
     318                                                break;
     319                                        }
     320                                }
     321                                for (k=0;k<3;k++){
     322                                        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])){
     323                                                *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
     324                                                break;
     325                                        }
     326                                }
     327                        }
     328                }
     329        }
     330        nsegs+=nriftsegs;
     331
     332        /*Assign output pointers: */
     333        *psegments=segments;
     334        *psegmentmarkerlist=segmentmarkerlist;
     335        *pnsegs=nsegs;
     336
     337        return noerr;
     338}/*}}}*/
     339int FindElement(int A,int B,int* index,int nel){/*{{{*/
     340
     341        int el=-1;
     342        for (int n=0;n<nel;n++){
     343                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))){
     344                        el=n;
     345                        break;
     346                }
     347        }
     348        return el;
     349}/*}}}*/
     350int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){/*{{{*/
     351
     352        /*Using segment markers, wring out the rift segments from the segments. Rift markers are
     353         *of the form 2+i where i=0 to number of rifts */
     354
     355        int noerr=1;
     356        int i,j,counter;
     357
     358        /*input: */
     359        int *segments          = NULL;
     360        int *segmentmarkerlist = NULL;
     361        int numsegs;
     362
     363        /*output: */
     364        int   new_numsegs;
     365        int  *riftsnumsegs       = NULL;
     366        int **riftssegments      = NULL;
     367        int  *new_segments       = NULL;
     368        int  *new_segmentmarkers = NULL;
     369
     370        /*intermediary: */
     371        int* riftsegment=NULL;
     372
     373        /*Recover input arguments: */
     374        segments          = *psegments;
     375        numsegs           = *pnumsegs;
     376        segmentmarkerlist = *psegmentmarkerlist;
     377
     378        /*First, figure out  how many segments will be left in 'segments': */
     379        counter=0;
     380        for (i=0;i<numsegs;i++){
     381                if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
     382        }
     383        /*Allocate new segments: */
     384        new_numsegs=counter;
     385        new_segments=xNew<int>(new_numsegs*3);
     386        new_segmentmarkers=xNew<int>(new_numsegs);
     387
     388        /*Copy new segments info : */
     389        counter=0;
     390        for (i=0;i<numsegs;i++){
     391                if (segmentmarkerlist[i]==1){
     392                        new_segments[3*counter+0]=segments[3*i+0];
     393                        new_segments[3*counter+1]=segments[3*i+1];
     394                        new_segments[3*counter+2]=segments[3*i+2];
     395                        new_segmentmarkers[counter]=segmentmarkerlist[i];
     396                        counter++;
     397                }
     398        }
     399
     400        /*Now deal with rift segments: */
     401        riftsnumsegs=xNew<int>(numrifts);
     402        riftssegments=xNew<int*>(numrifts);
     403        for (i=0;i<numrifts;i++){
     404                /*Figure out how many segments for rift i: */
     405                counter=0;
     406                for (j=0;j<numsegs;j++){
     407                        if (segmentmarkerlist[j]==2+i)counter++;
     408                }
     409                riftsnumsegs[i]=counter;
     410                riftsegment=xNew<int>(counter*3);
     411                /*Copy new segments info :*/
     412                counter=0;
     413                for (j=0;j<numsegs;j++){
     414                        if (segmentmarkerlist[j]==(2+i)){
     415                                riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods+1);
     416                                riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods+1);
     417                                riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel+1);
     418                                counter++;
     419                        }
     420                }
     421                *(riftssegments+i)=riftsegment;
     422        }
     423
     424        /*Free ressources: */
     425        xDelete<int>(segments);
     426
     427        /*Assign output pointers: */
     428        *psegments=new_segments;
     429        *psegmentmarkerlist=new_segmentmarkers;
     430        *pnumsegs=new_numsegs;
     431        *pnumrifts=numrifts;
     432        *priftssegments=riftssegments;
     433        *priftsnumsegs=riftsnumsegs;
     434        return noerr;
     435}/*}}}*/
     436int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){/*{{{*/
     437
     438        int noerr=1;
     439        int i,j,k;
     440
     441        /*output: */
     442        int  *riftsnumpairs = NULL;
     443        int **riftspairs    = NULL;
     444
     445        /*intermediary :*/
     446        int  numsegs;
     447        int* segments=NULL;
     448        int* pairs=NULL;
     449        int  node1,node2,node3,node4;
     450
     451        riftsnumpairs=xNew<int>(numrifts);
     452        riftspairs=xNew<int*>(numrifts);
     453        for (i=0;i<numrifts;i++){
     454                segments=riftssegments[i];
     455                numsegs =riftsnumsegments[i];
     456                riftsnumpairs[i]=numsegs;
     457                pairs=xNew<int>(2*numsegs);
     458                for (j=0;j<numsegs;j++){
     459                        pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
     460                        node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
     461                        /*Find element facing on other side of rift: */
     462                        for (k=0;k<numsegs;k++){
     463                                if (k==j)continue;
     464                                node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
     465                                /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
     466                                if (   ((x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2]))
     467                                    || ((x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1]))  ){
     468                                        /*We found the corresponding element: */
     469                                        pairs[2*j+1]=segments[3*k+2];
     470                                        break;
     471                                }
     472                        }
     473                }
     474                riftspairs[i]=pairs;
     475        }
     476
     477        /*Assign output pointers: */
     478        *priftsnumpairs=riftsnumpairs;
     479        *priftspairs=riftspairs;
     480        return noerr;
     481}/*}}}*/
     482int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){/*{{{*/
     483
     484        int i;
     485        int noerr=1;
     486
     487        /*output: */
     488        int riftflag=0;
     489        int numrifts=0;
     490
     491        int maxmark=1; //default marker for regular segments
     492
     493        /*Any marker >=2 indicates a certain rift: */
     494        numrifts=0;
     495        for (i=0;i<nsegs;i++){
     496                if (segmentmarkerlist[i]>maxmark){
     497                        numrifts++;
     498                        maxmark=segmentmarkerlist[i];
     499                }
     500        }
     501        if(numrifts)riftflag=1;
     502
     503        /*Assign output pointers:*/
     504        *priftflag=riftflag;
     505        *pnumrifts=numrifts;
     506        return noerr;
     507}/*}}}*/
     508int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){/*{{{*/
     509
     510        int noerr=1;
     511        int i,j,k,counter;
     512
     513        /*intermediary: */
     514        int *riftsegments = NULL;
     515        int *riftpairs    = NULL;
     516        int numsegs;
     517
     518        /*ordering and copy: */
     519        int *order             = NULL;
     520        int *riftsegments_copy = NULL;
     521        int *riftpairs_copy    = NULL;
     522
     523        /*node and element manipulation: */
     524        int node1,node2,node3,node4,temp_node,tip1,tip2,node;
     525        int el2;
     526        int already_ordered=0;
     527
     528        /*output: */
     529        int* riftstips=NULL;
     530
     531        /*Allocate byproduct of this routine, riftstips: */
     532        riftstips=xNew<int>(numrifts*2);
     533
     534        /*Go through all rifts: */
     535        for (i=0;i<numrifts;i++){
     536                riftsegments = riftssegments[i];
     537                riftpairs    = riftspairs[i];
     538                numsegs      = riftsnumsegments[i];
     539
     540                /*Allocate copy of riftsegments and riftpairs,
     541                 *as well as ordering vector: */
     542                riftsegments_copy=xNew<int>(numsegs*3);
     543                riftpairs_copy=xNew<int>(numsegs*2);
     544                order=xNew<int>(numsegs);
     545
     546                /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
     547                tip1=-1;
     548                tip2=-1;
     549
     550                for (j=0;j<numsegs;j++){
     551                        el2=*(riftpairs+2*j+1);
     552                        node1=*(riftsegments+3*j+0);
     553                        node2=*(riftsegments+3*j+1);
     554                        /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
     555                         *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
     556                        for (k=0;k<numsegs;k++){
     557                                if (*(riftsegments+3*k+2)==el2){
     558                                        node3=*(riftsegments+3*k+0);
     559                                        node4=*(riftsegments+3*k+1);
     560                                        break;
     561                                }
     562                        }
     563                        /* Make sure node3 faces node1 and node4 faces node2: */
     564                        _assert_(node1<nods+1 && node4<nods+1);
     565                        _assert_(node1>0 && node4>0);
     566                        if ((x[node1-1]==x[node4-1]) && (y[node1-1]==y[node4-1])){
     567                                /*Swap node3 and node4:*/
     568                                temp_node=node3;
     569                                node3=node4;
     570                                node4=temp_node;
     571                        }
     572
     573                        /*Figure out if a tip is on this element: */
     574                        if (node3==node1){
     575                                /*node1 is a tip*/
     576                                if (tip1==-1) {
     577                                        tip1=node1;
     578                                        continue;
     579                                }
     580                                if ((tip2==-1) && (node1!=tip1)){
     581                                        tip2=node1;
     582                                        break;
     583                                }
     584                        }
     585
     586                        if (node4==node2){
     587                                /*node2 is a tip*/
     588                                if (tip1==-1){
     589                                        tip1=node2;
     590                                        continue;
     591                                }
     592                                if ((tip2==-1) && (node2!=tip1)){
     593                                        tip2=node2;
     594                                        break;
     595                                }
     596                        }
     597                }
     598
     599                /*Record tips in riftstips: */
     600                *(riftstips+2*i+0)=tip1;
     601                *(riftstips+2*i+1)=tip2;
     602
     603                /*We have the two tips for this rift.  Go from tip1 to tip2, and figure out the order in which segments are sequential.
     604                 *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
     605                node=tip1;
     606                for (counter=0;counter<numsegs;counter++){
     607                        for (j=0;j<numsegs;j++){
     608                                node1=*(riftsegments+3*j+0);
     609                                node2=*(riftsegments+3*j+1);
     610
     611                                if ((node1==node) || (node2==node)){
     612                                        /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
     613                                        already_ordered=0;
     614                                        for (k=0;k<counter;k++){
     615                                                if(order[k]==j){
     616                                                        already_ordered=1;
     617                                                        break;
     618                                                }
     619                                        }
     620                                        if (!already_ordered){
     621                                                order[counter]=j;
     622                                                if(node1==node){
     623                                                        node=node2;
     624                                                }
     625                                                else if(node2==node){
     626                                                        node=node1;
     627                                                }
     628                                                break;
     629                                        }
     630                                }
     631                        }
     632                }
     633
     634                /*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
     635                for (j=0;j<numsegs;j++){
     636                        _assert_(order[j]<numsegs);
     637                        *(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
     638                        *(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
     639                        *(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
     640                        *(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
     641                        *(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
     642                }
     643
     644                for (j=0;j<numsegs;j++){
     645                        *(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
     646                        *(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
     647                        *(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
     648                        *(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
     649                        *(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
     650                }
     651
     652                xDelete<int>(order);
     653                xDelete<int>(riftsegments_copy);
     654                xDelete<int>(riftpairs_copy);
     655
     656        }
     657
     658        /*Assign output pointer:*/
     659        *priftstips=riftstips;
     660        return noerr;
     661}/*}}}*/
     662int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,/*{{{*/
     663                int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
     664
     665        int noerr=1;
     666        int i,j,k,k0;
     667
     668        double el1,el2,node1,node2,node3,node4;
     669        double temp_node;
     670
     671        /*output: */
     672        double **riftspenaltypairs    = NULL;
     673        double  *riftpenaltypairs     = NULL;
     674        int     *riftsnumpenaltypairs = NULL;
     675
     676        /*intermediary: */
     677        int numsegs;
     678        int* riftsegments=NULL;
     679        int* riftpairs=NULL;
     680        int counter;
     681        double normal[2];
     682        double length;
     683        int    k1,k2;
     684
     685        /*Allocate: */
     686        riftspenaltypairs=xNew<double*>(numrifts);
     687        riftsnumpenaltypairs=xNew<int>(numrifts);
     688
     689        for(i=0;i<numrifts;i++){
     690                numsegs=riftsnumsegs[i];
     691                riftsegments=riftssegments[i];
     692                riftpairs=riftspairs[i];
     693
     694                /*allocate riftpenaltypairs, and riftnumpenaltypairs: */
     695                if((numsegs/2-1)!=0)riftpenaltypairs=xNewZeroInit<double>((numsegs/2-1)*RIFTPENALTYPAIRSWIDTH);
     696
     697                /*Go through only one flank of the rifts, not counting the tips: */
     698                counter=0;
     699                for(j=0;j<(numsegs/2);j++){
     700                        el1=*(riftpairs+2*j+0);
     701                        el2=*(riftpairs+2*j+1);
     702                        node1=*(riftsegments+3*j+0);
     703                        node2=*(riftsegments+3*j+1);
     704                        /*Find segment index to recover node3 and node4, facing node1 and node2: */
     705                        k0=-1;
     706                        for(k=0;k<numsegs;k++){
     707                                if(*(riftsegments+3*k+2)==el2){
     708                                        k0=k;
     709                                        break;
     710                                }
     711                        }
     712                        node3=*(riftsegments+3*k0+0);
     713                        node4=*(riftsegments+3*k0+1);
     714
     715                        /* Make sure node3 faces node1 and node4 faces node2: */
     716                        if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
     717                                /*Swap node3 and node4:*/
     718                                temp_node=node3;
     719                                node3=node4;
     720                                node4=temp_node;
     721                        }       
     722                        /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
     723                         *this segment, and its length: */
     724                        normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
     725                        normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
     726                        length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
     727
     728                        /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
     729                         * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4,
     730                         * only once. We'll add the normals and the lengths : */
     731
     732                        if(node1!=node3){ //exclude tips from loads
     733                                k1=-1;
     734                                for(k=0;k<counter;k++){
     735                                        if( (*(riftpenaltypairs+k*7+0))==node1){
     736                                                k1=k;
     737                                                break;
     738                                        }
     739                                }
     740                                if(k1==-1){
     741                                        *(riftpenaltypairs+counter*7+0)=node1;
     742                                        *(riftpenaltypairs+counter*7+1)=node3;
     743                                        *(riftpenaltypairs+counter*7+2)=el1;
     744                                        *(riftpenaltypairs+counter*7+3)=el2;
     745                                        *(riftpenaltypairs+counter*7+4)=normal[0];
     746                                        *(riftpenaltypairs+counter*7+5)=normal[1];
     747                                        *(riftpenaltypairs+counter*7+6)=length/2;
     748                                        counter++;
     749                                }
     750                                else{
     751                                        *(riftpenaltypairs+k1*7+4)+=normal[0];
     752                                        *(riftpenaltypairs+k1*7+5)+=normal[1];
     753                                        *(riftpenaltypairs+k1*7+6)+=length/2;
     754                                }
     755                        }
     756                        if(node2!=node4){
     757                                k2=-1;
     758                                for(k=0;k<counter;k++){
     759                                        if( (*(riftpenaltypairs+k*7+0))==node2){
     760                                                k2=k;
     761                                                break;
     762                                        }
     763                                }
     764                                if(k2==-1){
     765                                        *(riftpenaltypairs+counter*7+0)=node2;
     766                                        *(riftpenaltypairs+counter*7+1)=node4;
     767                                        *(riftpenaltypairs+counter*7+2)=el1;
     768                                        *(riftpenaltypairs+counter*7+3)=el2;
     769                                        *(riftpenaltypairs+counter*7+4)=normal[0];
     770                                        *(riftpenaltypairs+counter*7+5)=normal[1];
     771                                        *(riftpenaltypairs+counter*7+6)=length/2;
     772                                        counter++;
     773                                }
     774                                else{
     775                                        *(riftpenaltypairs+k2*7+4)+=normal[0];
     776                                        *(riftpenaltypairs+k2*7+5)+=normal[1];
     777                                        *(riftpenaltypairs+k2*7+6)+=length/2;
     778                                }
     779                        }
     780                }
     781                /*Renormalize normals: */
     782                for(j=0;j<counter;j++){
     783                        double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
     784                        *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
     785                        *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
     786                }
     787
     788                riftspenaltypairs[i]=riftpenaltypairs;
     789                riftsnumpenaltypairs[i]=(numsegs/2-1);
     790        }
     791
     792        /*Assign output pointers: */
     793        *priftspenaltypairs=riftspenaltypairs;
     794        *priftsnumpenaltypairs=riftsnumpenaltypairs;
     795        return noerr;
     796}/*}}}*/
     797int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
     798
     799        int noerr=1;
     800        int i,j,k;
     801        int node1,node2,node3;
     802        int el;
     803        double  pair[2];
     804        int     pair_count=0;
     805        int     triple=0;
     806
     807        /*Recover input: */
     808        int    *index = *pindex;
     809        int     nel   = *pnel;
     810        double *x     = *px;
     811        double *y     = *py;
     812        int     nods  = *pnods;
     813
     814        for (i=0;i<num_seg;i++){
     815                node1=*(segments+3*i+0);
     816                node2=*(segments+3*i+1);
     817                /*Find all elements connected to [node1 node2]: */
     818                pair_count=0;
     819                for (j=0;j<nel;j++){
     820                        if (*(index+3*j+0)==node1){
     821                                if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
     822                                        pair[pair_count]=j;
     823                                        pair_count++;
     824                                }
     825                        }
     826                        if (*(index+3*j+1)==node1){
     827                                if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
     828                                        pair[pair_count]=j;
     829                                        pair_count++;
     830                                }
     831                        }
     832                        if (*(index+3*j+2)==node1){
     833                                if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
     834                                        pair[pair_count]=j;
     835                                        pair_count++;
     836                                }
     837                        }
     838                }
     839                /*Ok, we have pair_count elements connected to this segment. For each of these elements,
     840                 *figure out if the third node also belongs to a segment: */
     841                if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to  2 elements
     842                        continue;
     843                }
     844                else{
     845                        for (j=0;j<pair_count;j++){
     846                                el=(int)pair[j];
     847                                triple=0;
     848                                /*First find node3: */
     849                                if (*(index+3*el+0)==node1){
     850                                        if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
     851                                        else node3=*(index+3*el+1);
     852                                }
     853                                if (*(index+3*el+1)==node1){
     854                                        if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
     855                                        else node3=*(index+3*el+0);
     856                                }
     857                                if (*(index+3*el+2)==node1){
     858                                        if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
     859                                        else node3=*(index+3*el+0);
     860                                }
     861                                /*Ok, we have node3. Does node3 belong to a segment? : */
     862                                for (k=0;k<num_seg;k++){
     863                                        if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
     864                                                triple=1;
     865                                                break;
     866                                        }
     867                                }
     868                                if(triple==1){
     869                                        /*el is a corner element: we need to split it in 3 triangles: */
     870                                        x=xReNew<double>(x,nods,nods+1);
     871                                        y=xReNew<double>(y,nods,nods+1);
     872                                        x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
     873                                        y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
     874                                        index=xReNew<int>(index,nel*3,(nel+2*3));
     875                                        /*First, reassign element el: */
     876                                        *(index+3*el+0)=node1;
     877                                        *(index+3*el+1)=node2;
     878                                        *(index+3*el+2)=nods+1;
     879                                        /*Other two elements: */
     880                                        *(index+3*nel+0)=node2;
     881                                        *(index+3*nel+1)=node3;
     882                                        *(index+3*nel+2)=nods+1;
     883
     884                                        *(index+3*(nel+1)+0)=node3;
     885                                        *(index+3*(nel+1)+1)=node1;
     886                                        *(index+3*(nel+1)+2)=nods+1;
     887                                        /*we need  to change the segment elements corresponding to el: */
     888                                        for (k=0;k<num_seg;k++){
     889                                                if (*(segments+3*k+2)==(el+1)){
     890                                                        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;
     891                                                        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;
     892                                                        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;
     893                                                }
     894                                        }
     895
     896                                        nods=nods+1;
     897                                        nel=nel+2;
     898                                        i=0;
     899                                        break;
     900                                }
     901                        } //for (j=0;j<pair_count;j++)
     902                }
     903        }// for (i=0;i<num_seg;i++)
     904
     905        /*Assign output pointers: */
     906        *pindex=index;
     907        *pnel=nel;
     908        *px=x;
     909        *py=y;
     910        *pnods=nods;
     911        return noerr;
     912}/*}}}*/

  • ../trunk-jpl/src/c/shared/Triangle/SplitMeshForRifts.cpp

     
     1/*
     2 * SplitMeshForRifts.c:
     3 */
     4#include "./trimesh.h"
     5#include "../MemOps/MemOps.h"
     6
     7int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist){
     8
     9        /*Some notes on dimensions:
     10        index  of size nelx3
     11        x and y of size nodsx1
     12        segments of size nsegsx3*/
     13
     14        int i,j,k,l;
     15        int node;
     16        int el;
     17        int  nriftsegs;
     18        int* riftsegments=NULL;
     19        int* flags=NULL;
     20        int  NumGridElementListOnOneSideOfRift;
     21        int* GridElementListOnOneSideOfRift=NULL;
     22
     23        /*Recover input: */
     24        int     nel               = *pnel;
     25        int    *index             = *pindex;
     26        int     nods              = *pnods;
     27        double *x                 = *px;
     28        double *y                 = *py;
     29        int     nsegs             = *pnsegs;
     30        int    *segments          = *psegments;
     31        int    *segmentmarkerlist = *psegmentmarkerlist;
     32
     33        /*Establish list of segments that belong to a rift: */
     34        /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
     35        RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);
     36
     37        /*Go through all nodes of the rift segments, and start splitting the mesh: */
     38        flags=xNewZeroInit<int>(nods); //to make sure we don't split the same nodes twice!
     39        for (i=0;i<nriftsegs;i++){
     40                for (j=0;j<2;j++){
     41
     42                        node=riftsegments[4*i+j+2];
     43                        if(flags[node-1]){
     44                                /*This node was already split, skip:*/
     45                                continue;
     46                        }
     47                        else{
     48                                flags[node-1]=1;
     49                        }
     50
     51                        if(IsGridOnRift(riftsegments,nriftsegs,node)){
     52
     53                                DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
     54
     55                                /*Summary: we have for node, a list of elements
     56                                 * (GridElementListOnOneSideOfRift, of size
     57                                 * NumGridElementListOnOneSideOfRift) that all contain node
     58                                 *and that are on the same side of the rift. For all these
     59                                 elements, we clone node into another node, and we swap all
     60                                 instances of node in the triangulation *for those elements, to the
     61                                 new node.*/
     62
     63                                //create new node
     64                                x=xReNew<double>(x,nods,nods+1);
     65                                y=xReNew<double>(y,nods,nods+1);
     66                                x[nods]=x[node-1]; //matlab indexing
     67                                y[nods]=y[node-1]; //matlab indexing
     68
     69                                //augment number of nodes
     70                                nods++;
     71
     72                                //change elements owning this node
     73                                for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
     74                                        el=GridElementListOnOneSideOfRift[k];
     75                                        for (l=0;l<3;l++){
     76                                                if (index[3*el+l]==node) index[3*el+l]=nods; //again, matlab indexing.
     77                                        }
     78                                }
     79                        }// if(IsGridOnRift(riftsegments,nriftsegs,node))
     80                } //for(j=0;j<2;j++)
     81        } //for (i=0;i<nriftsegs;i++)
     82
     83        /*update segments: they got modified completely by adding new nodes.*/
     84        UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
     85
     86        /*Assign output pointers: */
     87        *pnel=nel;
     88        *pindex=index;
     89        *pnods=nods;
     90        *px=x;
     91        *py=y;
     92        *pnsegs=nsegs;
     93        *psegments=segments;
     94        *psegmentmarkerlist=segmentmarkerlist;
     95        return 1;
     96}

  • ../trunk-jpl/src/c/shared/Triangle/GridInsideHole.cpp

     
     1/*
     2 * GridInsideHole.c:
     3 * from a convex set of points, figure out a point that for sure lies inside the profile.
     4 */
     5
     6#include <math.h>
     7
     8#include "./trimesh.h"
     9#include "../Exp/exp.h"
     10
     11#undef M_PI
     12#define M_PI 3.141592653589793238462643
     13
     14int GridInsideHole(double* px0,double* py0,int n,double* x,double* y){
     15
     16        double flag=0.0;
     17        double xA,xB,xC,xD,xE;
     18        double yA,yB,yC,yD,yE;
     19
     20        /*Take first and last vertices: */
     21        xA=x[0];
     22        yA=y[0];
     23        xB=x[n-1];
     24        yB=y[n-1];
     25
     26        /*Figure out middle of segment [A B]: */
     27        xC=(xA+xB)/2;
     28        yC=(yA+yB)/2;
     29
     30        /*D and E are on each side of segment [A B], on the median line between segment [A  B],
     31         *at an angle of 10 degree (less than the minimum 30 enforced by the quality of the mesh: */
     32        xD=xC+tan(10./180.*M_PI)*(yC-yA);
     33        yD=yC+tan(10./180.*M_PI)*(xA-xC);
     34        xE=xC-tan(10./180.*M_PI)*(yC-yA);
     35        yE=yC-tan(10./180.*M_PI)*(xA-xC);
     36
     37        /*Either E or D is inside profile (x,y): */
     38        IsInPolySerial(&flag,&xD,&yD,1,x,y,n,2);
     39        /*FIXME: used to be 'flag' and not '!flag', check*/
     40        if(!flag){
     41                /*D is inside the poly: */
     42                *px0=xD;
     43                *py0=yD;
     44        }
     45        else{
     46                /*E is inside the poly: */
     47                *px0=xE;
     48                *py0=yE;
     49        }
     50        return 1;
     51}

  • ../trunk-jpl/src/c/shared/Triangle/triangle.h

     
     1/*!\file:  triangle.h
     2 * \brief
     3 */
     4
     5#ifndef _SHARED_TRIANGLE_H
     6#define _SHARED_TRIANGLE_H
     7
     8#include <stdio.h>
     9#include <math.h>
     10
     11//#define REAL double //took  it out because it may conflict with stdlib.h defines. put back if necessary
     12int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel);
     13int OrderSegments(int** psegments,int nseg, int* index,int nel);
     14int GridInsideHole(double* px0,double* py0,int n,double* x,double* y);
     15int FindElement(int A,int B,int* index,int nel);
     16int SplitMeshForRifts(int* pnel,int** pindex,int* pnods,double** px,double** py,int* pnsegs,int** psegments,int** psegmentmarkerlist);
     17int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
     18int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
     19int IsNeighbor(int el1,int el2,int* index);
     20int IsOnRift(int el,int nriftsegs,int* riftsegments);
     21void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index, int nsegs,int* segments);
     22int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift,int segmentnumber, int nriftsegs,int* riftsegments, int node,int* index,int nel);
     23int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
     24int FindElement(double A,double B,int* index,int nel);
     25int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
     26int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels);
     27int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels);
     28int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int**  riftssegments,
     29                int* riftsnumsegments,int** riftspairs,int* riftstips,double* x,double* y);
     30int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg);
     31int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y);
     32
     33#endif  /* _SHARED_TRIANGLE_H */

  • ../trunk-jpl/src/c/shared/Triangle/AssociateSegmentToElement.cpp

     
     1/*!\file:  AssociateSegmentToElement.cpp
     2 * \brief for each segment, look for the corresponding element.
     3 */
     4
     5#include "./trimesh.h"
     6
     7int AssociateSegmentToElement(int** psegments,int nseg,int* index,int nel){
     8
     9        /*node indices: */
     10        int A,B;
     11
     12        /*Recover segments: */
     13        int* segments=*psegments;
     14
     15        for(int i=0;i<nseg;i++){
     16                A=segments[3*i+0];
     17                B=segments[3*i+1];
     18                segments[3*i+2]=FindElement(A,B,index,nel)+1; //matlab indexing.
     19        }
     20
     21        /*Assign output pointers: */
     22        *psegments=segments;
     23        return 1;
     24}

  • ../trunk-jpl/src/c/shared/Triangle

  • ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.h 

    Property changes on: ../trunk-jpl/src/c/shared/Triangle
___________________________________________________________________
Added: svn:ignore
## -0,0 +1,2 ##
+.deps
+.dirstamp
     
    1 /*
    2  * TriMeshProcessRifts.h
    3  */
    4 
    5 #ifndef _TRIMESH_PROCESSRIFTS_H_
    6 #define _TRIMESH_PROCESSRIFTS_H_
    7 
    8 #ifdef HAVE_CONFIG_H
    9 #include <config.h>
    10 #else
    11 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    12 #endif
    13 
    14 /*For python modules: needs to come before header files inclusion*/
    15 #ifdef _HAVE_PYTHON_
    16 #define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
    17 #endif
    18 
    19 #include "../bindings.h"
    20 #include "../../c/main/globals.h"
    21 #include "../../c/modules/modules.h"
    22 #include "../../c/shared/shared.h"
    23 
    24 #undef __FUNCT__
    25 #define __FUNCT__  "TriMeshProcessRifts"
    26 
    27 #ifdef _HAVE_MATLAB_MODULES_
    28 /* serial input macros: */
    29 #define INDEXIN          prhs[0]
    30 #define XIN              prhs[1]
    31 #define YIN              prhs[2]
    32 #define SEGMENTSIN       prhs[3]
    33 #define SEGMENTMARKERSIN prhs[4]
    34 /* serial output macros: */
    35 #define INDEXOUT          (mxArray**)&plhs[0]
    36 #define XOUT              (mxArray**)&plhs[1]
    37 #define YOUT              (mxArray**)&plhs[2]
    38 #define SEGMENTSOUT       (mxArray**)&plhs[3]
    39 #define SEGMENTMARKERSOUT (mxArray**)&plhs[4]
    40 #define RIFTSTRUCT        (mxArray**)&plhs[5]
    41 #endif
    42 
    43 #ifdef _HAVE_PYTHON_MODULES_
    44 /* serial input macros: */
    45 #define INDEXIN          PyTuple_GetItem(args,0)
    46 #define XIN              PyTuple_GetItem(args,1)
    47 #define YIN              PyTuple_GetItem(args,2)
    48 #define SEGMENTSIN       PyTuple_GetItem(args,3)
    49 #define SEGMENTMARKERSIN PyTuple_GetItem(args,4)
    50 /* serial output macros: */
    51 #define INDEXOUT          output,0
    52 #define XOUT              output,1
    53 #define YOUT              output,2
    54 #define SEGMENTSOUT       output,3
    55 #define SEGMENTMARKERSOUT output,4
    56 #define RIFTSTRUCT        output,5
    57 #endif
    58 
    59 /* serial arg counts: */
    60 #undef NLHS
    61 #define NLHS  6
    62 #undef NRHS
    63 #define NRHS  5
    64 
    65 #endif

  • ../trunk-jpl/src/wrappers/TriMeshProcessRifts/TriMeshProcessRifts.cpp

     
    1 /*!\file:  TriMeshProcessRifts.cpp
    2  * \brief split a mesh where a rift (or fault) is present
    3  */
    4 
    5 #include "./TriMeshProcessRifts.h"
    6 
    7 void TriMeshProcessRiftsUsage(void){/*{{{*/
    8         _printf_("\n");
    9         _printf_("   usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(index1,x1,y1,segments1,segmentmarkers1) \n");
    10         _printf_("      where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
    11         _printf_("      index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
    12 }/*}}}*/
    13 WRAPPER(TriMeshProcessRifts_python){
    14 
    15         /* returned quantities: */
    16         RiftStruct *riftstruct = NULL;
    17 
    18         /* input: */
    19         int     nel,nods;
    20         int    *index          = NULL;
    21         double *x              = NULL;
    22         double *y              = NULL;
    23         int    *segments       = NULL;
    24         int    *segmentmarkers = NULL;
    25         int     num_seg;
    26 
    27         /*Boot module*/
    28         MODULEBOOT();
    29 
    30         /*checks on arguments on the matlab side: */
    31         CHECKARGUMENTS(NLHS,NRHS,&TriMeshProcessRiftsUsage);
    32 
    33         /*Fetch data */
    34         FetchData(&index,&nel,NULL,INDEXIN);
    35         FetchData(&x,&nods,XIN);
    36         FetchData(&y,NULL,YIN);
    37         FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
    38         FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
    39 
    40         /*call x layer*/
    41         TriMeshProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,&riftstruct);
    42 
    43         /*Output : */
    44         WriteData(INDEXOUT,index,nel,3);
    45         WriteData(XOUT,x,nods,1);
    46         WriteData(YOUT,y,nods,1);
    47         WriteData(SEGMENTSOUT,segments,num_seg,3);
    48         WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
    49         WriteData(RIFTSTRUCT,riftstruct);
    50 
    51         /*end module: */
    52         delete riftstruct;
    53         xDelete<int>(index);
    54         xDelete<double>(x);
    55         xDelete<double>(y);
    56         xDelete<int>(segments);
    57         xDelete<int>(segmentmarkers );
    58         MODULEEND();
    59 }

  • ../trunk-jpl/src/wrappers/TriMeshProcessRifts

  • ../trunk-jpl/src/wrappers/TriMesh/TriMesh.cpp 

    Property changes on: ../trunk-jpl/src/wrappers/TriMeshProcessRifts
___________________________________________________________________
Deleted: svn:ignore
## -1,2 +0,0 ##
-.deps
-.dirstamp
     
    1 /*
    2  * TriMesh: mesh a domain using an .exp file
    3  */
    4 
    5 #include "./TriMesh.h"
    6 
    7 void TriMeshUsage(void){/*{{{*/
    8         _printf_("\n");
    9         _printf_("   usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,area) \n");
    10         _printf_("      where: index,x,y defines a triangulation, segments is an array made \n");
    11         _printf_("      of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment, \n");
    12         _printf_("      outlinefilename an Argus domain outline file, \n");
    13         _printf_("      area is the maximum area desired for any element of the resulting mesh, \n");
    14         _printf_("\n");
    15 }/*}}}*/
    16 WRAPPER(TriMesh_python){
    17        
    18         /*intermediary: */
    19         double    area;
    20         Contours *domain = NULL;
    21         Contours *rifts  = NULL;
    22 
    23         /* output: */
    24         int    *index             = NULL;
    25         double *x                 = NULL;
    26         double *y                 = NULL;
    27         int    *segments          = NULL;
    28         int    *segmentmarkerlist = NULL;
    29         int     nel,nods,nsegs;
    30 
    31         /*Boot module: */
    32         MODULEBOOT();
    33 
    34         /*checks on arguments: */
    35         CHECKARGUMENTS(NLHS,NRHS,&TriMeshUsage);
    36 
    37         /*Fetch data needed for meshing: */
    38         FetchData(&domain,DOMAINOUTLINE);
    39         FetchData(&rifts,RIFTSOUTLINE);
    40         FetchData(&area,AREA);
    41 
    42         /*call x core: */
    43         TriMeshx(&index,&x,&y,&segments,&segmentmarkerlist,&nel,&nods,&nsegs,domain,rifts,area);
    44 
    45         /*write outputs: */
    46         WriteData(INDEX,index,nel,3);
    47         WriteData(X,x,nods);
    48         WriteData(Y,y,nods);
    49         WriteData(SEGMENTS,segments,nsegs,3);
    50         WriteData(SEGMENTMARKERLIST,segmentmarkerlist,nsegs);
    51 
    52         /*free ressources: */
    53         delete domain;
    54         delete rifts;
    55         xDelete<int>(index);
    56         xDelete<double>(x);
    57         xDelete<double>(y);
    58         xDelete<int>(segments);
    59         xDelete<int>(segmentmarkerlist);
    60 
    61         /*end module: */
    62         MODULEEND();
    63 }

  • ../trunk-jpl/src/wrappers/TriMesh/TriMesh.h

     
    1 /*
    2         TriMesh.h
    3 */
    4 
    5 #ifndef _TRIMESH_H
    6 #define _TRIMESH_H
    7 
    8 #ifdef HAVE_CONFIG_H
    9         #include <config.h>
    10 #else
    11         #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    12 #endif
    13 
    14 /*For python modules: needs to come before header files inclusion*/
    15 #ifdef _HAVE_PYTHON_
    16 #define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
    17 #endif
    18 
    19 #ifdef _HAVE_JAVASCRIPT_MODULES_
    20 #undef _DO_NOT_LOAD_GLOBALS_ /*only module where this needs to be undefined, so as to
    21                                                            not include IssmComm several times in the javascript Modle construct.*/
    22 #endif
    23 
    24 /*Header files: */
    25 #include "../bindings.h"
    26 #include "../../c/main/globals.h"
    27 #include "../../c/toolkits/toolkits.h"
    28 #include "../../c/modules/modules.h"
    29 #include "../../c/shared/shared.h"
    30 #include "../../c/shared/io/io.h"
    31 
    32 #undef __FUNCT__
    33 #define __FUNCT__  "TriMesh"
    34 
    35 #ifdef _HAVE_MATLAB_MODULES_
    36 /* serial input macros: */
    37 #define DOMAINOUTLINE  prhs[0]
    38 #define RIFTSOUTLINE   prhs[1]
    39 #define AREA           prhs[2]
    40 /* serial output macros: */
    41 #define INDEX             (mxArray**)&plhs[0]
    42 #define X                 (mxArray**)&plhs[1]
    43 #define Y                 (mxArray**)&plhs[2]
    44 #define SEGMENTS          (mxArray**)&plhs[3]
    45 #define SEGMENTMARKERLIST (mxArray**)&plhs[4]
    46 #endif
    47 
    48 #ifdef _HAVE_PYTHON_MODULES_
    49 /* serial input macros: */
    50 #define DOMAINOUTLINE PyTuple_GetItem(args,0)
    51 #define RIFTSOUTLINE  PyTuple_GetItem(args,1)
    52 #define AREA          PyTuple_GetItem(args,2)
    53 /* serial output macros: */
    54 #define INDEX             output,0
    55 #define X                 output,1
    56 #define Y                 output,2
    57 #define SEGMENTS          output,3
    58 #define SEGMENTMARKERLIST output,4
    59 #endif
    60 
    61 #ifdef _HAVE_JAVASCRIPT_MODULES_
    62 /* serial input macros: */
    63 #define DOMAINOUTLINE domainx,domainy,domainnods
    64 #define RIFTSOUTLINE  NULL,NULL,0
    65 #define AREA          areain
    66 /* serial output macros: */
    67 #define INDEX             pindex,pnel
    68 #define X                 px,pnods
    69 #define Y                 py,pnods
    70 #define SEGMENTS          psegments,pnsegs
    71 #define SEGMENTMARKERLIST psegmentmarkers,pnsegs
    72 #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)
    73 #define _DO_NOT_LOAD_GLOBALS_//we only load globals for TriMeshModule.js, not other modules!
    74 #endif
    75 
    76 
    77 /* serial arg counts: */
    78 #undef NLHS
    79 #define NLHS  5
    80 #undef NRHS
    81 #define NRHS  3
    82 
    83 #endif  /* _TRIMESH_H */

  • ../trunk-jpl/src/wrappers/TriMesh/TriMesh.js

     
    1 function TriMesh(md,domain,rifts, area){
    2 /*TriMesh
    3            usage: var array = TriMesh(domain,rifts,area);
    4               where: array is made of [index,x,y,segments,segmentmarkers]
    5                   and index,x,y defines a triangulation, segments is an array made
    6               of exterior segments to the mesh domain outline, segmentmarkers is an array
    7                   flagging each segment, domain a js array defining the domain outline  (sames for
    8                   rifts) and area is the maximum area desired for any element of the resulting mesh.
    9 
    10                   Ok, for now, we are not dealing with rifts. Also, the domain is made of only one
    11                   profile, this to avoid passing a double** pointer to js.
    12 */
    13 
    14         //Dynamic allocations: {{{
    15         //Retrieve domain arrays, and allocate on Module heap:
    16        
    17         //input
    18         var dx=new Float64Array(domain['x']); var nx=dx.length * dx.BYTES_PER_ELEMENT;
    19         var dxPtr= Module._malloc(nx); var domainxHeap = new Uint8Array(Module.HEAPU8.buffer,dxPtr,nx);
    20         domainxHeap.set(new Uint8Array(dx.buffer)); var domainx=domainxHeap.byteOffset;
    21 
    22         var dy=new Float64Array(domain['y']); var ny=dy.length * dy.BYTES_PER_ELEMENT;
    23         var dyPtr = Module._malloc(ny); var domainyHeap = new Uint8Array(Module.HEAPU8.buffer,dyPtr,ny);
    24         domainyHeap.set(new Uint8Array(dy.buffer)); var domainy=domainyHeap.byteOffset;
    25        
    26         //output
    27         var nel,indexlinear,index,nods,x,y;
    28         var pnel= Module._malloc(4);
    29         var pindex= Module._malloc(4);
    30         var pnods= Module._malloc(4);
    31         var px= Module._malloc(4);
    32         var py= Module._malloc(4);
    33         var psegments= Module._malloc(4);
    34         var psegmentmarkers= Module._malloc(4);
    35         var pnsegs= Module._malloc(4);
    36         //}}}
    37 
    38         //Declare TriMesh module:
    39         TriMeshModule = Module.cwrap('TriMeshModule','number',['number','number','number','number','number','number','number','number','number','number','number','number']);
    40        
    41         //Call TriMesh module:
    42         TriMeshModule(pindex,px,py,pnel,pnods,psegments,psegmentmarkers,pnsegs, domainx,domainy,dx.length,area);
    43        
    44         /*Dynamic copying from heap: {{{*/
    45         //recover mesh:
    46         nel = Module.getValue(pnel, 'i32');
    47         var indexptr = Module.getValue(pindex,'i32');
    48         indexlinear = Module.HEAPF64.slice(indexptr /8, indexptr/8 + nel*3);
    49         index = ListToMatrix(indexlinear,3);
    50 
    51         nods = Module.getValue(pnods, 'i32');
    52         var xptr = Module.getValue(px,'i32');
    53         var yptr = Module.getValue(py,'i32');
    54         x = Module.HEAPF64.slice(xptr /8, xptr/8 + nods);
    55         y = Module.HEAPF64.slice(yptr /8, yptr/8 + nods);
    56        
    57         nsegs = Module.getValue(pnsegs, 'i32');
    58         var segmentsptr = Module.getValue(psegments,'i32');
    59         segmentslinear = Module.HEAPF64.slice(segmentsptr /8, segmentsptr/8 + nsegs*3);
    60         segments = ListToMatrix(segmentslinear,3);
    61        
    62         var segmentmarkersptr = Module.getValue(psegmentmarkers,'i32');
    63         segmentmarkers = Module.HEAPF64.slice(segmentmarkersptr /8, segmentmarkersptr/8 + nsegs);
    64         /*}}}*/
    65 
    66         var return_array=[index,x,y,segments,segmentmarkers];
    67 
    68         /*Free ressources: */
    69         Module._free(pindex);
    70         Module._free(indexlinear);
    71         Module._free(px);
    72         Module._free(x);
    73         Module._free(py);
    74         Module._free(y);
    75         Module._free(pnel);
    76         Module._free(pnods);
    77         Module._free(psegments);
    78         Module._free(psegmentmarkers);
    79         Module._free(pnsegs);
    80 
    81         return return_array;
    82 }

  • ../trunk-jpl/src/wrappers/TriMesh

  • ../trunk-jpl/src/wrappers/Triangle/Triangle.js 

    Property changes on: ../trunk-jpl/src/wrappers/TriMesh
___________________________________________________________________
Deleted: svn:ignore
## -1,2 +0,0 ##
-.deps
-.dirstamp
     
     1function Triangle(md,domain,rifts, area){
     2/*Triangle
     3           usage: var array = Triangle(domain,rifts,area);
     4              where: array is made of [index,x,y,segments,segmentmarkers]
     5                  and index,x,y defines a triangulation, segments is an array made
     6              of exterior segments to the mesh domain outline, segmentmarkers is an array
     7                  flagging each segment, domain a js array defining the domain outline  (sames for
     8                  rifts) and area is the maximum area desired for any element of the resulting mesh.
     9
     10                  Ok, for now, we are not dealing with rifts. Also, the domain is made of only one
     11                  profile, this to avoid passing a double** pointer to js.
     12*/
     13
     14        //Dynamic allocations: {{{
     15        //Retrieve domain arrays, and allocate on Module heap:
     16       
     17        //input
     18        var dx=new Float64Array(domain['x']); var nx=dx.length * dx.BYTES_PER_ELEMENT;
     19        var dxPtr= Module._malloc(nx); var domainxHeap = new Uint8Array(Module.HEAPU8.buffer,dxPtr,nx);
     20        domainxHeap.set(new Uint8Array(dx.buffer)); var domainx=domainxHeap.byteOffset;
     21
     22        var dy=new Float64Array(domain['y']); var ny=dy.length * dy.BYTES_PER_ELEMENT;
     23        var dyPtr = Module._malloc(ny); var domainyHeap = new Uint8Array(Module.HEAPU8.buffer,dyPtr,ny);
     24        domainyHeap.set(new Uint8Array(dy.buffer)); var domainy=domainyHeap.byteOffset;
     25       
     26        //output
     27        var nel,indexlinear,index,nods,x,y;
     28        var pnel= Module._malloc(4);
     29        var pindex= Module._malloc(4);
     30        var pnods= Module._malloc(4);
     31        var px= Module._malloc(4);
     32        var py= Module._malloc(4);
     33        var psegments= Module._malloc(4);
     34        var psegmentmarkers= Module._malloc(4);
     35        var pnsegs= Module._malloc(4);
     36        //}}}
     37
     38        //Declare Triangle module:
     39        TriangleModule = Module.cwrap('TriangleModule','number',['number','number','number','number','number','number','number','number','number','number','number','number']);
     40       
     41        //Call Triangle module:
     42        TriangleModule(pindex,px,py,pnel,pnods,psegments,psegmentmarkers,pnsegs, domainx,domainy,dx.length,area);
     43       
     44        /*Dynamic copying from heap: {{{*/
     45        //recover mesh:
     46        nel = Module.getValue(pnel, 'i32');
     47        var indexptr = Module.getValue(pindex,'i32');
     48        indexlinear = Module.HEAPF64.slice(indexptr /8, indexptr/8 + nel*3);
     49        index = ListToMatrix(indexlinear,3);
     50
     51        nods = Module.getValue(pnods, 'i32');
     52        var xptr = Module.getValue(px,'i32');
     53        var yptr = Module.getValue(py,'i32');
     54        x = Module.HEAPF64.slice(xptr /8, xptr/8 + nods);
     55        y = Module.HEAPF64.slice(yptr /8, yptr/8 + nods);
     56       
     57        nsegs = Module.getValue(pnsegs, 'i32');
     58        var segmentsptr = Module.getValue(psegments,'i32');
     59        segmentslinear = Module.HEAPF64.slice(segmentsptr /8, segmentsptr/8 + nsegs*3);
     60        segments = ListToMatrix(segmentslinear,3);
     61       
     62        var segmentmarkersptr = Module.getValue(psegmentmarkers,'i32');
     63        segmentmarkers = Module.HEAPF64.slice(segmentmarkersptr /8, segmentmarkersptr/8 + nsegs);
     64        /*}}}*/
     65
     66        var return_array=[index,x,y,segments,segmentmarkers];
     67
     68        /*Free ressources: */
     69        Module._free(pindex);
     70        Module._free(indexlinear);
     71        Module._free(px);
     72        Module._free(x);
     73        Module._free(py);
     74        Module._free(y);
     75        Module._free(pnel);
     76        Module._free(pnods);
     77        Module._free(psegments);
     78        Module._free(psegmentmarkers);
     79        Module._free(pnsegs);
     80
     81        return return_array;
     82}

  • ../trunk-jpl/src/wrappers/Triangle/Triangle.cpp

     
     1/*
     2 * Triangle: mesh a domain using an .exp file
     3 */
     4
     5#include "./Triangle.h"
     6
     7void TriangleUsage(void){/*{{{*/
     8        _printf_("\n");
     9        _printf_("   usage: [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,area) \n");
     10        _printf_("      where: index,x,y defines a triangulation, segments is an array made \n");
     11        _printf_("      of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment, \n");
     12        _printf_("      outlinefilename an Argus domain outline file, \n");
     13        _printf_("      area is the maximum area desired for any element of the resulting mesh, \n");
     14        _printf_("\n");
     15}/*}}}*/
     16WRAPPER(Triangle_python){
     17       
     18        /*intermediary: */
     19        double    area;
     20        Contours *domain = NULL;
     21        Contours *rifts  = NULL;
     22
     23        /* output: */
     24        int    *index             = NULL;
     25        double *x                 = NULL;
     26        double *y                 = NULL;
     27        int    *segments          = NULL;
     28        int    *segmentmarkerlist = NULL;
     29        int     nel,nods,nsegs;
     30
     31        /*Boot module: */
     32        MODULEBOOT();
     33
     34        /*checks on arguments: */
     35        CHECKARGUMENTS(NLHS,NRHS,&TriangleUsage);
     36
     37        /*Fetch data needed for meshing: */
     38        FetchData(&domain,DOMAINOUTLINE);
     39        FetchData(&rifts,RIFTSOUTLINE);
     40        FetchData(&area,AREA);
     41
     42        /*call x core: */
     43        Trianglex(&index,&x,&y,&segments,&segmentmarkerlist,&nel,&nods,&nsegs,domain,rifts,area);
     44
     45        /*write outputs: */
     46        WriteData(INDEX,index,nel,3);
     47        WriteData(X,x,nods);
     48        WriteData(Y,y,nods);
     49        WriteData(SEGMENTS,segments,nsegs,3);
     50        WriteData(SEGMENTMARKERLIST,segmentmarkerlist,nsegs);
     51
     52        /*free ressources: */
     53        delete domain;
     54        delete rifts;
     55        xDelete<int>(index);
     56        xDelete<double>(x);
     57        xDelete<double>(y);
     58        xDelete<int>(segments);
     59        xDelete<int>(segmentmarkerlist);
     60
     61        /*end module: */
     62        MODULEEND();
     63}

  • ../trunk-jpl/src/wrappers/Triangle/Triangle.h

     
     1/*
     2        Triangle.h
     3*/
     4
     5#ifndef _TRIANGLE_H
     6#define _TRIANGLE_H
     7
     8#ifdef HAVE_CONFIG_H
     9        #include <config.h>
     10#else
     11        #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     12#endif
     13
     14/*For python modules: needs to come before header files inclusion*/
     15#ifdef _HAVE_PYTHON_
     16#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
     17#endif
     18
     19#ifdef _HAVE_JAVASCRIPT_MODULES_
     20#undef _DO_NOT_LOAD_GLOBALS_ /*only module where this needs to be undefined, so as to
     21                                                           not include IssmComm several times in the javascript Modle construct.*/
     22#endif
     23
     24/*Header files: */
     25#include "../bindings.h"
     26#include "../../c/main/globals.h"
     27#include "../../c/toolkits/toolkits.h"
     28#include "../../c/modules/modules.h"
     29#include "../../c/shared/shared.h"
     30#include "../../c/shared/io/io.h"
     31
     32#undef __FUNCT__
     33#define __FUNCT__  "Triangle"
     34
     35#ifdef _HAVE_MATLAB_MODULES_
     36/* serial input macros: */
     37#define DOMAINOUTLINE  prhs[0]
     38#define RIFTSOUTLINE   prhs[1]
     39#define AREA           prhs[2]
     40/* serial output macros: */
     41#define INDEX             (mxArray**)&plhs[0]
     42#define X                 (mxArray**)&plhs[1]
     43#define Y                 (mxArray**)&plhs[2]
     44#define SEGMENTS          (mxArray**)&plhs[3]
     45#define SEGMENTMARKERLIST (mxArray**)&plhs[4]
     46#endif
     47
     48#ifdef _HAVE_PYTHON_MODULES_
     49/* serial input macros: */
     50#define DOMAINOUTLINE PyTuple_GetItem(args,0)
     51#define RIFTSOUTLINE  PyTuple_GetItem(args,1)
     52#define AREA          PyTuple_GetItem(args,2)
     53/* serial output macros: */
     54#define INDEX             output,0
     55#define X                 output,1
     56#define Y                 output,2
     57#define SEGMENTS          output,3
     58#define SEGMENTMARKERLIST output,4
     59#endif
     60
     61#ifdef _HAVE_JAVASCRIPT_MODULES_
     62/* serial input macros: */
     63#define DOMAINOUTLINE domainx,domainy,domainnods
     64#define RIFTSOUTLINE  NULL,NULL,0
     65#define AREA          areain
     66/* serial output macros: */
     67#define INDEX             pindex,pnel
     68#define X                 px,pnods
     69#define Y                 py,pnods
     70#define SEGMENTS          psegments,pnsegs
     71#define SEGMENTMARKERLIST psegmentmarkers,pnsegs
     72#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)
     73#define _DO_NOT_LOAD_GLOBALS_//we only load globals for TriangleModule.js, not other modules!
     74#endif
     75
     76
     77/* serial arg counts: */
     78#undef NLHS
     79#define NLHS  5
     80#undef NRHS
     81#define NRHS  3
     82
     83#endif  /* _TRIANGLE_H */

  • ../trunk-jpl/src/wrappers/Triangle

  • ../trunk-jpl/src/wrappers/javascript/Makefile.am 

    Property changes on: ../trunk-jpl/src/wrappers/Triangle
___________________________________________________________________
Added: svn:ignore
## -0,0 +1,2 ##
+.deps
+.dirstamp
     
    66#define prefix (from http://www.gnu.org/software/autoconf/manual/autoconf-2.67/html_node/Defining-Directories.html)
    77AM_CPPFLAGS+=  -DISSM_PREFIX='"$(prefix)"'
    88
    9 js_scripts = ${ISSM_DIR}/src/wrappers/TriMesh/TriMesh.js  \
     9js_scripts = ${ISSM_DIR}/src/wrappers/Triangle/Triangle.js  \
    1010                         ${ISSM_DIR}/src/wrappers/NodeConnectivity/NodeConnectivity.js\
    1111                         ${ISSM_DIR}/src/wrappers/ContourToMesh/ContourToMesh.js\
    1212                         ${ISSM_DIR}/src/wrappers/ElementConnectivity/ElementConnectivity.js\
    … …  
    8080libISSMApi_la_LDFLAGS = -static
    8181endif
    8282
    83 IssmModule_SOURCES = ../TriMesh/TriMesh.cpp \
     83IssmModule_SOURCES = ../Triangle/Triangle.cpp \
    8484                                         ../NodeConnectivity/NodeConnectivity.cpp\
    8585                                         ../ContourToMesh/ContourToMesh.cpp\
    8686                                         ../ElementConnectivity/ElementConnectivity.cpp\
    … …  
    8888                                         ../IssmConfig/IssmConfig.cpp\
    8989                                         ../Issm/issm.cpp
    9090
    91 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
     91IssmModule_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
    9292IssmModule_LDADD = ${deps} $(TRIANGLELIB)  $(GSLLIB)
    9393#}}}

  • ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.cpp

     
     1/*!\file:  ProcessRifts.cpp
     2 * \brief split a mesh where a rift (or fault) is present
     3 */
     4
     5#include "./ProcessRifts.h"
     6
     7void ProcessRiftsUsage(void){/*{{{*/
     8        _printf_("\n");
     9        _printf_("   usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1) \n");
     10        _printf_("      where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
     11        _printf_("      index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
     12}/*}}}*/
     13WRAPPER(ProcessRifts_python){
     14
     15        /* returned quantities: */
     16        RiftStruct *riftstruct = NULL;
     17
     18        /* input: */
     19        int     nel,nods;
     20        int    *index          = NULL;
     21        double *x              = NULL;
     22        double *y              = NULL;
     23        int    *segments       = NULL;
     24        int    *segmentmarkers = NULL;
     25        int     num_seg;
     26
     27        /*Boot module*/
     28        MODULEBOOT();
     29
     30        /*checks on arguments on the matlab side: */
     31        CHECKARGUMENTS(NLHS,NRHS,&ProcessRiftsUsage);
     32
     33        /*Fetch data */
     34        FetchData(&index,&nel,NULL,INDEXIN);
     35        FetchData(&x,&nods,XIN);
     36        FetchData(&y,NULL,YIN);
     37        FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
     38        FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
     39
     40        /*call x layer*/
     41        ProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,&riftstruct);
     42
     43        /*Output : */
     44        WriteData(INDEXOUT,index,nel,3);
     45        WriteData(XOUT,x,nods,1);
     46        WriteData(YOUT,y,nods,1);
     47        WriteData(SEGMENTSOUT,segments,num_seg,3);
     48        WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
     49        WriteData(RIFTSTRUCT,riftstruct);
     50
     51        /*end module: */
     52        delete riftstruct;
     53        xDelete<int>(index);
     54        xDelete<double>(x);
     55        xDelete<double>(y);
     56        xDelete<int>(segments);
     57        xDelete<int>(segmentmarkers );
     58        MODULEEND();
     59}

  • ../trunk-jpl/src/wrappers/ProcessRifts/ProcessRifts.h

     
     1/*
     2 * ProcessRifts.h
     3 */
     4
     5#ifndef _PROCESSRIFTS_H_
     6#define _PROCESSRIFTS_H_
     7
     8#ifdef HAVE_CONFIG_H
     9#include <config.h>
     10#else
     11#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     12#endif
     13
     14/*For python modules: needs to come before header files inclusion*/
     15#ifdef _HAVE_PYTHON_
     16#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
     17#endif
     18
     19#include "../bindings.h"
     20#include "../../c/main/globals.h"
     21#include "../../c/modules/modules.h"
     22#include "../../c/shared/shared.h"
     23
     24#undef __FUNCT__
     25#define __FUNCT__  "ProcessRifts"
     26
     27#ifdef _HAVE_MATLAB_MODULES_
     28/* serial input macros: */
     29#define INDEXIN          prhs[0]
     30#define XIN              prhs[1]
     31#define YIN              prhs[2]
     32#define SEGMENTSIN       prhs[3]
     33#define SEGMENTMARKERSIN prhs[4]
     34/* serial output macros: */
     35#define INDEXOUT          (mxArray**)&plhs[0]
     36#define XOUT              (mxArray**)&plhs[1]
     37#define YOUT              (mxArray**)&plhs[2]
     38#define SEGMENTSOUT       (mxArray**)&plhs[3]
     39#define SEGMENTMARKERSOUT (mxArray**)&plhs[4]
     40#define RIFTSTRUCT        (mxArray**)&plhs[5]
     41#endif
     42
     43#ifdef _HAVE_PYTHON_MODULES_
     44/* serial input macros: */
     45#define INDEXIN          PyTuple_GetItem(args,0)
     46#define XIN              PyTuple_GetItem(args,1)
     47#define YIN              PyTuple_GetItem(args,2)
     48#define SEGMENTSIN       PyTuple_GetItem(args,3)
     49#define SEGMENTMARKERSIN PyTuple_GetItem(args,4)
     50/* serial output macros: */
     51#define INDEXOUT          output,0
     52#define XOUT              output,1
     53#define YOUT              output,2
     54#define SEGMENTSOUT       output,3
     55#define SEGMENTMARKERSOUT output,4
     56#define RIFTSTRUCT        output,5
     57#endif
     58
     59/* serial arg counts: */
     60#undef NLHS
     61#define NLHS  6
     62#undef NRHS
     63#define NRHS  5
     64
     65#endif

  • ../trunk-jpl/src/wrappers/ProcessRifts

  • ../trunk-jpl/src/wrappers/matlab/Makefile.am 

    Property changes on: ../trunk-jpl/src/wrappers/ProcessRifts
___________________________________________________________________
Added: svn:ignore
## -0,0 +1,2 ##
+.deps
+.dirstamp
     
    5757                                                 MeshProfileIntersection_matlab.la\
    5858                                                 PointCloudFindNeighbors_matlab.la\
    5959                                                 PropagateFlagsFromConnectivity_matlab.la\
    60                                                  TriMesh_matlab.la\
    61                                                  TriMeshProcessRifts_matlab.la\
     60                                                 Triangle_matlab.la\
     61                                                 ProcessRifts_matlab.la\
    6262                                                 Scotch_matlab.la
    6363
    6464if CHACO
    … …  
    232232ShpRead_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
    233233ShpRead_matlab_la_LIBADD = ${deps} $(SHAPELIBLIB) $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
    234234
    235 TriMesh_matlab_la_SOURCES = ../TriMesh/TriMesh.cpp
    236 TriMesh_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
    237 TriMesh_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(TRIANGLELIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
     235Triangle_matlab_la_SOURCES = ../Triangle/Triangle.cpp
     236Triangle_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
     237Triangle_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(TRIANGLELIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
    238238
    239 TriMeshProcessRifts_matlab_la_SOURCES = ../TriMeshProcessRifts/TriMeshProcessRifts.cpp
    240 TriMeshProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
    241 TriMeshProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
     239ProcessRifts_matlab_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp
     240ProcessRifts_matlab_la_CXXFLAGS = ${AM_CXXFLAGS}
     241ProcessRifts_matlab_la_LIBADD = ${deps} $(PETSCLIB) $(MPILIB) $(NEOPZLIB) $(GSLLIB) $(PROJ4LIB)
    242242#}}}

  • ../trunk-jpl/src/wrappers/python/Makefile.am

     
    4040                                                IssmConfig_python.la\
    4141                                                MeshProfileIntersection_python.la\
    4242                                                NodeConnectivity_python.la\
    43                                                 TriMesh_python.la\
    44                                                 TriMeshProcessRifts_python.la
     43                                                Triangle_python.la\
     44                                                ProcessRifts_python.la
    4545endif
    4646#}}}
    4747#Flags and libraries {{{
    … …  
    140140NodeConnectivity_python_la_CXXFLAGS = ${AM_CXXFLAGS}
    141141NodeConnectivity_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
    142142
    143 TriMesh_python_la_SOURCES = ../TriMesh/TriMesh.cpp
    144 TriMesh_python_la_CXXFLAGS = ${AM_CXXFLAGS}
    145 TriMesh_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)
     143Triangle_python_la_SOURCES = ../Triangle/Triangle.cpp
     144Triangle_python_la_CXXFLAGS = ${AM_CXXFLAGS}
     145Triangle_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(TRIANGLELIB) $(GSLLIB) $(PROJ4LIB)
    146146
    147 TriMeshProcessRifts_python_la_SOURCES = ../TriMeshProcessRifts/TriMeshProcessRifts.cpp
    148 TriMeshProcessRifts_python_la_CXXFLAGS = ${AM_CXXFLAGS}
    149 TriMeshProcessRifts_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
     147ProcessRifts_python_la_SOURCES = ../ProcessRifts/ProcessRifts.cpp
     148ProcessRifts_python_la_CXXFLAGS = ${AM_CXXFLAGS}
     149ProcessRifts_python_la_LIBADD = ${deps} $(MPILIB) $(PETSCLIB) $(GSLLIB) $(PROJ4LIB)
    150150#}}}

  • ../trunk-jpl/src/m/mesh/triangle.js

     
    11function triangle(md){
    22//TRIANGLE - create model mesh using the triangle package
    33//
    4 //   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
     4//   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
    55//   where md is a @model object, domainname is the name of an Argus domain outline file,
    66//   and resolution is a characteristic length for the mesh (same unit as the domain outline
    77//   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
    … …  
    3535        var area=Math.pow(resolution,2);
    3636
    3737        //Call mesher:
    38         var return_array=TriMesh(md, domain, rifts, area);
     38        var return_array=Triangle(md, domain, rifts, area);
    3939
    4040        //Plug into md:
    4141        md.mesh.elements=return_array[0];

  • ../trunk-jpl/src/m/mesh/triangle2dvertical.m

     
    11function md=triangle(md,domainname,resolution)
    22%TRIANGLE - create model mesh using the triangle package
    33%
    4 %   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
     4%   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
    55%   where md is a @model object, domainname is the name of an Argus domain outline file,
    66%   and resolution is a characteristic length for the mesh (same unit as the domain outline
    77%   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
    … …  
    2323
    2424area=resolution^2;
    2525
    26 %Mesh using TriMesh
    27 [elements,x,z,segments,segmentmarkers]=TriMesh(domainname,'',area);
     26%Mesh using Triangle
     27[elements,x,z,segments,segmentmarkers]=Triangle(domainname,'',area);
    2828
    2929%check that all the created nodes belong to at least one element
    3030orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));

  • ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py

     
    11import numpy as np
    2 from TriMeshProcessRifts import TriMeshProcessRifts
     2from ProcessRifts import ProcessRifts
    33from ContourToMesh import ContourToMesh
    44from meshprocessoutsiderifts import meshprocessoutsiderifts
    55from GetAreas import GetAreas
    … …  
    2121        """
    2222
    2323        #Call MEX file
    24         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)
     24        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)
    2525        md.mesh.elements=md.mesh.elements.astype(int)
    2626        md.mesh.x=md.mesh.x.reshape(-1)
    2727        md.mesh.y=md.mesh.y.reshape(-1)
    … …  
    2828        md.mesh.segments=md.mesh.segments.astype(int)
    2929        md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
    3030        if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
    31                 raise RuntimeError("TriMeshProcessRifts did not find any rift")
     31                raise RuntimeError("ProcessRifts did not find any rift")
    3232
    3333        #Fill in rest of fields:
    3434        numrifts=len(md.rifts.riftstruct)

  • ../trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m

     
    2424end
    2525
    2626%Call MEX file
    27 [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);
     27[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);
    2828if ~isstruct(md.rifts.riftstruct),
    29         error('TriMeshProcessRifts did not find any rift');
     29        error('ProcessRifts did not find any rift');
    3030end
    3131
    3232%Fill in rest of fields:

  • ../trunk-jpl/src/m/mesh/argusmesh.m

     
    88%      md=argusmesh(md,infile)
    99%
    1010%   Example:
    11 %     md=argusmesh(md,'TriMesh.exp')
     11%     md=argusmesh(md,'Domain.exp')
    1212
    1313%some argument check:
    1414if nargin~=2 | nargout~=1,

  • ../trunk-jpl/src/m/mesh/triangle.py

     
    11import os.path
    22import numpy as np
    33from mesh2d import mesh2d
    4 from TriMesh import TriMesh
     4from Triangle import Triangle
    55from NodeConnectivity import NodeConnectivity
    66from ElementConnectivity import ElementConnectivity
    77import MatlabFuncs as m
    … …  
    1010        """
    1111        TRIANGLE - create model mesh using the triangle package
    1212
    13            This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
     13           This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
    1414           where md is a @model object, domainname is the name of an Argus domain outline file,
    1515           and resolution is a characteristic length for the mesh (same unit as the domain outline
    1616           unit). Riftname is an optional argument (Argus domain outline) describing rifts.
    … …  
    4747        if not os.path.exists(domainname):
    4848                raise IOError("file '%s' not found" % domainname)
    4949
    50         #Mesh using TriMesh
     50        #Mesh using Triangle
    5151        md.mesh=mesh2d()
    52         md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=TriMesh(domainname,riftname,area)
     52        md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers=Triangle(domainname,riftname,area)
    5353        md.mesh.elements=md.mesh.elements.astype(int)
    5454        md.mesh.segments=md.mesh.segments.astype(int)
    5555        md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)

  • ../trunk-jpl/src/m/mesh/triangle.m

     
    11function md=triangle(md,domainname,varargin)
    22%TRIANGLE - create model mesh using the triangle package
    33%
    4 %   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
     4%   This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
    55%   where md is a @model object, domainname is the name of an Argus domain outline file,
    66%   and resolution is a characteristic length for the mesh (same unit as the domain outline
    77%   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
    … …  
    4141        error(['file "' domainname '" not found']);
    4242end
    4343
    44 %Mesh using TriMesh
    45 [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area);
     44%Mesh using Triangle
     45[elements,x,y,segments,segmentmarkers]=Triangle(domainname,riftname,area);
    4646
    4747%check that all the created nodes belong to at least one element
    4848removeorphans=1;

  • ../trunk-jpl/src/m/modules/TriMesh.m

     
    1 function [index,x,y,segments,segmentmarkers] = TriMesh(domainoutlinefilename,rifts,mesh_area);
    2 %TRIMESH - Mesh a domain using an .exp file
    3 %
    4 %   Usage:
    5 %     [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,mesh_area);
    6 %             
    7 %   index,x,y:  Defines a triangulation
    8 %   segments:   Array made of exterior segments to the mesh domain outline
    9 %   segmentmarkers:     Array flagging each segment
    10 %
    11 %   domainoutlinefilename:      Argus domain outline file
    12 %   mesh_area:  Maximum area desired for any element of the resulting mesh
    13 
    14 % Check usage
    15 if nargin~=3 && nargout~=5
    16         help TriMesh
    17         error('Wrong usage (see above)');
    18 end
    19 
    20 % Call mex module
    21 [index,x,y,segments,segmentmarkers]=TriMesh_matlab(domainoutlinefilename,rifts,mesh_area);

  • ../trunk-jpl/src/m/modules/TriMesh.py

     
    1 from TriMesh_python import TriMesh_python
    2 
    3 def TriMesh(domainoutlinefilename,rifts,mesh_area):
    4         """
    5         TRIMESH - Mesh a domain using an .exp file
    6 
    7            Usage:
    8                         [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,rifts,mesh_area);
    9 
    10            index,x,y: defines a triangulation
    11                 segments: An array made of exterior segments to the mesh domain outline
    12                 segmentmarkers: An array flagging each segment
    13 
    14            domainoutlinefilename: an Argus domain outline file
    15                 mesh_area: The maximum area desired for any element of the resulting mesh
    16         """
    17         # Call mex module
    18         index,x,y,segments,segmentmarkers=TriMesh_python(domainoutlinefilename,rifts,mesh_area)
    19         # Return
    20         return index,x,y,segments,segmentmarkers
    21 

  • ../trunk-jpl/src/m/modules/TriMeshProcessRifts.py

     
    1 from TriMeshProcessRifts_python import TriMeshProcessRifts_python
    2 
    3 def TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1):
    4         """
    5         TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
    6 
    7            Usage:
    8                    [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
    9 
    10            (index1,x1,y1,segments1,segmentmarkers1):    An initial triangulation.
    11            [index2,x2,y2,segments2,segmentmarkers2,rifts2]:     The resulting triangulation where rifts have been processed.
    12         """
    13         # Call mex module
    14         index2,x2,y2,segments2,segmentmarkers2,rifts2 = TriMeshProcessRifts_python(index1,x1,y1,segments1,segmentmarkers1)
    15         # Return
    16         return index2,x2,y2,segments2,segmentmarkers2,rifts2

  • ../trunk-jpl/src/m/modules/TriMeshProcessRifts.m

     
    1 function [index2,x2,y2,segments2,segmentmarkers2,rifts2] = TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
    2 %TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
    3 %
    4 %   Usage:
    5 %      [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
    6 %   
    7 %   (index1,x1,y1,segments1,segmentmarkers1):   An initial triangulation.
    8 %   [index2,x2,y2,segments2,segmentmarkers2,rifts2]:    The resulting triangulation where rifts have been processed.
    9 
    10 % Check usage
    11 if nargin~=5 && nargout~=6
    12         help TriMeshProcessRifts
    13         error('Wrong usage (see above)');
    14 end
    15 
    16 % Call mex module
    17 [index2,x2,y2,segments2,segmentmarkers2,rifts2] = TriMeshProcessRifts_matlab(index1,x1,y1,segments1,segmentmarkers1);

  • ../trunk-jpl/src/m/modules/ProcessRifts.m

     
     1function [index2,x2,y2,segments2,segmentmarkers2,rifts2] = ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
     2%TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
     3%
     4%   Usage:
     5%      [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
     6%   
     7%   (index1,x1,y1,segments1,segmentmarkers1):   An initial triangulation.
     8%   [index2,x2,y2,segments2,segmentmarkers2,rifts2]:    The resulting triangulation where rifts have been processed.
     9
     10% Check usage
     11if nargin~=5 && nargout~=6
     12        help ProcessRifts
     13        error('Wrong usage (see above)');
     14end
     15
     16% Call mex module
     17[index2,x2,y2,segments2,segmentmarkers2,rifts2] = ProcessRifts_matlab(index1,x1,y1,segments1,segmentmarkers1);

  • ../trunk-jpl/src/m/modules/Triangle.py

     
     1from Triangle_python import Triangle_python
     2
     3def Triangle(domainoutlinefilename,rifts,mesh_area):
     4        """
     5        TRIMESH - Mesh a domain using an .exp file
     6
     7           Usage:
     8                        [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,mesh_area);
     9
     10           index,x,y: defines a triangulation
     11                segments: An array made of exterior segments to the mesh domain outline
     12                segmentmarkers: An array flagging each segment
     13
     14           domainoutlinefilename: an Argus domain outline file
     15                mesh_area: The maximum area desired for any element of the resulting mesh
     16        """
     17        # Call mex module
     18        index,x,y,segments,segmentmarkers=Triangle_python(domainoutlinefilename,rifts,mesh_area)
     19        # Return
     20        return index,x,y,segments,segmentmarkers
     21

  • ../trunk-jpl/src/m/modules/Triangle.m

     
     1function [index,x,y,segments,segmentmarkers] = Triangle(domainoutlinefilename,rifts,mesh_area);
     2%TRIMESH - Mesh a domain using an .exp file
     3%
     4%   Usage:
     5%     [index,x,y,segments,segmentmarkers]=Triangle(domainoutlinefilename,rifts,mesh_area);
     6%             
     7%   index,x,y:  Defines a triangulation
     8%   segments:   Array made of exterior segments to the mesh domain outline
     9%   segmentmarkers:     Array flagging each segment
     10%
     11%   domainoutlinefilename:      Argus domain outline file
     12%   mesh_area:  Maximum area desired for any element of the resulting mesh
     13
     14% Check usage
     15if nargin~=3 && nargout~=5
     16        help Triangle
     17        error('Wrong usage (see above)');
     18end
     19
     20% Call mex module
     21[index,x,y,segments,segmentmarkers]=Triangle_matlab(domainoutlinefilename,rifts,mesh_area);

  • ../trunk-jpl/src/m/modules/ProcessRifts.py

     
     1from ProcessRifts_python import ProcessRifts_python
     2
     3def ProcessRifts(index1,x1,y1,segments1,segmentmarkers1):
     4        """
     5        TRIMESHPROCESSRIFTS - Split a mesh where a rift (or fault) is present
     6
     7           Usage:
     8                   [index2,x2,y2,segments2,segmentmarkers2,rifts2]=ProcessRifts(index1,x1,y1,segments1,segmentmarkers1);
     9
     10           (index1,x1,y1,segments1,segmentmarkers1):    An initial triangulation.
     11           [index2,x2,y2,segments2,segmentmarkers2,rifts2]:     The resulting triangulation where rifts have been processed.
     12        """
     13        # Call mex module
     14        index2,x2,y2,segments2,segmentmarkers2,rifts2 = ProcessRifts_python(index1,x1,y1,segments1,segmentmarkers1)
     15        # Return
     16        return index2,x2,y2,segments2,segmentmarkers2,rifts2
Note: See TracBrowser for help on using the repository browser.