Changeset 17844

Timestamp:

04/25/14 11:00:30 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: use edgeonboundary to figure out whether an edge should be constrained or not

Location:

issm/trunk-jpl/src/c/modules

Files:

• IoModelToConstraintsx/IoModelToConstraintsx.cpp (modified) (4 diffs)
• ModelProcessorx/CreateEdges.cpp (modified) (1 diff)
• ModelProcessorx/ModelProcessorx.h (modified) (1 diff)

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

@@ -44 +44 @@
         bool *my_edges = NULL;
         bool *my_faces = NULL;
+        bool *boundaryedge = NULL;
 
         switch(finite_element){
@@ -74 +75 @@
                                 FacesPartitioning(&my_faces,iomodel);
                         }
+                        EdgeOnBoundaryFlags(&boundaryedge,iomodel);
                         break;
                 case P2xP4Enum:
@@ -106 +108 @@
                                 }
                                 for(i=0;i<iomodel->numberofedges;i++){
-                                        if(my_edges[i]){
+                                        if(my_edges[i] && boundaryedge[i]){
                                                 v1 = iomodel->edges[3*i+0]-1;
                                                 v2 = iomodel->edges[3*i+1]-1;
@@ -466 +468 @@
         xDelete<bool>(my_edges);
         xDelete<bool>(my_faces);
+        xDelete<bool>(boundaryedge);
 }

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

@@ -144 +144 @@
         iomodel->numberofedges             = nbe;
 }/*}}}*/
-void EdgeOnBoundaryFlags(IoModel* iomodel,bool** pflags){
-
-        /*Get faces and edges*/
+void EdgeOnBoundaryFlags(bool** pflags,IoModel* iomodel){/*{{{*/
+
+        /*Intermediaries*/
+        bool isv1,isv2;
+        int  facenbv,v1,v2;
+        int  id_edge,id_element,id_face;
+        int  elementnbe;
+
+        /*Mesh dependent variables*/
+        switch(iomodel->meshelementtype){
+                case TriaEnum:  elementnbe = 3; break;
+                case TetraEnum: elementnbe = 6; break;
+                case PentaEnum: elementnbe = 9; break;
+                default:        _error_("mesh dimension not supported yet");
+        }
+
+        /*Get edges and allocate output*/
         if(!iomodel->edges) CreateEdges(iomodel);
-        if(!iomodel->faces) CreateFaces(iomodel);
-
-        /*Allocate output*/
-        bool* flags = xNew<bool>(iomodel->numberofedges);
-
-        /*Now, loop over the faces, whenever it is not connected to a second element, all edges are on boundary*/
-        for(int i=0;i<iomodel->numberoffaces;i++){
-
-                if(iomodel->faces[i*iomodel->facescols+4]==-1){
-
-                        /*Get all edges for this face*/
-                }
+        bool* flags = xNewZeroInit<bool>(iomodel->numberofedges);
+
+        if(iomodel->domaindim==2){
+
+                /*Count how many times an edge is found in elementtoedgeconnectivity*/
+                int* counter = xNewZeroInit<int>(iomodel->numberofedges);
+                for(int i=0;i<iomodel->numberofelements;i++){
+                        for(int j=0;j<elementnbe;j++){
+                                counter[iomodel->elementtoedgeconnectivity[elementnbe*i+j]] += 1;
+                        }
+                }
+
+                /*Now, loop over the egdes, whenever it is not connected to a second element, the edge is on boundary*/
+                for(int i=0;i<iomodel->numberofedges;i++){
+                        if(counter[i]==1) flags[i]=true;
+                }
+
+                /*Clean up*/
+                xDelete<int>(counter);
+        }
+        else if(iomodel->domaindim==3){
+
+                /*Get faces*/
+                if(!iomodel->faces) CreateFaces(iomodel);
+
+                /*Now, loop over the faces, whenever it is not connected to a second element, all edges are on boundary*/
+                for(int id_face=0;id_face<iomodel->numberoffaces;id_face++){
+
+                        if(iomodel->faces[id_face*iomodel->facescols+1]==-1){
+
+                                /*The face is connected to the element e only*/
+                                id_element = iomodel->faces[id_face*iomodel->facescols+0]-1;
+                                facenbv    = iomodel->faces[id_face*iomodel->facescols+3];
+
+                                /*Get all edges for this element*/
+                                for(int edge = 0; edge<elementnbe; edge++){
+
+                                        id_edge     = iomodel->elementtoedgeconnectivity[elementnbe*id_element+edge];
+                                        v1          = iomodel->edges[id_edge*3+0];
+                                        v2          = iomodel->edges[id_edge*3+1];
+
+                                        /*Test if v1 is in the face*/
+                                        isv1=false;
+                                        for(int i=0;i<facenbv;i++){
+                                                if(iomodel->faces[id_face*iomodel->facescols+4+i] == v1){
+                                                        isv1 = true; break;
+                                                }
+                                        }
+                                        if(!isv1) continue;
+
+                                        /*test if v2 is in the face*/
+                                        isv2=false;
+                                        for(int i=0;i<facenbv;i++){
+                                                if(iomodel->faces[id_face*iomodel->facescols+4+i] == v2){
+                                                        isv2 = true; break;
+                                                }
+                                        }
+
+                                        /*If v1 and v2 are found, this edge is on boundary*/
+                                        if(isv2) flags[id_edge] = true;
+                                }
+                        }
+                }
+        }
+        else{
+                _error_("dimension not supported");
         }
 
         /*Clean up and return*/
         *pflags = flags;
-}
+}/*}}}*/

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

@@ -33 +33 @@
 void CreateFaces(IoModel* iomodel);
 void CreateFaces3d(IoModel* iomodel);
-void EdgeOnBoundaryFlags(IoModel* iomodel,bool** pflags);
+void EdgeOnBoundaryFlags(bool** pflags,IoModel* iomodel);
 
 /*Connectivity*/