Changeset 17842


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

NEW: iomodel->faces now has all its vertices

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r17692 r17842  
    3737        this->numberoffaces=-1;
    3838        this->numberofedges=-1;
     39        this->facescols=-1;
    3940        this->elements=NULL;
    4041        this->faces=NULL;
    … …  
    8384        FetchData(&this->numberofelements,MeshNumberofelementsEnum);
    8485        FetchData(&this->elements,NULL,NULL,MeshElementsEnum);
     86        this->facescols                       = -1;
    8587        this->faces                           = NULL;
    8688        this->edges                           = NULL;

  • issm/trunk-jpl/src/c/classes/IoModel.h

    r17692 r17842  
    3838                int  numberoffaces;
    3939                int  numberofedges;
     40                int  facescols;
    4041                int *elements;
    4142                int *faces;

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

    r17692 r17842  
    3535        bool        isnan;
    3636        int         i,j,count,elementnbv,numfacevertices;
    37         int*        faceverticesid   = NULL;
    3837        IssmDouble  value;
    3938        IssmDouble *times            = NULL;
    … …  
    119118                                if(iomodel->meshelementtype==PentaEnum){
    120119                                        for(i=0;i<iomodel->numberoffaces;i++){
    121                                                 if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     120                                                if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    122121                                                        if(my_faces[i]){
    123                                                                 FaceGetVertexIndices(iomodel,&numfacevertices,&faceverticesid,i);
    124                                                                 isnan=0;
     122                                                                numfacevertices = iomodel->faces[i*iomodel->facescols+3];
     123                                                                value=0.;
    125124                                                                for(j=0;j<numfacevertices;j++){
    126                                                                         if(xIsNan<IssmDouble>(spcdata[faceverticesid[j]-1])) isnan=1;
     125                                                                        value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1];
    127126                                                                }
    128                                                                 if(isnan==0){
    129                                                                         value=0;
    130                                                                         for(j=0;j<numfacevertices;j++){
    131                                                                                 value=value+spcdata[faceverticesid[j]-1]/numfacevertices;
    132                                                                         }
    133                                                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,
     127                                                                value = value/reCast<IssmDouble>(numfacevertices);
     128                                                                if(!xIsNan<IssmDouble>(value)){
     129                                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,
     130                                                                                                        iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,
    134131                                                                                                        dof,value,analysis_type));
    135132                                                                        count++;
    136133                                                                }
    137                                                                 xDelete<int>(faceverticesid);
    138134                                                        }
    139135                                                }
    … …  
    179175                                }
    180176                                for(i=0;i<iomodel->numberoffaces;i++){
    181                                         if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     177                                        if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    182178                                                if(my_faces[i]){
    183                                                         FaceGetVertexIndices(iomodel,&numfacevertices,&faceverticesid,i);
    184                                                         isnan=0;
     179                                                        numfacevertices = iomodel->faces[i*iomodel->facescols+3];
     180                                                        value=0.;
    185181                                                        for(j=0;j<numfacevertices;j++){
    186                                                                 if(xIsNan<IssmDouble>(spcdata[faceverticesid[j]-1])) isnan=1;
    187                                                         }
    188                                                         if(isnan==0){
    189                                                                 value=0;
    190                                                                 for(j=0;j<numfacevertices;j++){
    191                                                                         value=value+spcdata[faceverticesid[j]-1]/numfacevertices;
    192                                                                 }
    193                                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1,
     182                                                                value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1];
     183                                                        }
     184                                                        value = value/reCast<IssmDouble>(numfacevertices);
     185                                                        if(!xIsNan<IssmDouble>(value)){
     186                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,
     187                                                                                                iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1,
    194188                                                                                                dof,value,analysis_type));
    195                                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2,
     189                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,
     190                                                                                                iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2,
    196191                                                                                                dof,value,analysis_type));
    197                                                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3,
     192                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,
     193                                                                                                iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3,
    198194                                                                                                dof,value,analysis_type));
    199195                                                                count=count+3;
    200196                                                        }
    201                                                         xDelete<int>(faceverticesid);
    202197                                                }
    203198                                        }

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

    r17692 r17842  
    55#include "../../classes/classes.h"
    66#include "../../shared/shared.h"
     7#include "./ModelProcessorx.h"
    78
    8 void CreateEdges(IoModel* iomodel){
     9void CreateEdges(IoModel* iomodel){/*{{{*/
    910
    1011        /*If edges are already present, exit*/
    … …  
    142143        iomodel->elementtoedgeconnectivity = element_edge_connectivity;
    143144        iomodel->numberofedges             = nbe;
     145}/*}}}*/
     146void EdgeOnBoundaryFlags(IoModel* iomodel,bool** pflags){
     147
     148        /*Get faces and edges*/
     149        if(!iomodel->edges) CreateEdges(iomodel);
     150        if(!iomodel->faces) CreateFaces(iomodel);
     151
     152        /*Allocate output*/
     153        bool* flags = xNew<bool>(iomodel->numberofedges);
     154
     155        /*Now, loop over the faces, whenever it is not connected to a second element, all edges are on boundary*/
     156        for(int i=0;i<iomodel->numberoffaces;i++){
     157
     158                if(iomodel->faces[i*iomodel->facescols+4]==-1){
     159
     160                        /*Get all edges for this face*/
     161                }
     162        }
     163
     164        /*Clean up and return*/
     165        *pflags = flags;
    144166}

  • 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 }/*}}}*/

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

    r17795 r17842  
    147147                                FacesPartitioning(&my_faces,iomodel);
    148148                                for(i=0;i<iomodel->numberoffaces;i++){
    149                                         if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     149                                        if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    150150                                                if(my_faces[i]){
    151151                                                        node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,approximation);
    … …  
    153153                                                }
    154154                                        }
    155                                         else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
     155                                        else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    156156                                                /*Nothing*/
    157157                                        }
    … …  
    196196                        id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges;
    197197                        for(i=0;i<iomodel->numberoffaces;i++){
    198                                 if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     198                                if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    199199                                        if(my_faces[i]){
    200200                                                node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
    … …  
    207207                                        counter=counter+3;
    208208                                }
    209                                 else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
     209                                else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    210210                                        /*Nothing*/
    211211                                }
    … …  
    312312                                FacesPartitioning(&my_faces,iomodel);
    313313                                for(i=0;i<iomodel->numberoffaces;i++){
    314                                         if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     314                                        if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    315315                                                if(my_faces[i]){
    316316                                                        node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
    … …  
    318318                                                }
    319319                                        }
    320                                         else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
     320                                        else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    321321                                                /*Nothing*/
    322322                                        }
    … …  
    377377                        id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges;
    378378                        for(i=0;i<iomodel->numberoffaces;i++){
    379                                 if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
     379                                if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    380380                                        if(my_faces[i]){
    381381                                                node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,FSvelocityEnum);
    … …  
    388388                                        counter=counter+3;
    389389                                }
    390                                 else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
     390                                else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    391391                                        /*Nothing*/
    392392                                }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r17383 r17842  
    3333void CreateFaces(IoModel* iomodel);
    3434void CreateFaces3d(IoModel* iomodel);
    35 void FaceGetVertexIndices(IoModel* iomodel,int* pnumvertices,int** pverticesid,int facenumber);
     35void EdgeOnBoundaryFlags(IoModel* iomodel,bool** pflags);
    3636
    3737/*Connectivity*/
Note: See TracChangeset for help on using the changeset viewer.