Changeset 17122

Timestamp:

01/16/14 11:39:30 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added face creation and partitioning

Location:

issm/trunk-jpl/src/c

Files:

• classes/IoModel.h (modified) (1 diff)
• modules/ModelProcessorx/CreateFaces.cpp (modified) (6 diffs)
• modules/ModelProcessorx/FacesPartitioning.cpp (modified) (1 diff)
• modules/ModelProcessorx/ModelProcessorx.h (modified) (1 diff)

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

@@ -40 +40 @@
                 int *edges;
                 int *elementtoedgeconnectivity;
+                int *elementtofaceconnectivity;
                 int *singlenodetoelementconnectivity;
                 int *numbernodetoelementconnectivity;

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

@@ -5 +5 @@
 #include "../../classes/classes.h"
 #include "../../shared/shared.h"
-
-void CreateFaces(IoModel* iomodel){
+#include "./ModelProcessorx.h"
+
+void CreateFaces(IoModel* iomodel){/*{{{*/
 
         /*If faces are already present, exit*/
         if(iomodel->faces) return;
 
-        /*Check Iomodel properties*/
-        if(iomodel->meshtype!=Mesh2DhorizontalEnum)             _error_("only 2d model are supported");
+        /*Some checks*/
         if(iomodel->numberofvertices<3) _error_("not enough elements in mesh");
         _assert_(iomodel->elements);
+
+        /*Check Iomodel properties*/
+        if(iomodel->meshtype==Mesh2DhorizontalEnum){
+                /*Keep going*/
+        }
+        else if(iomodel->meshtype==Mesh3DEnum){
+                CreateFaces3d(iomodel);
+                return;
+        }
+        else{
+                _error_("mesh dimension not supported yet");
+        }
 
         /*Intermediaries*/
@@ -31 +43 @@
         /*Initialize chain*/
         int* head_minv = xNew<int>(iomodel->numberofvertices);
-        int* next_edge = xNew<int>(maxnbf);
+        int* next_face = xNew<int>(maxnbf);
         for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
 
@@ -57 +69 @@
                         /*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/
                         _assert_(v1>=0 && v1<iomodel->numberofvertices);
-                        for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
+                        for(int e=head_minv[v1]; e!=-1; e=next_face[e]){
                                 if(facestemp[e*4+1]==v2){
                                         exist = true;
@@ -76 +88 @@
 
                                 /*Update chain*/
-                                next_edge[nbf] = head_minv[v1];
+                                next_face[nbf] = head_minv[v1];
                                 head_minv[v1]  = nbf;
 
@@ -87 +99 @@
         /*Clean up*/
         xDelete<int>(head_minv);
-        xDelete<int>(next_edge);
+        xDelete<int>(next_face);
 
         /*Create final faces*/
@@ -109 +121 @@
         iomodel->faces         = faces;
         iomodel->numberoffaces = nbf;
-}
+}/*}}}*/
+void CreateFaces3d(IoModel* iomodel){/*{{{*/
+
+        /*Intermediaries*/
+        bool exist;
+        int  i,j,k,v0,cols;
+        int  maxnbf,nbf,elementnbf,elementnbv,facemaxnbv;
+        int *elementfaces         = NULL;
+        int *elementfaces_markers = NULL;
+
+        /*Maximum number of faces*/
+        maxnbf = 6*iomodel->numberofelements;
+
+        /*Initialize intermediaries*/
+        int*  facestemp = xNew<int>(maxnbf*6);         /*format: [vertex1 vertex2 vertex3 element1 element2 marker]    */
+        for(i=0;i<maxnbf;i++) facestemp[i*6+4]=-1;     /*Initialize last column of faces as -1 (boundary face)         */
+
+        /*Initialize chain*/
+        int* head_minv = xNew<int>(iomodel->numberofvertices);
+        int* next_face = xNew<int>(maxnbf);
+        for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
+
+        /*Mesh specific face indexing per element*/
+        if(iomodel->meshtype==Mesh3DEnum){
+                elementnbv = 6; /*Number of vertices per element*/
+                elementnbf = 5; /*Number of faces per element*/
+                facemaxnbv = 4; /*Maximum number of vertices per face*/
+                cols       = facemaxnbv + 1;
+                elementfaces         = xNew<int>(elementnbf*cols);
+                elementfaces_markers = xNew<int>(elementnbf);
+                /*2 triangles*/
+                elementfaces_markers[0] = 1;
+                elementfaces_markers[1] = 1;
+                elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0;  elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
+                elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3;  elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
+                /*3 quads*/
+                elementfaces_markers[2] = 2;
+                elementfaces_markers[3] = 2;
+                elementfaces_markers[4] = 2;
+                elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 0;  elementfaces[cols*2+2] = 1; elementfaces[cols*2+3] = 4;  elementfaces[cols*2+4] = 3;
+                elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 1;  elementfaces[cols*3+2] = 2; elementfaces[cols*3+3] = 5;  elementfaces[cols*3+4] = 4;
+                elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 2;  elementfaces[cols*4+2] = 0; elementfaces[cols*4+3] = 3;  elementfaces[cols*4+4] = 5;
+        }
+        else{
+                _error_("mesh dimension not supported yet");
+        }
+        int *element_face_connectivity = xNewZeroInit<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
+
+        /*Initialize number of faces and list of vertex indices*/
+        nbf = 0;
+        int* v = xNew<int>(facemaxnbv);
+
+        for(i=0;i<iomodel->numberofelements;i++){
+                for(j=0;j<elementnbf;j++){
+
+                        /*Get indices of current face*/
+                        for(k=0;k<elementfaces[cols*j+0];k++){
+                                v[k] = iomodel->elements[i*elementnbv + elementfaces[cols*j+k+1]] - 1;
+                        }
+
+                        /*Sort list of vertices*/
+                        HeapSort(v,elementfaces[cols*j+0]);
+                        v0 = v[0];
+
+                        /*This face a priori has not been processed yet*/
+                        exist = false;
+
+                        /*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/
+                        _assert_(v0>=0 && v0<iomodel->numberofvertices);
+                        for(int e=head_minv[v0]; e!=-1; e=next_face[e]){
+                                if(facestemp[e*6+1]==v[1]+1 && facestemp[e*6+2]==v[2]+1){
+                                        exist = true;
+                                        facestemp[e*6+4]=i+1;
+                                        element_face_connectivity[i*elementnbf+j]=e;
+                                        break;
+                                }
+                        }
+
+                        /*If this face is new, add it to the lists*/
+                        if(!exist){
+                                _assert_(nbf<maxnbf);
+
+                                /*Update faces*/
+                                facestemp[nbf*6+0] = v[0]+1;
+                                facestemp[nbf*6+1] = v[1]+1;
+                                facestemp[nbf*6+2] = v[2]+1;
+                                facestemp[nbf*6+3] = i+1;
+                                facestemp[nbf*6+5] = elementfaces_markers[j];
+
+                                /*Update Connectivity*/
+                                element_face_connectivity[i*elementnbf+j]=nbf;
+
+                                /*Update chain*/
+                                next_face[nbf] = head_minv[v0];
+                                head_minv[v0]  = nbf;
+
+                                /*Increase number of faces*/
+                                nbf++;
+                        }
+                }
+        }
+
+        /*Clean up*/
+        xDelete<int>(head_minv);
+        xDelete<int>(next_face);
+
+        /*Create final faces (now that we have the correct size)*/
+        int* faces = xNew<int>(nbf*6); /*vertex1 vertex2 vertex3 element1 element2 marker*/
+        xMemCpy<int>(faces,facestemp,nbf*6);
+        xDelete<int>(facestemp);
+
+        /*Assign output pointers*/
+        iomodel->faces                     = faces;
+        iomodel->numberoffaces             = nbf;
+        iomodel->elementtofaceconnectivity = element_face_connectivity;
+}/*}}}*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/FacesPartitioning.cpp

@@ -11 +11 @@
 
         /*Intermediaries*/
