Ignore:
Timestamp:
02/04/13 08:01:04 (12 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 14308

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

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

    r13975 r14310  
    3636}/*}}}*/
    3737/*FUNCTION GridElementsList{{{*/
    38 int GridElementsList(int** pGridElements, int* pNumGridElements,int node,double * index,int nel){
     38int GridElementsList(int** pGridElements, int* pNumGridElements,int node,int* index,int nel){
    3939
    4040        /*From a node, recover all the elements that are connected to it: */
    … …  
    5858        for (i=0;i<nel;i++){
    5959                for (j=0;j<3;j++){
    60                         if ((int)*(index+3*i+j)==node){
     60                        if (index[3*i+j]==node){
    6161                                if (NumGridElements<=(current_size-1)){
    6262                                        GridElements[NumGridElements]=i;
    … …  
    9090}/*}}}*/
    9191/*FUNCTION IsNeighbor{{{*/
    92 int IsNeighbor(int el1,int el2,double* index){
     92int IsNeighbor(int el1,int el2,int* index){
    9393        /*From a triangulation held in index, figure out if elements 1 and 2 have two nodes in common: */
    9494        int i,j;
    … …  
    9696        for (i=0;i<3;i++){
    9797                for (j=0;j<3;j++){
    98                         if (*(index+3*el1+i)==*(index+3*el2+j))count++;
     98                        if (index[3*el1+i]==index[3*el2+j])count++;
    9999                }
    100100        }
    … …  
    118118}/*}}}*/
    119119/*FUNCTION RiftSegmentsFromSegments{{{*/
    120 void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
     120void RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel,int* index,int nsegs,int* segments){
    121121
    122122        int i,counter;
    … …  
    126126        int* riftsegments=NULL;
    127127        int* riftsegments_uncompressed=NULL;
    128         double element_nodes[3];
     128        int element_nodes[3];
    129129
    130130        /*Allocate segmentflags: */
    … …  
    143143                element_nodes[2]=*(index+3*el+2);
    144144
    145                 *(index+3*el+0)=-1;
    146                 *(index+3*el+1)=-1;
    147                 *(index+3*el+2)=-1;
     145                index[3*el+0]=-1;
     146                index[3*el+1]=-1;
     147                index[3*el+2]=-1;
    148148
    149149                el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel);
    150150
    151151                /*Restore index: */
    152                 *(index+3*el+0)=element_nodes[0];
    153                 *(index+3*el+1)=element_nodes[1];
    154                 *(index+3*el+2)=element_nodes[2];
     152                index[3*el+0]=element_nodes[0];
     153                index[3*el+1]=element_nodes[1];
     154                index[3*el+2]=element_nodes[2];
    155155
    156156                if (el2!=-1){
    157157                        /*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
    158                     *(riftsegments_uncompressed+5*i+0)=1;
    159                     *(riftsegments_uncompressed+5*i+1)=el;
    160                     *(riftsegments_uncompressed+5*i+2)=el2;
    161                     *(riftsegments_uncompressed+5*i+3)=(int)*(segments+3*i+0);
    162                     *(riftsegments_uncompressed+5*i+4)=(int)*(segments+3*i+1);
    163                         nriftsegs++;
     158                    riftsegments_uncompressed[5*i+0]=1;
     159                    riftsegments_uncompressed[5*i+1]=el;
     160                    riftsegments_uncompressed[5*i+2]=el2;
     161                    riftsegments_uncompressed[5*i+3]=segments[3*i+0];
     162                         riftsegments_uncompressed[5*i+4]=segments[3*i+1];
     163                         nriftsegs++;
    164164                }
    165165        }
    … …  
    169169        counter=0;
    170170        for (i=0;i<nsegs;i++){
    171                 if (*(riftsegments_uncompressed+5*i+0)){
    172                         *(riftsegments+counter*4+0)=*(riftsegments_uncompressed+5*i+1);
    173                         *(riftsegments+counter*4+1)=*(riftsegments_uncompressed+5*i+2);
    174                         *(riftsegments+counter*4+2)=*(riftsegments_uncompressed+5*i+3);
    175                         *(riftsegments+counter*4+3)=*(riftsegments_uncompressed+5*i+4);
     171                if (riftsegments_uncompressed[5*i+0]){
     172                        riftsegments[counter*4+0]=riftsegments_uncompressed[5*i+1];
     173                        riftsegments[counter*4+1]=riftsegments_uncompressed[5*i+2];
     174                        riftsegments[counter*4+2]=riftsegments_uncompressed[5*i+3];
     175                        riftsegments[counter*4+3]=riftsegments_uncompressed[5*i+4];
    176176                        counter++;
    177177                }
    178178        }
    179 
    180179        xDelete<int>(riftsegments_uncompressed);
    181180
    … …  
    185184}/*}}}*/
    186185/*FUNCTION DetermineGridElementListOnOneSideOfRift{{{*/
    187 int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,double* index,int nel){
     186int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int node,int* index,int nel){
    188187
    189188        int noerr=1;
    … …  
    249248}/*}}}*/
    250249/*FUNCTION UpdateSegments{{{*/
    251 int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){
     250int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel){
    252251
    253252        int noerr=1;
    … …  
    255254        int el1,el2;
    256255
    257         double *segments          = NULL;
    258         double *segmentmarkerlist = NULL;
    259         int     nsegs;
     256        int *segments          = NULL;
     257        int *segmentmarkerlist = NULL;
     258        int  nsegs;
    260259
    261260        /*Recover input: */
    … …  
    265264
    266265        /*Reallocate segments: */
    267         segments         =xReNew<double>(segments,         nsegs*3,(nsegs+nriftsegs)*3);
    268         segmentmarkerlist=xReNew<double>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
     266        segments         =xReNew<int>(segments,         nsegs*3,(nsegs+nriftsegs)*3);
     267        segmentmarkerlist=xReNew<int>(segmentmarkerlist,nsegs,(nsegs+nriftsegs));
    269268
    270269        /*First, update the existing segments to the new nodes :*/
    … …  
    278277                                 *we can only rely on the position (x,y) of the rift nodes to create a segment:*/
    279278                                for (k=0;k<3;k++){
    280                                         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])){
     279                                        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])){
    281280                                                *(segments+3*j+0)=*(index+el1*3+k); _assert_(segments[3*j+0]<nods+1);
    282281                                                break;
    … …  
    284283                                }
    285284                                for (k=0;k<3;k++){
    286                                         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])){
     285                                        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])){
    287286                                                *(segments+3*j+1)=*(index+el1*3+k); _assert_(segments[3*j+1]<nods+1);
    288287                                                break;
    … …  
    293292                                *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
    294293                                for (k=0;k<3;k++){
    295                                         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])){
     294                                        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])){
    296295                                                *(segments+3*(nsegs+i)+0)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+0]<nods+1);
    297296                                                break;
    … …  
    299298                                }
    300299                                for (k=0;k<3;k++){
    301                                         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])){
     300                                        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])){
    302301                                                *(segments+3*(nsegs+i)+1)=*(index+el2*3+k); _assert_(segments[3*(nsegs+i)+1]<nods+1);
    303302                                                break;
    … …  
    309308                                /*Let's update segments[j][:] using  element el2 and the corresponding rift segment: */
    310309                                for (k=0;k<3;k++){
    311                                         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                                        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])){
    312311                                                *(segments+3*j+0)=*(index+el2*3+k); _assert_(segments[3*j+0]<nods+1);
    313312                                                break;
    … …  
    315314                                }
    316315                                for (k=0;k<3;k++){
    317                                         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                                        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])){
    318317                                                *(segments+3*j+1)=*(index+el2*3+k);_assert_(segments[3*j+1]<nods+1);
    319318                                                break;
    … …  
    324323                                *(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
    325324                                for (k=0;k<3;k++){
    326                                         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])){
     325                                        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])){
    327326                                                *(segments+3*(nsegs+i)+0)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+0]<nods+1);
    328327                                                break;
    … …  
    330329                                }
    331330                                for (k=0;k<3;k++){
    332                                         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])){
     331                                        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])){
    333332                                                *(segments+3*(nsegs+i)+1)=*(index+el1*3+k);_assert_(segments[3*(nsegs+i)+1]<nods+1);
    334333                                                break;
    … …  
    348347}/*}}}*/
    349348/*FUNCTION FindElement{{{*/
    350 int FindElement(double A,double B,double* index,int nel){
    351 
    352         int n;
     349int FindElement(int A,int B,int* index,int nel){
     350
    353351        int el=-1;
    354         for (n=0;n<nel;n++){
    355                 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))){
     352        for (int n=0;n<nel;n++){
     353                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))){
    356354                        el=n;
    357355                        break;
    … …  
    361359}/*}}}*/
    362360/*FUNCTION SplitRiftSegments{{{*/
    363 int SplitRiftSegments(double** psegments,double** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,double*** priftssegments,int numrifts,int nods,int nel){
     361int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nel){
    364362
    365363        /*Using segment markers, wring out the rift segments from the segments. Rift markers are
    … …  
    370368
    371369        /*input: */
    372         double *segments          = NULL;
    373         double *segmentmarkerlist = NULL;
     370        int *segments          = NULL;
     371        int *segmentmarkerlist = NULL;
    374372        int numsegs;
    375373
    376374        /*output: */
    377         int      new_numsegs;
    378         int     *riftsnumsegs       = NULL;
    379         double **riftssegments      = NULL;
    380         double  *new_segments       = NULL;
    381         double  *new_segmentmarkers = NULL;
     375        int   new_numsegs;
     376        int  *riftsnumsegs       = NULL;
     377        int **riftssegments      = NULL;
     378        int  *new_segments       = NULL;
     379        int  *new_segmentmarkers = NULL;
    382380
    383381        /*intermediary: */
    384         double* riftsegment=NULL;
     382        int* riftsegment=NULL;
    385383
    386384        /*Recover input arguments: */
    387         segments=*psegments;
    388         numsegs=*pnumsegs;
    389         segmentmarkerlist=*psegmentmarkerlist;
     385        segments          = *psegments;
     386        numsegs           = *pnumsegs;
     387        segmentmarkerlist = *psegmentmarkerlist;
    390388
    391389        /*First, figure out  how many segments will be left in 'segments': */
    … …  
    396394        /*Allocate new segments: */
    397395        new_numsegs=counter;
    398         new_segments=xNew<double>(new_numsegs*3);
    399         new_segmentmarkers=xNew<double>(new_numsegs);
     396        new_segments=xNew<int>(new_numsegs*3);
     397        new_segmentmarkers=xNew<int>(new_numsegs);
    400398
    401399        /*Copy new segments info : */
    … …  
    413411        /*Now deal with rift segments: */
    414412        riftsnumsegs=xNew<int>(numrifts);
    415         riftssegments=xNew<double*>(numrifts);
     413        riftssegments=xNew<int*>(numrifts);
    416414        for (i=0;i<numrifts;i++){
    417415                /*Figure out how many segments for rift i: */
    … …  
    421419                }
    422420                riftsnumsegs[i]=counter;
    423                 riftsegment=xNew<double>(counter*3);
     421                riftsegment=xNew<int>(counter*3);
    424422                /*Copy new segments info :*/
    425423                counter=0;
    … …  
    436434
    437435        /*Free ressources: */
    438         xDelete<double>(segments);
     436        xDelete<int>(segments);
    439437
    440438        /*Assign output pointers: */
    … …  
    448446}/*}}}*/
    449447/*FUNCTION PairRiftElements{{{*/
    450 int PairRiftElements(int** priftsnumpairs, double*** priftspairs,int numrifts,int* riftsnumsegments, double** riftssegments,double* x,double* y){
     448int PairRiftElements(int** priftsnumpairs,int*** priftspairs,int numrifts,int* riftsnumsegments,int** riftssegments,double* x,double* y){
    451449
    452450        int noerr=1;
    … …  
    454452
    455453        /*output: */
    456         int     *riftsnumpairs = NULL;
    457         double **riftspairs    = NULL;
     454        int  *riftsnumpairs = NULL;
     455        int **riftspairs    = NULL;
    458456
    459457        /*intermediary :*/
    460         int     numsegs;
    461         double* segments=NULL;
    462         double* pairs=NULL;
    463         int     node1,node2,node3,node4;
     458        int  numsegs;
     459        int* segments=NULL;
     460        int* pairs=NULL;
     461        int  node1,node2,node3,node4;
    464462
    465463        riftsnumpairs=xNew<int>(numrifts);
    466         riftspairs=xNew<double*>(numrifts);
     464        riftspairs=xNew<int*>(numrifts);
    467465        for (i=0;i<numrifts;i++){
    468466                segments=riftssegments[i];
    469                 numsegs=riftsnumsegments[i];
     467                numsegs =riftsnumsegments[i];
    470468                riftsnumpairs[i]=numsegs;
    471                 pairs=xNew<double>(2*numsegs);
     469                pairs=xNew<int>(2*numsegs);
    472470                for (j=0;j<numsegs;j++){
    473                         *(pairs+2*j+0)=*(segments+3*j+2); //retrieve element to which this segment belongs.
    474                         node1=(int)*(segments+3*j+0)-1; node2=(int)*(segments+3*j+1)-1;
     471                        pairs[2*j+0]=segments[3*j+2]; //retrieve element to which this segment belongs.
     472                        node1=segments[3*j+0]-1; node2=segments[3*j+1]-1;
    475473                        /*Find element facing on other side of rift: */
    476474                        for (k=0;k<numsegs;k++){
    477475                                if (k==j)continue;
    478                                 node3=(int)*(segments+3*k+0)-1; node4=(int)*(segments+3*k+1)-1;
     476                                node3=segments[3*k+0]-1; node4=segments[3*k+1]-1;
    479477                                /*We are trying to find 2 elements, where position of node3 == position of node1, and position of node4 == position of node2*/
    480478                                if (   (x[node3]==x[node1]) && (y[node3]==y[node1]) && (x[node4]==x[node2]) && (y[node4]==y[node2])
    481479                                    || (x[node3]==x[node2]) && (y[node3]==y[node2]) && (x[node4]==x[node1]) && (y[node4]==y[node1])  ){
    482480                                        /*We found the corresponding element: */
    483                                         *(pairs+2*j+1)=*(segments+3*k+2);
     481                                        pairs[2*j+1]=segments[3*k+2];
    484482                                        break;
    485483                                }
    … …  
    495493}/*}}}*/
    496494/*FUNCTION RemoveRifts{{{*/
    497 int RemoveRifts(double** pindex,double** px,double** py,int* pnods,double** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,double** rifts1segments,double** rifts1pairs,int nel){
     495int RemoveRifts(int** pindex,double** px,double** py,int* pnods,int** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,int** rifts1segments,int** rifts1pairs,int nel){
    498496
    499497        int noerr=1;
    500498        int i,j,k,counter,counter1,counter2;
    501499
    502         /*input: */
    503         double* index=NULL;
    504         double* x=NULL;
    505         double* y=NULL;
    506         int     nods;
    507         double* segments=NULL;
    508         int     numsegs;
    509 
    510500        /*intermediary: */
    511         double *riftsegments = NULL;
    512         double *riftpairs    = NULL;
     501        int    *riftsegments = NULL;
     502        int    *riftpairs    = NULL;
    513503        int     node1,node2,node3,node4,temp_node;
    514         double  el2;
     504        int     el2;
    515505        int     newnods; //temporary # node counter.
    516506        double  xmin,ymin;
    … …  
    523513
    524514        /*Recover input: */
    525         index=*pindex;
    526         x=*px;
    527         y=*py;
    528         nods=*pnods;;
    529         segments=*psegments;
    530         numsegs=*pnumsegs;
     515        int    *index    = *pindex;
     516        double *x        = *px;
     517        double *y        = *py;
     518        int     nods     = *pnods;     ;
     519        int    *segments = *psegments;
     520        int     numsegs  = *pnumsegs;
    531521
    532522        /*initialize newnods : */
    … …  
    558548                riftpairs=rifts1pairs[i];
    559549                for (j=0;j<rifts1numsegs[i];j++){
    560                         el2=*(riftpairs+2*j+1);
     550                        el2=riftpairs[2*j+1];
    561551                        node1=(int)*(riftsegments+3*j+0);
    562552                        node2=(int)*(riftsegments+3*j+1);
    … …  
    565555                        for (k=0;k<rifts1numsegs[i];k++){
    566556                                if (*(riftsegments+3*k+2)==el2){
    567                                         node3=(int)*(riftsegments+3*k+0);
    568                                         node4=(int)*(riftsegments+3*k+1);
     557                                        node3=*(riftsegments+3*k+0);
     558                                        node4=*(riftsegments+3*k+1);
    569559                                        break;
    570560                                }
    … …  
    677667}/*}}}*/
    678668/*FUNCTION IsRiftPresent{{{*/
    679 int IsRiftPresent(int* priftflag,int* pnumrifts, double* segmentmarkerlist,int nsegs){
     669int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){
    680670
    681671        int i;
    … …  
    686676        int numrifts=0;
    687677
    688         double maxmark=1; //default marker for regular segments
     678        int maxmark=1; //default marker for regular segments
    689679
    690680        /*Any marker >=2 indicates a certain rift: */
    … …  
    696686                }
    697687        }
    698         if (numrifts)riftflag=1;
     688        if(numrifts)riftflag=1;
    699689
    700690        /*Assign output pointers:*/
    … …  
    704694}/*}}}*/
    705695/*FUNCTION OrderRifts{{{*/
    706 int OrderRifts(double** priftstips, double** riftssegments,double** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){
     696int OrderRifts(int** priftstips,int** riftssegments,int** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y,int nods,int nels){
    707697
    708698        int noerr=1;
    … …  
    710700
    711701        /*intermediary: */
    712         double *riftsegments = NULL;
    713         double *riftpairs    = NULL;
     702        int *riftsegments = NULL;
     703        int *riftpairs    = NULL;
    714704        int numsegs;
    715705
    716706        /*ordering and copy: */
    717         int    *order             = NULL;
    718         double *riftsegments_copy = NULL;
    719         double *riftpairs_copy    = NULL;
     707        int *order             = NULL;
     708        int *riftsegments_copy = NULL;
     709        int *riftpairs_copy    = NULL;
    720710
    721711        /*node and element manipulation: */
    722         int     node1,node2,node3,node4,temp_node,tip1,tip2,node;
    723         double el2;
    724         int     already_ordered=0;
     712        int node1,node2,node3,node4,temp_node,tip1,tip2,node;
     713        int el2;
     714        int already_ordered=0;
    725715
    726716        /*output: */
    727         double* riftstips=NULL;
     717        int* riftstips=NULL;
    728718
    729719        /*Allocate byproduct of this routine, riftstips: */
    730         riftstips=xNew<double>(numrifts*2);
     720        riftstips=xNew<int>(numrifts*2);
    731721
    732722        /*Go through all rifts: */
    733723        for (i=0;i<numrifts;i++){
    734                 riftsegments=riftssegments[i];
    735                 riftpairs=riftspairs[i];
    736                 numsegs=riftsnumsegments[i];
     724                riftsegments = riftssegments[i];
     725                riftpairs    = riftspairs[i];
     726                numsegs      = riftsnumsegments[i];
    737727
    738728                /*Allocate copy of riftsegments and riftpairs,
    739729                 *as well as ordering vector: */
    740                 riftsegments_copy=xNew<double>(numsegs*3);
    741                 riftpairs_copy=xNew<double>(numsegs*2);
     730                riftsegments_copy=xNew<int>(numsegs*3);
     731                riftpairs_copy=xNew<int>(numsegs*2);
    742732                order=xNew<int>(numsegs);
    743733
    … …  
    748738                for (j=0;j<numsegs;j++){
    749739                        el2=*(riftpairs+2*j+1);
    750                         node1=(int)*(riftsegments+3*j+0);
    751                         node2=(int)*(riftsegments+3*j+1);
     740                        node1=*(riftsegments+3*j+0);
     741                        node2=*(riftsegments+3*j+1);
    752742                        /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and
    753743                         *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
    754744                        for (k=0;k<numsegs;k++){
    755745                                if (*(riftsegments+3*k+2)==el2){
    756                                         node3=(int)*(riftsegments+3*k+0);
    757                                         node4=(int)*(riftsegments+3*k+1);
     746                                        node3=*(riftsegments+3*k+0);
     747                                        node4=*(riftsegments+3*k+1);
    758748                                        break;
    759749                                }
    … …  
    796786
    797787                /*Record tips in riftstips: */
    798                 *(riftstips+2*i+0)=(double)tip1;
    799                 *(riftstips+2*i+1)=(double)tip2;
     788                *(riftstips+2*i+0)=tip1;
     789                *(riftstips+2*i+1)=tip2;
    800790
    801791                /*We have the two tips for this rift.  Go from tip1 to tip2, and figure out the order in which segments are sequential.
    … …  
    804794                for (counter=0;counter<numsegs;counter++){
    805795                        for (j=0;j<numsegs;j++){
    806                                 node1=(int)*(riftsegments+3*j+0);
    807                                 node2=(int)*(riftsegments+3*j+1);
     796                                node1=*(riftsegments+3*j+0);
     797                                node2=*(riftsegments+3*j+1);
    808798
    809799                                if ((node1==node) || (node2==node)){
    … …  
    849839
    850840                xDelete<int>(order);
    851                 xDelete<double>(riftsegments_copy);
    852                 xDelete<double>(riftpairs_copy);
     841                xDelete<int>(riftsegments_copy);
     842                xDelete<int>(riftpairs_copy);
    853843
    854844        }
    … …  
    859849}/*}}}*/
    860850/*FUNCTION PenaltyPairs{{{*/
    861 int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,double** riftssegments,
    862                 int* riftsnumsegs,double** riftspairs,double* riftstips,double* x,double* y){
     851int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,int** riftssegments,
     852                int* riftsnumsegs,int** riftspairs,int* riftstips,double* x,double* y){
    863853
    864854        int noerr=1;
    … …  
    875865        /*intermediary: */
    876866        int numsegs;
    877         double* riftsegments=NULL;
    878         double* riftpairs=NULL;
     867        int* riftsegments=NULL;
     868        int* riftpairs=NULL;
    879869        int counter;
    880870        double normal[2];
    … …  
    980970                /*Renormalize normals: */
    981971                for(j=0;j<counter;j++){
    982                         double magnitude=sqrt(pow( *(riftpenaltypairs+j*7+4),2) + pow( *(riftpenaltypairs+j*7+5),2) );
     972                        double magnitude=sqrt(pow( double(riftpenaltypairs[j*7+4]),2) + pow( double(riftpenaltypairs[j*7+5]),2) );
    983973                        *(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
    984974                        *(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
    … …  
    993983        *priftsnumpenaltypairs=riftsnumpenaltypairs;
    994984        return noerr;
    995 }
    996 
    997 /******************************************************************************************************************************
    998                                    RemoveCorners
    999 ******************************************************************************************************************************/
    1000 
    1001 int RemoveCornersFromRifts(double** pindex,int* pnel,double** px,double** py,int* pnods, double* segments,double* segmentmarkers,int num_seg){
     985}/*}}}*/
     986int RemoveCornersFromRifts(int** pindex,int* pnel,double** px,double** py,int* pnods,int* segments,int* segmentmarkers,int num_seg){/*{{{*/
    1002987
    1003988        int noerr=1;
    1004989        int i,j,k;
    1005         double node1,node2,node3;
     990        int node1,node2,node3;
    1006991        int el;
    1007 
    1008         /*input: */
    1009         double* index=NULL;
    1010         int     nel;
    1011         double* x=NULL;
    1012         double* y=NULL;
    1013         int     nods;
    1014992        double  pair[2];
    1015993        int     pair_count=0;
    … …  
    1017995
    1018996        /*Recover input: */
    1019         index=*pindex;
    1020         nel=*pnel;
    1021         x=*px;
    1022         y=*py;
    1023         nods=*pnods;
     997        int    *index = *pindex;
     998        int     nel   = *pnel;
     999        double *x     = *px;
     1000        double *y     = *py;
     1001        int     nods  = *pnods;
    10241002
    10251003        for (i=0;i<num_seg;i++){
    … …  
    10831061                                        x[nods]=(x[(int)node1-1]+x[(int)node2-1]+x[(int)node3-1])/3;
    10841062                                        y[nods]=(y[(int)node1-1]+y[(int)node2-1]+y[(int)node3-1])/3;
    1085                                         index=xReNew<double>(index,nel*3,(nel+2*3));
     1063                                        index=xReNew<int>(index,nel*3,(nel+2*3));
    10861064                                        /*First, reassign element el: */
    10871065                                        *(index+3*el+0)=node1;
    … …  
    10981076                                        /*we need  to change the segment elements corresponding to el: */
    10991077                                        for (k=0;k<num_seg;k++){
    1100                                                 if (*(segments+3*k+2)==(double)(el+1)){
    1101                                                         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);
    1102                                                         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);
    1103                                                         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);
     1078                                                if (*(segments+3*k+2)==(el+1)){
     1079                                                        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;
     1080                                                        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;
     1081                                                        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;
    11041082                                                }
    11051083                                        }
    … …  
    11211099        *pnods=nods;
    11221100        return noerr;
    1123 }
     1101}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.