Changeset 13588


Ignore:
Timestamp:
10/10/12 18:56:32 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: creating x layer (local for now)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp

    r13586 r13588  
    1111        _printLine_("      index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.");
    1212}/*}}}*/
     13void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg,
     14                        int  *pnumrifts, int  **priftsnumsegments, double ***priftssegments,
     15                        int  **priftsnumpairs, double ***priftspairs, double **priftstips,
     16                        double ***priftspenaltypairs, int  **priftsnumpenaltypairs);
    1317WRAPPER(TriMeshProcessRifts){
    1418
    15         /*Matlab arrays: */
     19        /* returned quantities: */
     20        int      numrifts;
     21        int     *riftsnumsegments     = NULL;
     22        double **riftssegments        = NULL;
     23        int     *riftsnumpairs        = NULL;
     24        double **riftspairs           = NULL;
     25        double  *riftstips            = NULL;
     26        double **riftspenaltypairs    = NULL;
     27        int     *riftsnumpenaltypairs = NULL;
     28
     29        /* input: */
     30        int     nel,nods;
     31        double *index          = NULL;
     32        double *x              = NULL;
     33        double *y              = NULL;
     34        double *segments       = NULL;
     35        double *segmentmarkers = NULL;
     36        int     num_seg;
     37
     38        /*Boot module*/
     39        MODULEBOOT();
     40
     41        /*checks on arguments on the matlab side: */
     42        CHECKARGUMENTS(NLHS,NRHS,&TriMeshProcessRiftsUsage);
     43
     44        /*Fetch data */
     45        FetchData(&index,&nel,NULL,INDEXIN);
     46        FetchData(&x,&nods,XIN);
     47        FetchData(&y,NULL,YIN);
     48        FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
     49        FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
     50
     51        /*call x layer*/
     52        TriMeshProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,
     53                                &numrifts,&riftsnumsegments,&riftssegments,&riftsnumpairs,&riftspairs,&riftstips,&riftspenaltypairs,&riftsnumpenaltypairs);
     54
     55        /*Output : */
     56        WriteData(INDEXOUT,index,nel,3);
     57        WriteData(XOUT,x,nods,1);
     58        WriteData(YOUT,y,nods,1);
     59        WriteData(SEGMENTSOUT,segments,num_seg,3);
     60        WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
     61
    1662        mxArray *pmxa_array  = NULL;
    1763        mxArray *pmxa_array2 = NULL;
    1864        mxArray *pmxa_array3 = NULL;
    19         int i,j,k,counter;
    20        
    21         /* returned quantities: */
    22         int      out_numrifts;
    23         int     *out_riftsnumsegments     = NULL;
    24         double **out_riftssegments        = NULL;
    25         int     *out_riftsnumpairs        = NULL;
    26         double **out_riftspairs           = NULL;
    27         double  *out_riftstips            = NULL;
    28         double **out_riftspenaltypairs    = NULL;
    29         int     *out_riftsnumpenaltypairs = NULL;
    30 
    31         /*empty rifts structure: */
    32         double* pNaN=NULL;
    33         const   char*   fnames[10];
    34         mwSize     ndim=2;
    35         mwSize          dimensions[2] = {1,1};
    36         double* pair=NULL;
    37        
    38         /* input: */
    39         double *index_in  = NULL;
    40         int     nel;
    41         double *x_in  = NULL; //copy of matlab vector
    42         int     nods;
    43         double *y_in               = NULL; //copy of matlab vector
    44         double *segments_in        = NULL;
    45         double *segmentmarkers_in  = NULL;
    46 
    47         /* state: */
    48         double *state = NULL;
    49         int     num_seg;
    50 
    51         /*rifts: */
    52         int     riftflag;
    53         int     numrifts;
    54 
    55         /*Boot module*/
    56         MODULEBOOT();
    57 
    58         /*checks on arguments on the matlab side: */
    59         CHECKARGUMENTS(NLHS,NRHS,&TriMeshProcessRiftsUsage);
    60 
    61         /*Fetch data */
    62         FetchData(&index_in,&nel,NULL,INDEXIN);
    63         FetchData(&x_in,&nods,XIN);
    64         FetchData(&y_in,NULL,YIN);
    65         FetchData(&segments_in,&num_seg,NULL,SEGMENTSIN);
    66         FetchData(&segmentmarkers_in,NULL,SEGMENTMARKERSIN);
    67 
    68         /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
    69          *all the nodes of this element belong to the segments (tends to happen when there are corners: */
    70         RemoveCornersFromRifts(&index_in,&nel,&x_in,&y_in,&nods,segments_in,segmentmarkers_in,num_seg);
    71 
    72         /*Figure out if we have rifts, and how many: */
    73         IsRiftPresent(&riftflag,&numrifts,segmentmarkers_in,num_seg);
    74        
    75         if(riftflag){   
    76                 SplitMeshForRifts(&nel,&index_in,&nods,&x_in,&y_in,&num_seg,&segments_in,&segmentmarkers_in);
     65        const char *fnames[10];
     66        mwSize      ndim          = 2;
     67        mwSize      dimensions[2] = {1,1};
     68        double     *pair          = NULL;
     69
     70
     71        /*Create a structure rifts where if i is a rift number, we have the following fields rifts(i).segments and rifts(i).numsegs: */
     72        fnames[0] = "numsegs";
     73        fnames[1] = "segments";
     74        fnames[2] = "pairs";
     75        fnames[3] = "tips";
     76        fnames[4] = "penaltypairs";
     77        fnames[5] = "fill";
     78        fnames[6] = "friction";
     79        fnames[7] = "fraction";
     80        fnames[8] = "fractionincrement";
     81        fnames[9] = "state";
     82
     83        dimensions[0]=numrifts;
     84        pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
     85
     86        for(int i=0;i<numrifts;i++){
     87
     88                /*Segments: */
     89                WriteData(&pmxa_array3,riftssegments[i],riftsnumsegments[i],3);
     90                mxSetField(pmxa_array,i,"segments",pmxa_array3);
     91                mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)riftsnumsegments[i]));
     92
     93                /*Element pairs: */
     94                WriteData(&pmxa_array3,riftspairs[i],riftsnumpairs[i],2);
     95                mxSetField(pmxa_array,i,"pairs",pmxa_array3);
     96
     97                /*Tips: */
     98                pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
     99                mxSetM(pmxa_array2,1);
     100                pair=(double*)mxMalloc(2*sizeof(double));
     101                pair[0]=*(riftstips+2*i+0);
     102                pair[1]=*(riftstips+2*i+1);
     103                mxSetN(pmxa_array2,2);
     104                mxSetPr(pmxa_array2,pair);
     105                mxSetField(pmxa_array,i,"tips",pmxa_array2);
     106
     107                /*Penalty pairs: */
     108                WriteData(&pmxa_array3,riftspenaltypairs[i],riftsnumpenaltypairs[i],7);
     109                mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
     110
     111                /*Friction fraction, fractionincrement  and fill: */
     112                mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
     113                mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum)); //default is ice
     114                mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
     115                mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1));
     116
     117                /*State: */
     118                double *state=(double*)mxMalloc(riftsnumpenaltypairs[i]*sizeof(double));
     119                for(int j=0;j<riftsnumpenaltypairs[i];j++) state[j]=FreeEnum;
     120                pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
     121                mxSetM(pmxa_array2,1);
     122                mxSetN(pmxa_array2,riftsnumpenaltypairs[i]);
     123                mxSetPr(pmxa_array2,state);
     124                mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
     125
     126                mxSetField(pmxa_array,i,"state",pmxa_array3);
    77127        }
    78128
    79         /*Order segments so that their normals point outside the domain: */
    80         OrderSegments(&segments_in,num_seg, index_in,nel);
    81        
    82         if(riftflag){
    83                
    84                 /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
    85                  *segmentmarkerlist:*/
    86                 SplitRiftSegments(&segments_in,&segmentmarkers_in,&num_seg,&out_numrifts,&out_riftsnumsegments,&out_riftssegments,numrifts,nods,nel);
    87 
    88                 /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
    89                 PairRiftElements(&out_riftsnumpairs,&out_riftspairs,out_numrifts,out_riftsnumsegments,out_riftssegments,x_in,y_in);
    90                
    91                 /*Order rifts so that they start from one tip, go to the other tip, and back: */
    92                 OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in,nods,nel);
    93 
    94                 /*Create penalty pairs, used by Imp: */
    95                 PenaltyPairs(&out_riftspenaltypairs,&out_riftsnumpenaltypairs,numrifts,out_riftssegments,out_riftsnumsegments,out_riftspairs,out_riftstips,x_in,y_in);
    96         }
    97 
    98 
    99         /*Output : */
    100         WriteData(INDEXOUT,index_in,nel,3);
    101         WriteData(XOUT,x_in,nods,1);
    102         WriteData(YOUT,y_in,nods,1);
    103         WriteData(SEGMENTSOUT,segments_in,num_seg,3);
    104         WriteData(SEGMENTMARKERSOUT,segmentmarkers_in,num_seg,1);
    105 
    106         if(riftflag){
    107                 /*Create a structure rifts where if i is a rift number, we have the following fields rifts(i).segments and rifts(i).numsegs: */
    108                 fnames[0] = "numsegs";
    109                 fnames[1] = "segments";
    110                 fnames[2] = "pairs";
    111                 fnames[3] = "tips";
    112                 fnames[4] = "penaltypairs";
    113                 fnames[5] = "fill";
    114                 fnames[6] = "friction";
    115                 fnames[7] = "fraction";
    116                 fnames[8] = "fractionincrement";
    117                 fnames[9] = "state";
    118 
    119                 dimensions[0]=out_numrifts;
    120                 pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
    121                
    122                 for (i=0;i<out_numrifts;i++){
    123 
    124                         /*Segments: */
    125                         WriteData(&pmxa_array3,out_riftssegments[i],out_riftsnumsegments[i],3);
    126                         mxSetField(pmxa_array,i,"segments",pmxa_array3);
    127                         mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)out_riftsnumsegments[i]));
    128 
    129                         /*Element pairs: */
    130                         WriteData(&pmxa_array3,out_riftspairs[i],out_riftsnumpairs[i],2);
    131                         mxSetField(pmxa_array,i,"pairs",pmxa_array3);
    132 
    133                         /*Tips: */
    134                         pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
    135                         mxSetM(pmxa_array2,1);
    136                         pair=(double*)mxMalloc(2*sizeof(double));
    137                         pair[0]=*(out_riftstips+2*i+0);
    138                         pair[1]=*(out_riftstips+2*i+1);
    139                         mxSetN(pmxa_array2,2);
    140                         mxSetPr(pmxa_array2,pair);
    141                         mxSetField(pmxa_array,i,"tips",pmxa_array2);
    142 
    143                         /*Penalty pairs: */
    144                         WriteData(&pmxa_array3,out_riftspenaltypairs[i],out_riftsnumpenaltypairs[i],7);
    145                         mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
    146 
    147                         /*Friction fraction, fractionincrement  and fill: */
    148                         mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
    149                         mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum)); //default is ice
    150                         mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
    151                         mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1));
    152 
    153                         /*State: */
    154                         state=(double*)mxMalloc(out_riftsnumpenaltypairs[i]*sizeof(double));
    155                         for(j=0;j<out_riftsnumpenaltypairs[i];j++)state[j]=FreeEnum;
    156                         pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
    157                         mxSetM(pmxa_array2,1);
    158                         mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
    159                         mxSetPr(pmxa_array2,state);
    160                         mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
    161                        
    162                         mxSetField(pmxa_array,i,"state",pmxa_array3);
    163                 }
    164         }
    165         else{
    166                 /*output NaN :*/
    167                 pNaN=(double*)mxMalloc(sizeof(double));
    168                 *pNaN=NAN;
    169                 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
    170                 mxSetM(pmxa_array,1);
    171                 mxSetN(pmxa_array,1);
    172                 mxSetPr(pmxa_array,pNaN);
    173                
    174         }
    175129        plhs[5]=pmxa_array;
    176130
    … …  
    178132        MODULEEND();
    179133}
     134void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg,
     135                        int  *pnumrifts, int  **priftsnumsegments, double ***priftssegments,
     136                        int  **priftsnumpairs, double ***priftspairs, double **priftstips,
     137                        double ***priftspenaltypairs, int  **priftsnumpenaltypairs){
     138
     139        /*Output*/
     140        int      numrifts,numrifts0;
     141        int     *riftsnumsegments     = NULL;
     142        double **riftssegments        = NULL;
     143        int     *riftsnumpairs        = NULL;
     144        double **riftspairs           = NULL;
     145        double  *riftstips            = NULL;
     146        double **riftspenaltypairs    = NULL;
     147        int     *riftsnumpenaltypairs = NULL;
     148
     149        /*Recover initial mesh*/
     150        int     nel            = *pnel;
     151        double *index          = *pindex;
     152        double *x              = *px;
     153        double *y              = *py;
     154        int     nods           = *pnods;
     155        double *segments       = *psegments;
     156        double *segmentmarkers = *psegmentmarkers;
     157        int     num_seg        = *pnum_seg;
     158
     159        /*Intermediary*/
     160        int     riftflag;
     161
     162        /*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie:
     163         *all the nodes of this element belong to the segments (tends to happen when there are corners: */
     164        RemoveCornersFromRifts(&index,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
     165
     166        /*Figure out if we have rifts, and how many: */
     167        IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
     168
     169        if(!riftflag) _error_("No rift present in mesh");
     170
     171        /*Split mesh*/
     172        SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
     173
     174        /*Order segments so that their normals point outside the domain: */
     175        OrderSegments(&segments,num_seg, index,nel);
     176
     177        /*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the
     178         *segmentmarkerlist:*/
     179        SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
     180
     181        /*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
     182        PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
     183
     184        /*Order rifts so that they start from one tip, go to the other tip, and back: */
     185        OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
     186
     187        /*Create penalty pairs, used by Imp: */
     188        PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
     189
     190        /*Assign output pointers for mesh*/
     191        *pnel            = nel;
     192        *pindex          = index;
     193        *px              = x;
     194        *py              = y;
     195        *pnods           = nods;
     196        *psegments       = segments;
     197        *psegmentmarkers = segmentmarkers;
     198        *pnum_seg        = num_seg;
     199
     200        /*Assign output pointers for rifts*/
     201        *pnumrifts             = numrifts;
     202        *priftsnumsegments     = riftsnumsegments;
     203        *priftssegments        = riftssegments;
     204        *priftsnumpairs        = riftsnumpairs;
     205        *priftspairs           = riftspairs;
     206        *priftstips            = riftstips;
     207        *priftspenaltypairs    = riftspenaltypairs;
     208        *priftsnumpenaltypairs = riftsnumpenaltypairs;
     209}
Note: See TracChangeset for help on using the changeset viewer.