Changeset 12080


Ignore:
Timestamp:
04/20/12 13:49:38 (13 years ago)
Author:
Mathieu Morlighem
Message:

Reorganizing and adding assertion for safety

Location:
issm/trunk-jpl/src/c/shared/TriMesh
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r8301 r12080  
    5757                for (j=0;j<2;j++){
    5858       
    59                         node=*(riftsegments+4*i+j+2);
     59                        node=riftsegments[4*i+j+2];
    6060                        if(flags[node-1]){
    6161                                /*This node was already split, skip:*/
    … …  
    9494
    9595        /*update segments: they got modified completely by adding new nodes.*/
    96         UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs);
     96        UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);
    9797
    9898        /*Assign output pointers: */

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

    r12070 r12080  
    1111
    1212#define RIFTPENALTYPAIRSWIDTH 8
     13/*FUNCTION IsGridOnRift{{{*/
    1314int IsGridOnRift(int* riftsegments, int nriftsegs, int node){
    1415
    … …  
    3233                return 0;
    3334        }
    34 }
    35                                
    36 
     35}/*}}}*/
     36/*FUNCTION GridElementsList{{{*/
    3737int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel){
    3838
    … …  
    8787        *pNumGridElements=NumGridElements;
    8888        return noerr;
    89 }
    90 
    91 
     89}/*}}}*/
     90/*FUNCTION IsNeighbor{{{*/
    9291int IsNeighbor(int el1,int el2,double* index){
    9392        /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
    … …  
    105104                return 0;
    106105        }
    107 }
    108                                                        
    109 
     106}/*}}}*/
     107/*FUNCTION IsOnRift{{{*/
    110108int IsOnRift(int el,int nriftsegs,int* riftsegments){
    111109        /*From a list of elements segments, figure out if el belongs to it: */
    … …  
    117115        }
    118116        return 0;
    119 }
    120 
    121 
    122 /******************************************************************************************************************************
    123                                    RiftSegmentsFromSegments
    124 ******************************************************************************************************************************/
    125 
     117}/*}}}*/
     118/*FUNCTION RiftSegmentsFromSegments{{{*/
    126119void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
    127120       
    … …  
    189182        *priftsegments=riftsegments;
    190183        *pnriftsegs=nriftsegs;
    191 }
    192 
    193 /******************************************************************************************************************************
    194                                    DetermineGridElementListOnOneSideOfRift
    195 ******************************************************************************************************************************/
    196 
     184}/*}}}*/
     185/*FUNCTION DetermineGridElementListOnOneSideOfRift{{{*/
    197186int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel){
    198187
    … …  
    257246        *pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
    258247        return noerr;
    259 }
    260 
    261 /******************************************************************************************************************************
    262                                    UpdateSegments
    263 ******************************************************************************************************************************/
    264 
    265 int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs){
     248}/*}}}*/
     249/*FUNCTION UpdateSegments{{{*/
     250int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){
    266251
    267252        int noerr=1;
    … …  
    284269        /*First, update the existing segments to the new nodes :*/
    285270        for (i=0;i<nriftsegs;i++){
    286                 el1=*(riftsegments+4*i+0);
    287                 el2=*(riftsegments+4*i+1);
     271                el1=riftsegments[4*i+0];
     272                el2=riftsegments[4*i+1];
    288273                for (j=0;j<nsegs;j++){
    289                         if (*(segments+3*j+2)==(el1+1)){
     274                        if (segments[3*j+2]==(el1+1)){
    290275                                /*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
    291276                                 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts,
    … …  
    293278                                for (k=0;k<3;k++){
    294279                                        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])){
    295                                                 *(segments+3*j+0)=*(index+el1*3+k);
     280                                                *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods);
    296281                                                break;
    297282                                        }
    … …  
    299284                                for (k=0;k<3;k++){
    300285                                        if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+1)-1])  && (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
    301                                                 *(segments+3*j+1)=*(index+el1*3+k);
     286                                                *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods);
    302287                                                break;
    303288                                        }
    … …  
    308293                                for (k=0;k<3;k++){
    309294                                        if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+0)-1]) && (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
    310                                                 *(segments+3*(nsegs+i)+0)=*(index+el2*3+k);
     295                                                *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods);
    311296                                                break;
    312297                                        }
    … …  
    314299                                for (k=0;k<3;k++){
    315300                                        if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+1)-1]) && (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
    316                                                 *(segments+3*(nsegs+i)+1)=*(index+el2*3+k);
     301                                                *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods);
    317302                                                break;
    318303                                        }
    … …  
    324309                                for (k=0;k<3;k++){
    325310                                        if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+0)-1]) && (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
    326                                                 *(segments+3*j+0)=*(index+el2*3+k);
     311                                                *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods);
    327312                                                break;
    328313                                        }
    … …  
    330315                                for (k=0;k<3;k++){
    331316                                        if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+1)-1]) && (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
    332                                                 *(segments+3*j+1)=*(index+el2*3+k);
     317                                                *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods);
    333318                                                break;
    334319                                        }
    … …  
    339324                                for (k=0;k<3;k++){
    340325                                        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])){
    341                                                 *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);
     326                                                *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods);
    342327                                                break;
    343328                                        }
    … …  
    345330                                for (k=0;k<3;k++){
    346331                                        if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+1)-1]) && (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
    347                                                 *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);
     332                                                *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods);
    348333                                                break;
    349334                                        }
    … …  
    360345       
    361346        return noerr;
    362 }
    363 
    364 /******************************************************************************************************************************
    365                                    pnpoly
    366 ******************************************************************************************************************************/
     347}/*}}}*/
     348/*FUNCTION pnpoly{{{*/
    367349int pnpoly(int npol, double *xp, double *yp, double x, double y) {
    368350        int i, j, c = 0;
    … …  
    374356        }
    375357        return c;
    376 }
    377 
    378 /******************************************************************************************************************************
    379                                    IsInPoly
    380 ******************************************************************************************************************************/
    381 //void IsInPoly(double* in,double* xc,double* yc,int numnodes,double* x,double* y,int nods){
    382 //
    383 //      int i;
    384 //      double x0,y0;
    385 //
    386 //      /*Go through all nodes of the mesh:*/
    387 //      for (i=0;i<nods;i++){
    388 //              if (in[i]){
    389 //                      /*this node already is inside one of the contours, continue*/
    390 //                      continue;
    391 //              }
    392 //              /*pick up node: */
    393 //              x0=x[i];
    394 //              y0=y[i];
    395 //              if (pnpoly(numnodes,xc,yc,x0,y0)){
    396 //                      in[i]=1;
    397 //              }
    398 //      }
    399 //}
    400 
    401 /******************************************************************************************************************************
    402                                    FindElement
    403 ******************************************************************************************************************************/
    404 
     358}/*}}}*/
     359/*FUNCTION FindElement{{{*/
    405360int FindElement(double A,double B,double* index,int nel){
    406361
    … …  
    414369        }
    415370        return el;
    416 }
    417 /******************************************************************************************************************************
    418                                    SplitRiftSegments
    419 ******************************************************************************************************************************/
    420 
    421 int SplitRiftSegments(double** psegments,double** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,double*** priftssegments,int numrifts){
     371}/*}}}*/
     372/*FUNCTION SplitRiftSegments{{{*/
     373int SplitRiftSegments(double** psegments,double** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,double*** priftssegments,int numrifts,int nods,int nel){
    422374
    423375        /*Using segment markers, wring out the rift segments from the segments. Rift markers are
    … …  
    461413        for (i=0;i<numsegs;i++){
    462414                if (segmentmarkerlist[i]==1){
    463                         *(new_segments+3*counter+0)=*(segments+3*i+0);
    464                         *(new_segments+3*counter+1)=*(segments+3*i+1);
    465                         *(new_segments+3*counter+2)=*(segments+3*i+2);
     415                        new_segments[3*counter+0]=segments[3*i+0];
     416                        new_segments[3*counter+1]=segments[3*i+1];
     417                        new_segments[3*counter+2]=segments[3*i+2];
    466418                        new_segmentmarkers[counter]=segmentmarkerlist[i];
    467419                        counter++;
    … …  
    484436                for (j=0;j<numsegs;j++){
    485437                        if (segmentmarkerlist[j]==(2+i)){
    486                                 *(riftsegment+3*counter+0)=*(segments+3*j+0);
    487                                 *(riftsegment+3*counter+1)=*(segments+3*j+1);
    488                                 *(riftsegment+3*counter+2)=*(segments+3*j+2);
     438                                riftsegment[3*counter+0]=segments[3*j+0];_assert_(riftsegment[3*counter+0]<nods);
     439                                riftsegment[3*counter+1]=segments[3*j+1];_assert_(riftsegment[3*counter+1]<nods);
     440                                riftsegment[3*counter+2]=segments[3*j+2];_assert_(riftsegment[3*counter+2]<nel);
    489441                                counter++;
    490442                        }
    … …  
    504456        *priftsnumsegs=riftsnumsegs;
    505457        return noerr;
    506 }
    507 
    508 /******************************************************************************************************************************
    509                                    PairRiftElements
    510 ******************************************************************************************************************************/
    511 
     458}/*}}}*/
     459/*FUNCTION PairRiftElements{{{*/
    512460int PairRiftElements(int** priftsnumpairs, double*** priftspairs,int numrifts,int* riftsnumsegments, double** riftssegments,double* x,double* y){
    513461
    … …  
    558506
    559507        return noerr;
    560 }
    561 
    562 
    563 /******************************************************************************************************************************
    564                                    RemoveRifts
    565 ******************************************************************************************************************************/
    566 
    567 double dabs(double x){
    568         if (x<0)x=-x;
    569         return x;
    570 }
     508}/*}}}*/
     509/*FUNCTION RemoveRifts{{{*/
    571510int RemoveRifts(double** pindex,double** px,double** py,int* pnods,double** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,double** rifts1segments,double** rifts1pairs,int nel){
    572511
    … …  
    615554                if (y[i]<ymin)ymin=y[i];
    616555        }
    617         xmin=xmin-dabs(xmin);
    618         ymin=ymin-dabs(ymin);
     556        xmin=xmin-fabs(xmin);
     557        ymin=ymin-fabs(ymin);
    619558
    620559        /*Initialize two arrays, one for nodes that are going to be merged, the other with corresponding nodes being merge into: */
    … …  
    751690
    752691        return noerr;
    753 }
    754 
    755 /******************************************************************************************************************************
    756                                    IsRiftPresent
    757 ******************************************************************************************************************************/
    758 
     692}/*}}}*/
     693/*FUNCTION IsRiftPresent{{{*/
    759694int IsRiftPresent(int* priftflag,int* pnumrifts, double* segmentmarkerlist,int nsegs){
    760695
    … …  
    783718
    784719        return noerr;
    785 }
    786 
    787 /******************************************************************************************************************************
    788                                    OrderRifts
    789 ******************************************************************************************************************************/
    790 
     720}/*}}}*/
     721/*FUNCTION OrderRifts{{{*/
    791722int OrderRifts(double** priftstips, double** riftssegments,double** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){
    792723       
    … …  
    944875        *priftstips=riftstips;
    945876        return noerr;
    946 }
    947 
    948 /******************************************************************************************************************************
    949                                    PenaltyPairs
    950 ******************************************************************************************************************************/
    951 
     877}/*}}}*/
     878/*FUNCTION PenaltyPairs{{{*/
    952879int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,double** riftssegments,
    953880                int* riftsnumsegs,double** riftspairs,double* riftstips,double* x,double* y){

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

    r12070 r12080  
    66#define _SHARED_TRIMESH_H
    77
    8 
    98#include <stdio.h>
    109#include <math.h>
    1110
    12 
    13 
    1411//#define REAL double //took  it out because it may conflict with stdlib.h defines. put back if necessary
    15 
    1612int AssociateSegmentToElement(double** psegments,int nseg, double* index,int nel);
    1713int OrderSegments(double** psegments,int nseg, double* index,int nel);
    18                
    1914int GridInsideHole(double* px0,double* py0,int n,double* x,double* y);
    2015int FindElement(double A,double B,double* index,int nel);
    21 
    2216int SplitMeshForRifts(int* pnel,double** pindex,int* pnods,double** px,double** py,int* pnsegs,double** psegments,double** psegmentmarkerlist);
    23 
    2417int IsGridOnRift(int* riftsegments, int nriftsegs, int node);
    2518int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel);
    … …  
    2821void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments);
    2922int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel);
    30 int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs);
     23int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
    3124int pnpoly(int npol, double *xp, double *yp, double x, double y);
    3225int FindElement(double A,double B,double* index,int nel);
    3326int RemoveRifts(double** pindex,double** px,double** py,int* pnods,double** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,double** rifts1segments,double** rifts1pairs,int nel);
    3427int IsRiftPresent(int* priftflag,int* pnumrifts, double* segmentmarkerlist,int nsegs);
    35 int SplitRiftSegments(double** psegments,double** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,double*** priftssegments,int numrifts);
     28int SplitRiftSegments(double** psegments,double** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,double*** priftssegments,int numrifts,int nods,int nels);
    3629int OrderRifts(double** priftstips, double** riftssegments,double** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels);
    3730int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,double**  riftssegments,
    3831                int* riftsnumsegments,double** riftspairs,double* riftstips,double* x,double* y);
    39 
    4032int RemoveCornersFromRifts(double** pindex,int* pnel,double** px,double** py,int* pnods, double* segments,double* segmentmarkers,int num_seg);
    4133int PairRiftElements(int** priftsnumpairs, double*** priftspairs,int numrifts,int* riftsnumsegments, double** riftssegments,double* x,double* y);
    4234
    43 
    4435#endif  /* _SHARED_TRIMESH_H */
Note: See TracChangeset for help on using the changeset viewer.