Ignore:
Timestamp:
04/25/14 09:27:03 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: iomodel->faces now has all its vertices

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp

    r17700 r17842  
    126126        /*Intermediaries*/
    127127        bool exist;
    128         int  i,j,k,v0,cols;
    129         int  maxnbf,nbf,elementnbf,elementnbv,facemaxnbv;
     128        int  i,j,k,v0,cols,facescols;
     129        int  maxnbf,nbf,elementnbf,elementnbv,facenbv,facemaxnbv;
    130130        int *elementfaces         = NULL;
    131131        int *elementfaces_markers = NULL;
    132132
    133         /*Maximum number of faces*/
    134         maxnbf = 5*iomodel->numberofelements;
     133        /*Mesh specific face indexing per element*/
     134        switch(iomodel->meshelementtype){
     135                case PentaEnum:
     136                        elementnbv = 6; /*Number of vertices per element*/
     137                        elementnbf = 5; /*Number of faces per element*/
     138                        facemaxnbv = 4; /*Maximum number of vertices per face*/
     139                        cols       = facemaxnbv + 1;
     140                        elementfaces         = xNew<int>(elementnbf*cols);
     141                        elementfaces_markers = xNew<int>(elementnbf);
     142                        /*2 triangles*/
     143                        elementfaces_markers[0] = 1;
     144                        elementfaces_markers[1] = 1;
     145                        elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0;  elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
     146                        elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3;  elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
     147                        /*3 quads*/
     148                        elementfaces_markers[2] = 2;
     149                        elementfaces_markers[3] = 2;
     150                        elementfaces_markers[4] = 2;
     151                        elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 1;  elementfaces[cols*2+2] = 2; elementfaces[cols*2+3] = 5;  elementfaces[cols*2+4] = 4;
     152                        elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 2;  elementfaces[cols*3+2] = 0; elementfaces[cols*3+3] = 3;  elementfaces[cols*3+4] = 5;
     153                        elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 0;  elementfaces[cols*4+2] = 1; elementfaces[cols*4+3] = 4;  elementfaces[cols*4+4] = 3;
     154                        break;
     155                default:
     156                _error_("mesh "<< EnumToStringx(iomodel->meshelementtype) <<" not supported");
     157        }
     158
     159        /*Allocate connectivity*/
     160        int *element_face_connectivity = xNew<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
     161
     162        /*Maximum number of faces for initial allocation*/
     163        maxnbf     = elementnbf*iomodel->numberofelements;
     164        facescols  = 4+facemaxnbv; _assert_(facescols>6);
    135165
    136166        /*Initialize intermediaries*/
    137         int*  facestemp = xNew<int>(maxnbf*6);         /*format: [vertex1 vertex2 vertex3 element1 element2 marker]    */
    138         for(i=0;i<maxnbf;i++) facestemp[i*6+4]=-1;     /*Initialize last column of faces as -1 (boundary face)         */
     167        int* facestemp = xNew<int>(maxnbf*facescols);        /*format: [element1 element2 marker nbv vertex1 vertex2 vertex3 ...]    */
     168        for(i=0;i<maxnbf;i++) facestemp[i*facescols+1]=-1;   /*Initialize second column of faces as -1 (boundary face)               */
    139169
    140170        /*Initialize chain*/
    … …  
    143173        for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
    144174
    145         /*Mesh specific face indexing per element*/
    146         if(iomodel->meshelementtype==PentaEnum){
    147                 elementnbv = 6; /*Number of vertices per element*/
    148                 elementnbf = 5; /*Number of faces per element*/
    149                 facemaxnbv = 4; /*Maximum number of vertices per face*/
    150                 cols       = facemaxnbv + 1;
    151                 elementfaces         = xNew<int>(elementnbf*cols);
    152                 elementfaces_markers = xNew<int>(elementnbf);
    153                 /*2 triangles*/
    154                 elementfaces_markers[0] = 1;
    155                 elementfaces_markers[1] = 1;
    156                 elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0;  elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
    157                 elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3;  elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
    158                 /*3 quads*/
    159                 elementfaces_markers[2] = 2;
    160                 elementfaces_markers[3] = 2;
    161                 elementfaces_markers[4] = 2;
    162                 elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 1;  elementfaces[cols*2+2] = 2; elementfaces[cols*2+3] = 5;  elementfaces[cols*2+4] = 4;
    163                 elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 2;  elementfaces[cols*3+2] = 0; elementfaces[cols*3+3] = 3;  elementfaces[cols*3+4] = 5;
    164                 elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 0;  elementfaces[cols*4+2] = 1; elementfaces[cols*4+3] = 4;  elementfaces[cols*4+4] = 3;
    165         }
    166         else{
    167                 _error_("mesh dimension not supported yet");
    168         }
    169         int *element_face_connectivity = xNew<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
    170 
    171175        /*Initialize number of faces and list of vertex indices*/
    172176        nbf = 0;
    … …  
    177181
    178182                        /*Get indices of current face*/
    179                         for(k=0;k<elementfaces[cols*j+0];k++){
     183                        facenbv = elementfaces[cols*j+0];
     184                        for(k=0;k<facenbv;k++){
    180185                                v[k] = iomodel->elements[i*elementnbv + elementfaces[cols*j+k+1]] - 1;
    181186                        }
    … …  
    183188                        /*Sort list of vertices*/
    184189                        HeapSort(v,elementfaces[cols*j+0]);
    185                         v0 = v[0];
     190                        v0 = v[0]; _assert_(v0>=0 && v0<iomodel->numberofvertices);
    186191
    187192                        /*This face a priori has not been processed yet*/
    188193                        exist = false;
    189194
    190                         /*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/
    191                         _assert_(v0>=0 && v0<iomodel->numberofvertices);
    192                         for(int e=head_minv[v0]; e!=-1; e=next_face[e]){
    193                                 if(facestemp[e*6+1]==v[1]+1 && facestemp[e*6+2]==v[2]+1){
     195                        /*Go through all processed faces connected to v0 and check whether we have seen this face yet*/
     196                        for(int f=head_minv[v0]; f!=-1; f=next_face[f]){
     197                                if(facestemp[f*facescols+5]==v[1]+1 && facestemp[f*facescols+6]==v[2]+1){
    194198                                        exist = true;
    195                                         facestemp[e*6+4]=i+1;
    196                                         element_face_connectivity[i*elementnbf+j]=e;
     199                                        facestemp[f*facescols+1]=i+1;
     200                                        element_face_connectivity[i*elementnbf+j]=f;
    197201                                        break;
    198202                                }
    … …  
    204208
    205209                                /*Update faces*/
    206                                 facestemp[nbf*6+0] = v[0]+1;
    207                                 facestemp[nbf*6+1] = v[1]+1;
    208                                 facestemp[nbf*6+2] = v[2]+1;
    209                                 facestemp[nbf*6+3] = i+1;
    210                                 facestemp[nbf*6+5] = elementfaces_markers[j];
     210                                facestemp[nbf*facescols+0] = i+1;
     211                                facestemp[nbf*facescols+2] = elementfaces_markers[j];
     212                                facestemp[nbf*facescols+3] = facenbv;
     213                                for(k=0;k<facenbv;k++) facestemp[nbf*facescols+4+k] = v[k]+1;
    211214
    212215                                /*Update Connectivity*/
    … …  
    231234
    232235        /*Create final faces (now that we have the correct size)*/
    233         int* faces = xNew<int>(nbf*6); /*vertex1 vertex2 vertex3 element1 element2 marker*/
    234         xMemCpy<int>(faces,facestemp,nbf*6);
     236        int* faces = xNew<int>(nbf*facescols);
     237        xMemCpy<int>(faces,facestemp,nbf*facescols);
    235238        xDelete<int>(facestemp);
    236239
    … …  
    238241        iomodel->faces                     = faces;
    239242        iomodel->numberoffaces             = nbf;
     243        iomodel->facescols                 = facescols;
    240244        iomodel->elementtofaceconnectivity = element_face_connectivity;
    241245}/*}}}*/
    242 void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber){/*{{{*/
    243 
    244         int numbervertices;
    245         if(iomodel->domaintype==Domain3DEnum){
    246                 if((iomodel->faces[6*facenumber+5])==1){
    247                         numbervertices=3;
    248                 }
    249                 else if((iomodel->faces[6*facenumber+5])==2){
    250                         numbervertices=4;
    251                 }
    252                 else _error_("face marker not supported yet ");
    253         }
    254         else _error_("mesh type not supported yet");
    255 
    256         int* facevertices = xNew<int>(numbervertices);
    257         if(numbervertices==3){
    258                 for(int i=0;i<3;i++) facevertices[i]=iomodel->faces[6*facenumber+i]-1;
    259         }
    260         else if(numbervertices==4){
    261                 int  i,j,k,cols,faceid;
    262                 int  maxnbf,nbf,elementnbf,elementnbv,facemaxnbv;
    263                 int *elementfaces         = NULL;
    264                 int *elementfaces_markers = NULL;
    265                 int elementid=iomodel->faces[6*facenumber+3];
    266                 int counter=0;
    267 
    268                 /*Recreate element properties*/
    269                 elementnbv = 6; /*Number of vertices per element*/
    270                 elementnbf = 5; /*Number of faces per element*/
    271                 facemaxnbv = 4; /*Maximum number of vertices per face*/
    272                 cols       = facemaxnbv + 1;
    273                 elementfaces         = xNew<int>(elementnbf*cols);
    274                 elementfaces_markers = xNew<int>(elementnbf);
    275                 /*2 triangles*/
    276                 elementfaces_markers[0] = 1;
    277                 elementfaces_markers[1] = 1;
    278                 elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0;  elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
    279                 elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3;  elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
    280                 /*3 quads*/
    281                 elementfaces_markers[2] = 2;
    282                 elementfaces_markers[3] = 2;
    283                 elementfaces_markers[4] = 2;
    284                 elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 1;  elementfaces[cols*2+2] = 2; elementfaces[cols*2+3] = 5;  elementfaces[cols*2+4] = 4;
    285                 elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 2;  elementfaces[cols*3+2] = 0; elementfaces[cols*3+3] = 3;  elementfaces[cols*3+4] = 5;
    286                 elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 0;  elementfaces[cols*4+2] = 1; elementfaces[cols*4+3] = 4;  elementfaces[cols*4+4] = 3;
    287 
    288                 for(faceid=2;faceid<5;faceid++){
    289                         counter=0;
    290                         for(j=0;j<3;j++){
    291                                 for(k=1;k<5;k++){
    292                                         if(iomodel->elements[(elementid-1)*6+elementfaces[cols*faceid+k]] == iomodel->faces[6*facenumber+j]) counter++;
    293                                 }
    294                         }
    295                         if(counter==3) break;
    296                 }
    297                 if(counter!=3) _error_("face not found in element");
    298 
    299                 for(j=0;j<4;j++){
    300                         facevertices[j]=iomodel->elements[(elementid-1)*6+elementfaces[cols*faceid+j+1]];
    301                 }
    302 
    303                 /*Delete*/
    304                 xDelete<int>(elementfaces);
    305                 xDelete<int>(elementfaces_markers);
    306         }
    307 
    308         /*Output results*/
    309         *pverticesid=facevertices;
    310         *pnumvertices=numbervertices;
    311 }/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.