Changeset 8301


Ignore:
Timestamp:
05/16/11 15:34:39 (14 years ago)
Author:
Mathieu Morlighem
Message:

moved all grids to nodes

Location:
issm/trunk/src/c/shared
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/shared/Exp/DomainOutlineRead.cpp

    r6412 r8301  
    11/*!\file:  DomainOutlineRead.cpp
    2  * \brief DomainOutlineRead.c: read the grid coordinates defined in a domain
     2 * \brief DomainOutlineRead.c: read the vertex coordinates defined in a domain
    33 * outline from Argus (.exp file). The first contour in the file is for
    44 * the outside domain outline. The following contours represent holes in
     
    1111#include "../Exceptions/exceptions.h"
    1212
    13 int DomainOutlineRead(int* pnprof,int** pprofngrids,double*** ppprofx,double*** ppprofy,char* domainname){
     13int DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,char* domainname){
    1414
    1515       
     
    2626        /*output: */
    2727        int nprof; //number of profiles in the domainname file
    28         int* profngrids=NULL; //array holding the number of grids for the nprof profiles
     28        int* profnvertices=NULL; //array holding the number of vertices for the nprof profiles
    2929        double** pprofx=NULL; //array of profiles x coordinates
    3030        double** pprofy=NULL; //array of profiles y coordinates
     
    6161       
    6262        /*Allocate and initialize all the profiles: */
    63         profngrids=(int*)xmalloc(nprof*sizeof(int));
     63        profnvertices=(int*)xmalloc(nprof*sizeof(int));
    6464        pprofx=(double**)xmalloc(nprof*sizeof(double*));
    6565        pprofy=(double**)xmalloc(nprof*sizeof(double*));
     
    8080                fscanf(fid,"%s %s %s %s\n",chardummy,chardummy,chardummy,chardummy);
    8181               
    82                 /*Get number of profile grids: */
     82                /*Get number of profile vertices: */
    8383                fscanf(fid,"%u %s\n",&n,chardummy);
    8484       
     
    8686                fscanf(fid,"%s %s %s %s %s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
    8787
    88                 /*Allocate grids: */
     88                /*Allocate vertices: */
    8989                x=(double*)xmalloc(n*sizeof(double));
    9090                y=(double*)xmalloc(n*sizeof(double));
    9191               
    9292
    93                 /*Read grids: */
     93                /*Read vertices: */
    9494                for (i=0;i<n;i++){
    9595                        fscanf(fid,"%lf %lf\n",x+i,y+i);
     
    102102
    103103                /*Assign pointers: */
    104                 profngrids[counter]=n;
     104                profnvertices[counter]=n;
    105105                pprofx[counter]=x;
    106106                pprofy[counter]=y;
     
    120120        /*Assign output pointers: */
    121121        *pnprof=nprof;
    122         *pprofngrids=profngrids;
     122        *pprofnvertices=profnvertices;
    123123        *ppprofx=pprofx;
    124124        *ppprofy=pprofy;
  • issm/trunk/src/c/shared/Exp/IsInPoly.cpp

    r8270 r8301  
    1515
    1616/*IsInPoly {{{1*/
    17 int IsInPoly(Vec in,double* xc,double* yc,int numgrids,double* x,double* y,int i0,int i1, int edgevalue){
     17int IsInPoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
    1818
    1919        int i;
     
    2626
    2727        /*Get extrema*/
    28         for (i=1;i<numgrids;i++){
     28        for (i=1;i<numvertices;i++){
    2929                if(xc[i]<xmin) xmin=xc[i];
    3030                if(xc[i]>xmax) xmax=xc[i];
     
    3333        }
    3434
    35         /*Go through all grids of the mesh:*/
     35        /*Go through all vertices of the mesh:*/
    3636        for (i=i0;i<i1;i++){
    3737
     
    3939                VecGetValues(in,1,&i,&value);
    4040                if (value){
    41                         /*this grid already is inside one of the contours, continue*/
     41                        /*this vertex already is inside one of the contours, continue*/
    4242                        continue;
    4343                }
    4444
    45                 /*pick up grid (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
     45                /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
    4646                x0=x[i]; y0=y[i];
    4747                if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
     
    4949                }
    5050                else{
    51                         value=pnpoly(numgrids,xc,yc,x0,y0,edgevalue);
     51                        value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
    5252                }
    5353                VecSetValues(in,1,&i,&value,INSERT_VALUES);
     
    5656}/*}}}*/
    5757/*IsOutsidePoly {{{1*/
    58 int IsOutsidePoly(Vec in,double* xc,double* yc,int numgrids,double* x,double* y,int i0,int i1, int edgevalue){
     58int IsOutsidePoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
    5959
    6060        int i,j;
     
    6767
    6868        /*Get extrema*/
    69         for (i=1;i<numgrids;i++){
     69        for (i=1;i<numvertices;i++){
    7070                if(xc[i]<xmin) xmin=xc[i];
    7171                if(xc[i]>xmax) xmax=xc[i];
     
    7474        }
    7575
    76         /*Go through all grids of the mesh:*/
     76        /*Go through all vertices of the mesh:*/
    7777        for (i=i0;i<i1;i++){
    7878
     
    8080                VecGetValues(in,1,&i,&value);
    8181                if (value){
    82                         /*this grid already is inside one of the contours, continue*/
     82                        /*this vertex already is inside one of the contours, continue*/
    8383                        continue;
    8484                }
    8585
    86                 /*pick up grid (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
     86                /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
    8787                x0=x[i]; y0=y[i];
    8888                if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
     
    9090                }
    9191                else{
    92                         value=1-pnpoly(numgrids,xc,yc,x0,y0,edgevalue);
     92                        value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
    9393                }
    9494                VecSetValues(in,1,&i,&value,INSERT_VALUES);
  • issm/trunk/src/c/shared/Exp/IsInPolySerial.cpp

    r5016 r8301  
    99
    1010
    11 int IsInPolySerial(double* in,double* xc,double* yc,int numgrids,double* x,double* y,int nods, int edgevalue){
     11int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue){
    1212
    1313        int i,j;
    1414        double x0,y0;
    1515
    16         /*Go through all grids of the mesh:*/
     16        /*Go through all vertices of the mesh:*/
    1717        for (i=0;i<nods;i++){
    1818                if (in[i]){
    19                         /*this grid already is inside one of the contours, continue*/
     19                        /*this vertex already is inside one of the contours, continue*/
    2020                        continue;
    2121                }
    22                 /*pick up grid: */
     22                /*pick up vertex: */
    2323                x0=x[i];
    2424                y0=y[i];
    25                 if (pnpoly(numgrids,xc,yc,x0,y0,edgevalue)){
     25                if (pnpoly(numvertices,xc,yc,x0,y0,edgevalue)){
    2626                        in[i]=1;
    2727                }
  • issm/trunk/src/c/shared/Exp/exp.h

    r8270 r8301  
    99#include "../../toolkits/toolkits.h"
    1010
    11 int IsInPoly(Vec in,double* xc,double* yc,int numgrids,double* x,double* y,int i0,int i1, int edgevalue);
    12 int IsOutsidePoly(Vec in,double* xc,double* yc,int numgrids,double* x,double* y,int i0,int i1, int edgevalue);
    13 int IsInPolySerial(double* in,double* xc,double* yc,int numgrids,double* x,double* y,int nods, int edgevalue);
    14 int DomainOutlineRead(int* pnprof,int** pprofngrids,double*** ppprofx,double*** ppprofy,char* domainname);
     11int IsInPoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue);
     12int IsOutsidePoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue);
     13int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue);
     14int DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,char* domainname);
    1515int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue);
    1616
  • issm/trunk/src/c/shared/TriMesh/AssociateSegmentToElement.cpp

    r3332 r8301  
    1414        double* segments=NULL;
    1515
    16         /*grid indices: */
     16        /*node indices: */
    1717        double A,B;
    1818
  • issm/trunk/src/c/shared/TriMesh/GridInsideHole.cpp

    r3332 r8301  
    1818        double yA,yB,yC,yD,yE;
    1919
    20         /*Take first and last grids: */
     20        /*Take first and last vertices: */
    2121        xA=x[0];
    2222        yA=y[0];
  • issm/trunk/src/c/shared/TriMesh/OrderSegments.cpp

    r3332 r8301  
    1414        double* segments=NULL;
    1515
    16         /*grid indices: */
     16        /*vertex indices: */
    1717        double A,B;
    1818        /*element indices: */
  • issm/trunk/src/c/shared/TriMesh/SplitMeshForRifts.cpp

    r5016 r8301  
    1717       
    1818        int i,j,k,l;
    19         int grid;
     19        int node;
    2020        int el;
    2121
     
    5050        /*Establish list of segments that belong to a rift: */
    5151        RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments); /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,
    52                                                                                                                                                                    first grid and second sgrid)*/
     52                                                                                                                                                                   first node and second snode)*/
    5353
    54         /*Go through all grids of the rift segments, and start splitting the mesh: */
    55         flags=(int*)xcalloc(nods,sizeof(int)); //to make sure we don't split the same grids twice!
     54        /*Go through all nodes of the rift segments, and start splitting the mesh: */
     55        flags=(int*)xcalloc(nods,sizeof(int)); //to make sure we don't split the same nodes twice!
    5656        for (i=0;i<nriftsegs;i++){
    5757                for (j=0;j<2;j++){
    5858       
    59                         grid=*(riftsegments+4*i+j+2);
    60                         if(flags[grid-1]){
    61                                 /*This grid was already split, skip:*/
     59                        node=*(riftsegments+4*i+j+2);
     60                        if(flags[node-1]){
     61                                /*This node was already split, skip:*/
    6262                                continue;
    6363                        }
    6464                        else{
    65                                 flags[grid-1]=1;
     65                                flags[node-1]=1;
    6666                        }
    6767
    68                         if(IsGridOnRift(riftsegments,nriftsegs,grid)){
     68                        if(IsGridOnRift(riftsegments,nriftsegs,node)){
    6969                       
    70                                 DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,grid,index,nel);
     70                                DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
    7171                       
    72                                 /*Summary: we have for grid, a list of elements (GridElementListOnOneSideOfRift, of size NumGridElementListOnOneSideOfRift) that all contain grid
    73                                  *and that are on the same side of the rift. For all these elements, we clone grid into another grid, and we swap all instances of grid in the triangulation
    74                                  *for those elements, to the new grid.*/
     72                                /*Summary: we have for node, a list of elements (GridElementListOnOneSideOfRift, of size NumGridElementListOnOneSideOfRift) that all contain node
     73                                 *and that are on the same side of the rift. For all these elements, we clone node into another node, and we swap all instances of node in the triangulation
     74                                 *for those elements, to the new node.*/
    7575                               
    76                                 //augment number of grids
     76                                //augment number of nodes
    7777                                nods=nods+1;
    78                                 //create new grid
     78                                //create new node
    7979                                x=(double*)xrealloc(x,nods*sizeof(double));
    8080                                y=(double*)xrealloc(y,nods*sizeof(double));
    81                                 x[nods-1]=x[grid-1]; //matlab indexing
    82                                 y[nods-1]=y[grid-1]; //matlab indexing
     81                                x[nods-1]=x[node-1]; //matlab indexing
     82                                y[nods-1]=y[node-1]; //matlab indexing
    8383
    84                                 //change elements owning this grid
     84                                //change elements owning this node
    8585                                for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
    8686                                        el=GridElementListOnOneSideOfRift[k];
    8787                                        for (l=0;l<3;l++){
    88                                                 if (*(index+3*el+l)==grid)*(index+3*el+l)=nods; //again, matlab indexing.
     88                                                if (*(index+3*el+l)==node)*(index+3*el+l)=nods; //again, matlab indexing.
    8989                                        }
    9090                                }
    91                         }// if(IsGridOnRift(riftsegments,nriftsegs,grid))
     91                        }// if(IsGridOnRift(riftsegments,nriftsegs,node))
    9292                } //for(j=0;j<2;j++)
    9393        } //for (i=0;i<nriftsegs;i++)
    9494
    95         /*update segments: they got modified completely by adding new grids.*/
     95        /*update segments: they got modified completely by adding new nodes.*/
    9696        UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs);
    9797
  • issm/trunk/src/c/shared/TriMesh/TriMeshUtils.cpp

    r6748 r8301  
    1111
    1212#define RIFTPENALTYPAIRSWIDTH 8
    13 int IsGridOnRift(int* riftsegments, int nriftsegs, int grid){
    14 
    15         /*Does this grid belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
     13int IsGridOnRift(int* riftsegments, int nriftsegs, int node){
     14
     15        /*Does this node belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift,
    1616         *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).*/
    1717       
     
    2323        for (i=0;i<nriftsegs;i++){
    2424                for (j=0;j<2;j++){
    25                         if ((*(riftsegments+4*i+2+j))==grid) count++;
     25                        if ((*(riftsegments+4*i+2+j))==node) count++;
    2626                }
    2727        }
     
    3535                               
    3636
    37 int GridElementsList(int** pGridElements, int* pNumGridElements,int grid,double * index,int nel){
    38 
    39         /*From a grid, recover all the elements that are connected to it: */
     37int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel){
     38
     39        /*From a node, recover all the elements that are connected to it: */
    4040        int i,j;
    4141        int noerr=1;
     
    4848
    4949        /*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own
    50          * the grid. We start by allocating GridElements with that size, and realloc
     50         * the node. We start by allocating GridElements with that size, and realloc
    5151         * more if needed.*/
    5252
     
    5757        for (i=0;i<nel;i++){
    5858                for (j=0;j<3;j++){
    59                         if ((int)*(index+3*i+j)==grid){
     59                        if ((int)*(index+3*i+j)==node){
    6060                                if (NumGridElements<=(current_size-1)){
    6161                                        GridElements[NumGridElements]=i;
     
    9191
    9292int IsNeighbor(int el1,int el2,double* index){
    93         /*From a triangulation held in index, figure out if elements 1 and 2 have two grids in common: */
     93        /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
    9494        int i,j;
    9595        int count=0;
     
    132132        int* riftsegments=NULL;
    133133        int* riftsegments_uncompressed=NULL;
    134         double element_grids[3];
     134        double element_nodes[3];
    135135
    136136        /*Allocate segmentflags: */
     
    142142        for (i=0;i<nsegs;i++){
    143143                el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
    144                 /*Temporarily set grids belonging to the segments to -1 in the triangulation index, and
     144                /*Temporarily set nodes belonging to the segments to -1 in the triangulation index, and
    145145                 *then  proceed to find another element that owns the segment. If we don't find it, we know
    146146                 *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
    147                 element_grids[0]=*(index+3*el+0);
    148                 element_grids[1]=*(index+3*el+1);
    149                 element_grids[2]=*(index+3*el+2);
     147                element_nodes[0]=*(index+3*el+0);
     148                element_nodes[1]=*(index+3*el+1);
     149                element_nodes[2]=*(index+3*el+2);
    150150
    151151                *(index+3*el+0)=-1;
     
    156156
    157157                /*Restore index: */
    158                 *(index+3*el+0)=element_grids[0];
    159                 *(index+3*el+1)=element_grids[1];
    160                 *(index+3*el+2)=element_grids[2];
     158                *(index+3*el+0)=element_nodes[0];
     159                *(index+3*el+1)=element_nodes[1];
     160                *(index+3*el+2)=element_nodes[2];
    161161
    162162                if (el2!=-1){
     
    195195******************************************************************************************************************************/
    196196
    197 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int grid,double* index,int nel){
     197int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel){
    198198
    199199        int noerr=1;
     
    208208        int* GridElementListOnOneSideOfRift=NULL;
    209209
    210         /*Build a list of all the elements connected to this grid: */
    211         GridElementsList(&GridElements,&NumGridElements,grid,index,nel);
     210        /*Build a list of all the elements connected to this node: */
     211        GridElementsList(&GridElements,&NumGridElements,node,index,nel);
    212212
    213213        /*Figure out the list of elements  that are on the same side of the rift. To do so, we start from one
     
    282282        segmentmarkerlist=(double*)xrealloc(segmentmarkerlist,(nsegs+nriftsegs)*sizeof(double));
    283283
    284         /*First, update the existing segments to the new grids :*/
     284        /*First, update the existing segments to the new nodes :*/
    285285        for (i=0;i<nriftsegs;i++){
    286286                el1=*(riftsegments+4*i+0);
     
    290290                                /*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
    291291                                 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
    292                                  *we can only rely on the position (x,y) of the rift grids to create a segment:*/
     292                                 *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
    293293                                for (k=0;k<3;k++){
    294294                                        if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+0)-1]) && (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
     
    379379                                   IsInPoly
    380380******************************************************************************************************************************/
    381 int IsInPoly(double* in,double* xc,double* yc,int numgrids,double* x,double* y,int nods){
     381int IsInPoly(double* in,double* xc,double* yc,int numnodes,double* x,double* y,int nods){
    382382
    383383        int i,j;
    384384        double x0,y0;
    385385
    386         /*Go through all grids of the mesh:*/
     386        /*Go through all nodes of the mesh:*/
    387387        for (i=0;i<nods;i++){
    388388                if (in[i]){
    389                         /*this grid already is inside one of the contours, continue*/
     389                        /*this node already is inside one of the contours, continue*/
    390390                        continue;
    391391                }
    392                 /*pick up grid: */
     392                /*pick up node: */
    393393                x0=x[i];
    394394                y0=y[i];
    395                 if (pnpoly(numgrids,xc,yc,x0,y0)){
     395                if (pnpoly(numnodes,xc,yc,x0,y0)){
    396396                        in[i]=1;
    397397                }
     
    524524        double* segments=NULL;
    525525        double* pairs=NULL;
    526         int     grid1,grid2,grid3,grid4;
     526        int     node1,node2,node3,node4;
    527527
    528528        riftsnumpairs=(int*)xmalloc(numrifts*sizeof(int));
     
    535535                for (j=0;j<numsegs;j++){
    536536                        *(pairs+2*j+0)=*(segments+3*j+2); //retrieve element to which this segment belongs.
    537                         grid1=(int)*(segments+3*j+0)-1; grid2=(int)*(segments+3*j+1)-1;
     537                        node1=(int)*(segments+3*j+0)-1; node2=(int)*(segments+3*j+1)-1;
    538538                        /*Find element facing on other side of rift: */
    539539                        for (k=0;k<numsegs;k++){
    540540                                if (k==j)continue;
    541                                 grid3=(int)*(segments+3*k+0)-1; grid4=(int)*(segments+3*k+1)-1;
    542                                 /*We are trying to find 2 elements, where position of grid3 == position of grid1, and position of grid4 == position of grid2*/
    543                                 if (   (x[grid3]==x[grid1]) && (y[grid3]==y[grid1]) && (x[grid4]==x[grid2]) && (y[grid4]==y[grid2])
    544                                     || (x[grid3]==x[grid2]) && (y[grid3]==y[grid2]) && (x[grid4]==x[grid1]) && (y[grid4]==y[grid1])  ){
     541                                node3=(int)*(segments+3*k+0)-1; node4=(int)*(segments+3*k+1)-1;
     542                                /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
     543                                if (   (x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2])
     544                                    || (x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1])  ){
    545545                                        /*We found the corresponding element: */
    546546                                        *(pairs+2*j+1)=*(segments+3*k+2);
     
    585585        double* riftsegments=NULL;
    586586        double* riftpairs=NULL;
    587         int     grid1,grid2,grid3,grid4,temp_grid;
     587        int     node1,node2,node3,node4,temp_node;
    588588        double  el1,el2;
    589         int     newnods; //temporary # grid counter.
     589        int     newnods; //temporary # node counter.
    590590        double  xmin,ymin;
    591591        double* xreal=NULL;
    592592        double* yreal=NULL;
    593         int* grids=NULL;
    594         int* merginggrids=NULL;
     593        int* nodes=NULL;
     594        int* mergingnodes=NULL;
    595595        int     max_size;
    596596        int     redundant;
     
    608608        newnods=nods;
    609609
    610         /*Figure out a unique value to flag x and y for grid removal: */
     610        /*Figure out a unique value to flag x and y for node removal: */
    611611        xmin=x[0];
    612612        ymin=y[0];
     
    618618        ymin=ymin-dabs(ymin);
    619619
    620         /*Initialize two arrays, one for grids that are going to be merged, the other with corresponding grids being merge into: */
     620        /*Initialize two arrays, one for nodes that are going to be merged, the other with corresponding nodes being merge into: */
    621621        max_size=0;
    622622        for (i=0;i<numrifts1;i++){
    623623                max_size+=rifts1numsegs[i];
    624624        }
    625         grids=(int*)xmalloc(max_size*sizeof(int));
    626         merginggrids=(int*)xmalloc(max_size*sizeof(int));
    627 
    628         /*Go through the rifts segments, and identify which grid we are going to merge with its counterpart on the other side
    629          *of the rift. The way we identify this grid is by looking at the element pairs, and the corresponding grids: */
     625        nodes=(int*)xmalloc(max_size*sizeof(int));
     626        mergingnodes=(int*)xmalloc(max_size*sizeof(int));
     627
     628        /*Go through the rifts segments, and identify which node we are going to merge with its counterpart on the other side
     629         *of the rift. The way we identify this node is by looking at the element pairs, and the corresponding nodes: */
    630630        counter=0;
    631631        for (i=0;i<numrifts1;i++){
     
    635635                        el1=*(riftpairs+2*j+0);
    636636                        el2=*(riftpairs+2*j+1);
    637                         grid1=(int)*(riftsegments+3*j+0);
    638                         grid2=(int)*(riftsegments+3*j+1);
    639                         /*Summary, el1 and el2 are facing one another across the rift. grid1 and grid2 belong to el1 and
    640                          *are located on the rift. Find grid3 and grid4, grids belonging to el2 and located on the rift: */
     637                        node1=(int)*(riftsegments+3*j+0);
     638                        node2=(int)*(riftsegments+3*j+1);
     639                        /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
     640                         *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
    641641                        for (k=0;k<rifts1numsegs[i];k++){
    642642                                if (*(riftsegments+3*k+2)==el2){
    643                                         grid3=(int)*(riftsegments+3*k+0);
    644                                         grid4=(int)*(riftsegments+3*k+1);
     643                                        node3=(int)*(riftsegments+3*k+0);
     644                                        node4=(int)*(riftsegments+3*k+1);
    645645                                        break;
    646646                                }
    647647                        }
    648                         /* Make sure grid3 faces grid1 and grid4 faces grid2: */
    649                         if ((x[grid1]==x[grid4]) && (y[grid1]==y[grid4])){
    650                                 /*Swap grid3 and grid4:*/
    651                                 temp_grid=grid3;
    652                                 grid3=grid4;
    653                                 grid4=temp_grid;
    654                         }
    655                         /* Is any of these two grid pairs on the tip of a rift, in which case, we don't include it in grids or merginggrids: */
    656                         if ((grid1==grid3) || (grid2==grid4)){
    657                                 if(grid1!=grid3){
    658                                         /*Add grid1 and grid3 to grids and merginggrids if they have not already been added: */
     648                        /* Make sure node3 faces node1 and node4 faces node2: */
     649                        if ((x[node1]==x[node4]) && (y[node1]==y[node4])){
     650                                /*Swap node3 and node4:*/
     651                                temp_node=node3;
     652                                node3=node4;
     653                                node4=temp_node;
     654                        }
     655                        /* Is any of these two node pairs on the tip of a rift, in which case, we don't include it in nodes or mergingnodes: */
     656                        if ((node1==node3) || (node2==node4)){
     657                                if(node1!=node3){
     658                                        /*Add node1 and node3 to nodes and mergingnodes if they have not already been added: */
    659659                                        redundant=0;
    660660                                        for (k=0;k<counter;k++){
    661                                                 if ((merginggrids[k]==grid1) || (grids[k]==grid1))redundant=1;
     661                                                if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;
    662662                                        }
    663663                                        if(!redundant){
    664                                                 /*Ok, add grid1 to grids, and grid3 to merginggrids: */
    665                                                 grids[counter]=grid1;
    666                                                 merginggrids[counter]=grid3;
     664                                                /*Ok, add node1 to nodes, and node3 to mergingnodes: */
     665                                                nodes[counter]=node1;
     666                                                mergingnodes[counter]=node3;
    667667                                                counter++;
    668668                                        }
    669669                                }
    670                                 if(grid2!=grid4){
    671                                         /*Add grid2 and grid4 to grids and merginggrids if they have not already been added: */
     670                                if(node2!=node4){
     671                                        /*Add node2 and node4 to nodes and mergingnodes if they have not already been added: */
    672672                                        redundant=0;
    673673                                        for (k=0;k<counter;k++){
    674                                                 if ((merginggrids[k]==grid2) || (grids[k]==grid2))redundant=1;
     674                                                if ((mergingnodes[k]==node2) || (nodes[k]==node2))redundant=1;
    675675                                        }
    676676                                        if(!redundant){
    677                                                 /*Ok, add grid2 to grids, and grid4 to merginggrids: */
    678                                                 grids[counter]=grid2;
    679                                                 merginggrids[counter]=grid4;
     677                                                /*Ok, add node2 to nodes, and node4 to mergingnodes: */
     678                                                nodes[counter]=node2;
     679                                                mergingnodes[counter]=node4;
    680680                                                counter++;
    681681                                        }
     
    683683                        }
    684684                        else{
    685                                 /*Check that grid1 is not already present in the merginggrids: */
     685                                /*Check that node1 is not already present in the mergingnodes: */
    686686                                redundant=0;
    687687                                for (k=0;k<counter;k++){
    688                                         if ((merginggrids[k]==grid1) || (grids[k]==grid1))redundant=1;
     688                                        if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;
    689689                                }
    690690                                if(!redundant){
    691                                         /*Ok, add grid1 to grids, and grid3 to merginggrids: */
    692                                         grids[counter]=grid1;
    693                                         merginggrids[counter]=grid3;
     691                                        /*Ok, add node1 to nodes, and node3 to mergingnodes: */
     692                                        nodes[counter]=node1;
     693                                        mergingnodes[counter]=node3;
    694694                                        counter++;
    695695                                }
    696                                 /*Check that grid2 is not already present in the merginggrids: */
     696                                /*Check that node2 is not already present in the mergingnodes: */
    697697                                redundant=0;
    698698                                for (k=0;k<counter;k++){
    699                                         if ((merginggrids[k]==grid1) || (grids[k]==grid1))redundant=1;
     699                                        if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;
    700700                                }
    701701                                if(!redundant){
    702                                         /*Ok, add grid2 to grids, and grid4 to merginggrids: */
    703                                         grids[counter]=grid2;
    704                                         merginggrids[counter]=grid4;
     702                                        /*Ok, add node2 to nodes, and node4 to mergingnodes: */
     703                                        nodes[counter]=node2;
     704                                        mergingnodes[counter]=node4;
    705705                                        counter++;
    706706                                }
     
    709709        }
    710710
    711         /*Ok, we have counter pairs of grids (grids and merginggrids): start merging grids in the triangulation: */
     711        /*Ok, we have counter pairs of nodes (nodes and mergingnodes): start merging nodes in the triangulation: */
    712712        newnods=nods;
    713713        for (i=0;i<counter;i++){
    714                 grid1=grids[i];
    715                 grid3=merginggrids[i];
    716                 /*Merge grid3 into grid1: */
    717                 x[grid3]=xmin; //flag  for later removal from x
    718                 y[grid3]=ymin; //flag  for later removal from y
     714                node1=nodes[i];
     715                node3=mergingnodes[i];
     716                /*Merge node3 into node1: */
     717                x[node3]=xmin; //flag  for later removal from x
     718                y[node3]=ymin; //flag  for later removal from y
    719719                newnods--;
    720720                for (k=0;k<nel;k++){
    721                         if (*(index+3*k+0)==grid3)*(index+3*k+0)=grid1;
    722                         if (*(index+3*k+1)==grid3)*(index+3*k+1)=grid1;
    723                         if (*(index+3*k+2)==grid3)*(index+3*k+2)=grid1;
     721                        if (*(index+3*k+0)==node3)*(index+3*k+0)=node1;
     722                        if (*(index+3*k+1)==node3)*(index+3*k+1)=node1;
     723                        if (*(index+3*k+2)==node3)*(index+3*k+2)=node1;
    724724                }
    725725        }
     
    805805        double* riftpairs_copy=NULL;
    806806
    807         /*grid and element manipulation: */
    808         int     grid1,grid2,grid3,grid4,temp_grid,tip1,tip2,grid;
     807        /*node and element manipulation: */
     808        int     node1,node2,node3,node4,temp_node,tip1,tip2,node;
    809809        double  el1,el2;
    810810        int     already_ordered=0;
     
    830830                order=(int*)xmalloc(numsegs*sizeof(int));
    831831
    832                 /*First find the tips, using the pairs. If a pair of elements has one grid in common, this grid is a rift tip: */
     832                /*First find the tips, using the pairs. If a pair of elements has one node in common, this node is a rift tip: */
    833833                tip1=-1;
    834834                tip2=-1;
     
    837837                        el1=*(riftpairs+2*j+0);
    838838                        el2=*(riftpairs+2*j+1);
    839                         grid1=(int)*(riftsegments+3*j+0);
    840                         grid2=(int)*(riftsegments+3*j+1);
    841                         /*Summary, el1 and el2 are facing one another across the rift. grid1 and grid2 belong to el1 and
    842                          *are located on the rift. Find grid3 and grid4, grids belonging to el2 and located on the rift: */
     839                        node1=(int)*(riftsegments+3*j+0);
     840                        node2=(int)*(riftsegments+3*j+1);
     841                        /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
     842                         *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
    843843                        for (k=0;k<numsegs;k++){
    844844                                if (*(riftsegments+3*k+2)==el2){
    845                                         grid3=(int)*(riftsegments+3*k+0);
    846                                         grid4=(int)*(riftsegments+3*k+1);
     845                                        node3=(int)*(riftsegments+3*k+0);
     846                                        node4=(int)*(riftsegments+3*k+1);
    847847                                        break;
    848848                                }
    849849                        }
    850                         /* Make sure grid3 faces grid1 and grid4 faces grid2: */
    851                         if ((x[grid1]==x[grid4]) && (y[grid1]==y[grid4])){
    852                                 /*Swap grid3 and grid4:*/
    853                                 temp_grid=grid3;
    854                                 grid3=grid4;
    855                                 grid4=temp_grid;
     850                        /* Make sure node3 faces node1 and node4 faces node2: */
     851                        if ((x[node1]==x[node4]) && (y[node1]==y[node4])){
     852                                /*Swap node3 and node4:*/
     853                                temp_node=node3;
     854                                node3=node4;
     855                                node4=temp_node;
    856856                        }
    857857
    858858                        /*Figure out if a tip is on this element: */
    859                         if (grid3==grid1){
    860                                 /*grid1 is a tip*/
     859                        if (node3==node1){
     860                                /*node1 is a tip*/
    861861                                if (tip1==-1) {
    862                                         tip1=grid1;
     862                                        tip1=node1;
    863863                                        continue;
    864864                                }
    865                                 if ((tip2==-1) && (grid1!=tip1)){
    866                                         tip2=grid1;
     865                                if ((tip2==-1) && (node1!=tip1)){
     866                                        tip2=node1;
    867867                                        break;
    868868                                }
    869869                        }
    870870               
    871                         if (grid4==grid2){
    872                                 /*grid2 is a tip*/
     871                        if (node4==node2){
     872                                /*node2 is a tip*/
    873873                                if (tip1==-1){
    874                                         tip1=grid2;
     874                                        tip1=node2;
    875875                                        continue;
    876876                                }
    877                                 if ((tip2==-1) && (grid2!=tip1)){
    878                                         tip2=grid2;
     877                                if ((tip2==-1) && (node2!=tip1)){
     878                                        tip2=node2;
    879879                                        break;
    880880                                }
     
    889889                /*We have the two tips for this rift.  Go from tip1 to tip2, and figure out the order in which segments are sequential.
    890890                 *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
    891                 grid=tip1;
     891                node=tip1;
    892892                for (counter=0;counter<numsegs;counter++){
    893893                        for (j=0;j<numsegs;j++){
    894                                 grid1=(int)*(riftsegments+3*j+0);
    895                                 grid2=(int)*(riftsegments+3*j+1);
     894                                node1=(int)*(riftsegments+3*j+0);
     895                                node2=(int)*(riftsegments+3*j+1);
    896896                               
    897                                 if ((grid1==grid) || (grid2==grid)){
    898                                         /*Ok, this segment is connected to grid, plug its index into order, unless we already plugged it before: */
     897                                if ((node1==node) || (node2==node)){
     898                                        /*Ok, this segment is connected to node, plug its index into order, unless we already plugged it before: */
    899899                                        already_ordered=0;
    900900                                        for (k=0;k<counter;k++){
     
    906906                                        if (!already_ordered){
    907907                                                order[counter]=j;
    908                                                 if(grid1==grid){
    909                                                         grid=grid2;
     908                                                if(node1==node){
     909                                                        node=node2;
    910910                                                }
    911                                                 else if(grid2==grid){
    912                                                         grid=grid1;
     911                                                else if(node2==node){
     912                                                        node=node1;
    913913                                                }
    914914                                                break;
     
    958958        int i,j,k,k0;
    959959
    960         double el1,el2,grid1,grid2,grid3,grid4;
    961         double tip1,tip2,temp_grid;
     960        double el1,el2,node1,node2,node3,node4;
     961        double tip1,tip2,temp_node;
    962962
    963963        /*output: */
     
    994994                        el1=*(riftpairs+2*j+0);
    995995                        el2=*(riftpairs+2*j+1);
    996                         grid1=*(riftsegments+3*j+0);
    997                         grid2=*(riftsegments+3*j+1);
    998                         /*Find segment index to recover grid3 and grid4, facing grid1 and grid2: */
     996                        node1=*(riftsegments+3*j+0);
     997                        node2=*(riftsegments+3*j+1);
     998                        /*Find segment index to recover node3 and node4, facing node1 and node2: */
    999999                        k0=-1;
    10001000                        for(k=0;k<numsegs;k++){
     
    10041004                                }
    10051005                        }
    1006                         grid3=*(riftsegments+3*k0+0);
    1007                         grid4=*(riftsegments+3*k0+1);
    1008 
    1009                         /* Make sure grid3 faces grid1 and grid4 faces grid2: */
    1010                         if ((x[(int)grid1-1]==x[(int)grid4-1]) && (y[(int)grid1-1]==y[(int)grid4-1])){
    1011                                 /*Swap grid3 and grid4:*/
    1012                                 temp_grid=grid3;
    1013                                 grid3=grid4;
    1014                                 grid4=temp_grid;
     1006                        node3=*(riftsegments+3*k0+0);
     1007                        node4=*(riftsegments+3*k0+1);
     1008
     1009                        /* Make sure node3 faces node1 and node4 faces node2: */
     1010                        if ((x[(int)node1-1]==x[(int)node4-1]) && (y[(int)node1-1]==y[(int)node4-1])){
     1011                                /*Swap node3 and node4:*/
     1012                                temp_node=node3;
     1013                                node3=node4;
     1014                                node4=temp_node;
    10151015                        }       
    1016                         /*Ok, we have grid1 facing grid3, and grid2 facing grid4. Compute the normal to
     1016                        /*Ok, we have node1 facing node3, and node2 facing node4. Compute the normal to
    10171017                         *this segment, and its length: */
    1018                         normal[0]=cos(atan2(x[(int)grid1-1]-x[(int)grid2-1],y[(int)grid2-1]-y[(int)grid1-1]));
    1019                         normal[1]=sin(atan2(x[(int)grid1-1]-x[(int)grid2-1],y[(int)grid2-1]-y[(int)grid1-1]));
    1020                         length=sqrt(pow(x[(int)grid2-1]-x[(int)grid1-1],(double)2)+pow(y[(int)grid2-1]-y[(int)grid1-1],(double)2));
    1021 
    1022                         /*Be careful here, we want penalty loads on each grid, not on each segment. This means we cannot plug grid1,
    1023                          * grid2, grid3 and grid4 directly into riftpenaltypairs. We need to include grid1, grid2, grid3 and grid4,
     1018                        normal[0]=cos(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
     1019                        normal[1]=sin(atan2(x[(int)node1-1]-x[(int)node2-1],y[(int)node2-1]-y[(int)node1-1]));
     1020                        length=sqrt(pow(x[(int)node2-1]-x[(int)node1-1],(double)2)+pow(y[(int)node2-1]-y[(int)node1-1],(double)2));
     1021
     1022                        /*Be careful here, we want penalty loads on each node, not on each segment. This means we cannot plug node1,
     1023                         * node2, node3 and node4 directly into riftpenaltypairs. We need to include node1, node2, node3 and node4,
    10241024                         * only once. We'll add the normals and the lengths : */
    10251025
    1026                         if(grid1!=grid3){ //exclude tips from loads
     1026                        if(node1!=node3){ //exclude tips from loads
    10271027                                k1=-1;
    10281028                                for(k=0;k<counter;k++){
    1029                                         if( (*(riftpenaltypairs+k*7+0))==grid1){
     1029                                        if( (*(riftpenaltypairs+k*7+0))==node1){
    10301030                                                k1=k;
    10311031                                                break;
     
    10331033                                }
    10341034                                if(k1==-1){
    1035                                         *(riftpenaltypairs+counter*7+0)=grid1;
    1036                                         *(riftpenaltypairs+counter*7+1)=grid3;
     1035                                        *(riftpenaltypairs+counter*7+0)=node1;
     1036                                        *(riftpenaltypairs+counter*7+1)=node3;
    10371037                                        *(riftpenaltypairs+counter*7+2)=el1;
    10381038                                        *(riftpenaltypairs+counter*7+3)=el2;
     
    10481048                                }
    10491049                        }
    1050                         if(grid2!=grid4){
     1050                        if(node2!=node4){
    10511051                                k2=-1;
    10521052                                for(k=0;k<counter;k++){
    1053                                         if( (*(riftpenaltypairs+k*7+0))==grid2){
     1053                                        if( (*(riftpenaltypairs+k*7+0))==node2){
    10541054                                                k2=k;
    10551055                                                break;
     
    10571057                                }
    10581058                                if(k2==-1){
    1059                                         *(riftpenaltypairs+counter*7+0)=grid2;
    1060                                         *(riftpenaltypairs+counter*7+1)=grid4;
     1059                                        *(riftpenaltypairs+counter*7+0)=node2;
     1060                                        *(riftpenaltypairs+counter*7+1)=node4;
    10611061                                        *(riftpenaltypairs+counter*7+2)=el1;
    10621062                                        *(riftpenaltypairs+counter*7+3)=el2;
     
    11011101        int noerr=1;
    11021102        int i,j,k;
    1103         double grid1,grid2,grid3;
     1103        double node1,node2,node3;
    11041104        int el;
    11051105
     
    11231123
    11241124        for (i=0;i<num_seg;i++){
    1125                 grid1=*(segments+3*i+0);
    1126                 grid2=*(segments+3*i+1);
    1127                 /*Find all elements connected to [grid1 grid2]: */
     1125                node1=*(segments+3*i+0);
     1126                node2=*(segments+3*i+1);
     1127                /*Find all elements connected to [node1 node2]: */
    11281128                pair_count=0;
    11291129                for (j=0;j<nel;j++){
    1130                         if (*(index+3*j+0)==grid1){
    1131                                 if ((*(index+3*j+1)==grid2) || (*(index+3*j+2)==grid2)){
     1130                        if (*(index+3*j+0)==node1){
     1131                                if ((*(index+3*j+1)==node2) || (*(index+3*j+2)==node2)){
    11321132                                        pair[pair_count]=j;
    11331133                                        pair_count++;
    11341134                                }
    11351135                        }
    1136                         if (*(index+3*j+1)==grid1){
    1137                                 if ((*(index+3*j+0)==grid2) || (*(index+3*j+2)==grid2)){
     1136                        if (*(index+3*j+1)==node1){
     1137                                if ((*(index+3*j+0)==node2) || (*(index+3*j+2)==node2)){
    11381138                                        pair[pair_count]=j;
    11391139                                        pair_count++;
    11401140                                }
    11411141                        }
    1142                         if (*(index+3*j+2)==grid1){
    1143                                 if ((*(index+3*j+0)==grid2) || (*(index+3*j+1)==grid2)){
     1142                        if (*(index+3*j+2)==node1){
     1143                                if ((*(index+3*j+0)==node2) || (*(index+3*j+1)==node2)){
    11441144                                        pair[pair_count]=j;
    11451145                                        pair_count++;
     
    11481148                }
    11491149                /*Ok, we have pair_count elements connected to this segment. For each of these elements,
    1150                  *figure out if the third grid also belongs to a segment: */
     1150                 *figure out if the third node also belongs to a segment: */
    11511151                if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to  2 elements
    11521152                        continue;
     
    11561156                                el=(int)pair[j];
    11571157                                triple=0;
    1158                                 /*First find grid3: */
    1159                                 if (*(index+3*el+0)==grid1){
    1160                                         if (*(index+3*el+1)==grid2)grid3=*(index+3*el+2);
    1161                                         else grid3=*(index+3*el+1);
    1162                                 }
    1163                                 if (*(index+3*el+1)==grid1){
    1164                                         if (*(index+3*el+0)==grid2)grid3=*(index+3*el+2);
    1165                                         else grid3=*(index+3*el+0);
    1166                                 }
    1167                                 if (*(index+3*el+2)==grid1){
    1168                                         if (*(index+3*el+0)==grid2)grid3=*(index+3*el+1);
    1169                                         else grid3=*(index+3*el+0);
    1170                                 }
    1171                                 /*Ok, we have grid3. Does grid3 belong to a segment? : */
     1158                                /*First find node3: */
     1159                                if (*(index+3*el+0)==node1){
     1160                                        if (*(index+3*el+1)==node2)node3=*(index+3*el+2);
     1161                                        else node3=*(index+3*el+1);
     1162                                }
     1163                                if (*(index+3*el+1)==node1){
     1164                                        if (*(index+3*el+0)==node2)node3=*(index+3*el+2);
     1165                                        else node3=*(index+3*el+0);
     1166                                }
     1167                                if (*(index+3*el+2)==node1){
     1168                                        if (*(index+3*el+0)==node2)node3=*(index+3*el+1);
     1169                                        else node3=*(index+3*el+0);
     1170                                }
     1171                                /*Ok, we have node3. Does node3 belong to a segment? : */
    11721172                                for (k=0;k<num_seg;k++){
    1173                                         if ((grid3==*(segments+3*k+0)) || (grid3==*(segments+3*k+1))){
     1173                                        if ((node3==*(segments+3*k+0)) || (node3==*(segments+3*k+1))){
    11741174                                                triple=1;
    11751175                                                break;
     
    11801180                                        x=(double*)xrealloc(x,(nods+1)*sizeof(double));
    11811181                                        y=(double*)xrealloc(y,(nods+1)*sizeof(double));
    1182                                         x[nods]=(x[(int)grid1-1]+x[(int)grid2-1]+x[(int)grid3-1])/3;
    1183                                         y[nods]=(y[(int)grid1-1]+y[(int)grid2-1]+y[(int)grid3-1])/3;
     1182                                        x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
     1183                                        y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
    11841184
    11851185                                        index=(double*)xrealloc(index,(nel+2)*3*sizeof(double));
    11861186                                        /*First, reassign element el: */
    1187                                         *(index+3*el+0)=grid1;
    1188                                         *(index+3*el+1)=grid2;
     1187                                        *(index+3*el+0)=node1;
     1188                                        *(index+3*el+1)=node2;
    11891189                                        *(index+3*el+2)=nods+1;
    11901190                                        /*Other two elements: */
    1191                                         *(index+3*nel+0)=grid2;
    1192                                         *(index+3*nel+1)=grid3;
     1191                                        *(index+3*nel+0)=node2;
     1192                                        *(index+3*nel+1)=node3;
    11931193                                        *(index+3*nel+2)=nods+1;
    11941194
    1195                                         *(index+3*(nel+1)+0)=grid3;
    1196                                         *(index+3*(nel+1)+1)=grid1;
     1195                                        *(index+3*(nel+1)+0)=node3;
     1196                                        *(index+3*(nel+1)+1)=node1;
    11971197                                        *(index+3*(nel+1)+2)=nods+1;
    11981198                                        /*we need  to change the segment elements corresponding to el: */
    11991199                                        for (k=0;k<num_seg;k++){
    12001200                                                if (*(segments+3*k+2)==(double)(el+1)){
    1201                                                         if ( ((*(segments+3*k+0)==grid1) && (*(segments+3*k+1)==grid2)) || ((*(segments+3*k+0)==grid2) && (*(segments+3*k+1)==grid1))) *(segments+3*k+2)=(double)(el+1);
    1202                                                         if ( ((*(segments+3*k+0)==grid2) && (*(segments+3*k+1)==grid3)) || ((*(segments+3*k+0)==grid3) && (*(segments+3*k+1)==grid2))) *(segments+3*k+2)=(double)(nel+1);
    1203                                                         if ( ((*(segments+3*k+0)==grid3) && (*(segments+3*k+1)==grid1)) || ((*(segments+3*k+0)==grid1) && (*(segments+3*k+1)==grid3))) *(segments+3*k+2)=(double)(nel+2);
     1201                                                        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)=(double)(el+1);
     1202                                                        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)=(double)(nel+1);
     1203                                                        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)=(double)(nel+2);
    12041204                                                }
    12051205                                        }
  • issm/trunk/src/c/shared/TriMesh/trimesh.h

    r1439 r8301  
    2222int SplitMeshForRifts(int* pnel,double** pindex,int* pnods,double** px,double** py,int* pnsegs,double** psegments,double** psegmentmarkerlist);
    2323
    24 int IsGridOnRift(int* riftsegments, int nriftsegs, int grid);
    25 int GridElementsList(int** pGridElements, int* pNumGridElements,int grid,double * index,int nel);
     24int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
     25int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
    2626int IsNeighbor(int el1,int el2,double* index);
    2727int IsOnRift(int el,int nriftsegs,int* riftsegments);
    2828int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments);
    29 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int grid,double* index,int nel);
     29int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel);
    3030int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs);
    3131int pnpoly(int npol, double *xp, double *yp, double x, double y);
Note: See TracChangeset for help on using the changeset viewer.