-        int  el1,el2;
-        bool my_face;
-
-        /*Check Iomodel properties*/
-        if(iomodel->meshtype!=Mesh2DhorizontalEnum) _error_("only 2d model are supported");
+        int elementnbf;
 
         /*Get faces and elements*/
         CreateFaces(iomodel);
+        _assert_(iomodel->elementtofaceconnectivity);
 
+        /*Mesh dependent variables*/
+        if(iomodel->meshtype==Mesh2DhorizontalEnum){
+                elementnbf = 3;
+        }
+        else if(iomodel->meshtype==Mesh2DverticalEnum){
+                elementnbf = 3;
+        }
+        else if(iomodel->meshtype==Mesh3DEnum){
+                elementnbf = 5;
+        }
+        else{
+                _error_("mesh dimension not supported yet");
+        }
         /*output: */
-        bool* my_faces=xNew<bool>(iomodel->numberoffaces);
+        bool* my_faces=xNewZeroInit<bool>(iomodel->numberoffaces);
 
-        for(int i=0;i<iomodel->numberoffaces;i++){
-
-                /*Get left and right elements*/
-                el1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
-                el2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
-
-                /*Check whether we should include this face (el2 is -2 for boundary faces)*/
-                my_face = iomodel->my_elements[el1];
-                if(!my_face && el2>=0){
-                        my_face = iomodel->my_elements[el2];
+        for(int i=0;i<iomodel->numberofelements;i++){
+                if(iomodel->my_elements[i]){
+                        for(int j=0;j<elementnbf;j++){
+                                _assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] >= 0);
+                                _assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] <  iomodel->numberoffaces);
+                                my_faces[iomodel->elementtofaceconnectivity[i*elementnbf+j]] = true;
+                        }
                 }
-
-                my_faces[i] = my_face;
         }
 

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

@@ -30 +30 @@
 
 /*Mesh properties*/
+void CreateEdges(IoModel* iomodel);
 void CreateFaces(IoModel* iomodel);
-void CreateEdges(IoModel* iomodel);
+void CreateFaces3d(IoModel* iomodel);
 
 /*Connectivity*/