source: issm/oecreview/Archive/16554-17801/ISSM-17134-17135.diff@ 17802

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

@@ -159 +159 @@
                 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;
+                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;
+                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;
+                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;
         }
         else{
                 _error_("mesh dimension not supported yet");

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

@@ -51 +51 @@
                 elementedges[2*3+0] = 1; elementedges[2*3+1] = 2; elementedges_markers[3] = 1;
                 elementedges[2*4+0] = 2; elementedges[2*4+1] = 0; elementedges_markers[4] = 1;
                 elementedges[2*5+0] = 0; elementedges[2*5+1] = 1; elementedges_markers[5] = 1;
-                elementedges[2*6+0] = 4; elementedges[2*6+1] = 5; elementedges_markers[6] = 3;
-                elementedges[2*7+0] = 5; elementedges[2*7+1] = 3; elementedges_markers[7] = 3;
-                elementedges[2*8+0] = 3; elementedges[2*8+1] = 4; elementedges_markers[8] = 3;
+                elementedges[2*6+0] = 4; elementedges[2*6+1] = 5; elementedges_markers[6] = 1;
+                elementedges[2*7+0] = 5; elementedges[2*7+1] = 3; elementedges_markers[7] = 1;
+                elementedges[2*8+0] = 3; elementedges[2*8+1] = 4; elementedges_markers[8] = 1;
         }
         else{
                 _error_("mesh dimension not supported yet");

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

@@ -13 +13 @@
         /*Intermediaries*/
         int   i,j,counter,vnodes,lid=0;
         int   id0 = iomodel->nodecounter;
+        bool *my_faces = NULL;
         bool *my_edges = NULL;
         bool *my_nodes = NULL;
         Node *node     = NULL;
@@ -143 +144 @@
                         break;
                 case P2xP4Enum:
                         EdgesPartitioning(&my_edges,iomodel);
+                        FacesPartitioning(&my_faces,iomodel);
                         for(i=0;i<iomodel->numberofvertices;i++){
                                 if(iomodel->my_vertices[i]){
                                         nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,approximation));
@@ -152 +154 @@
                         for(i=0;i<iomodel->numberofedges;i++){
                                 if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
                                         if(my_edges[i]){
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
+                                                node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
+                                                node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
+                                                node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
                                         }
                                         counter=counter+3;
                                 }
-                                else if(iomodel->edges[i*3+2]==1){/*Basal edges*/
+                                else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
                                         if(my_edges[i]){
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
+                                                node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
+                                        }
+                                        counter=counter+1;
+                                }
+                                else{
+                                        _error_("not supported");
+                                }
+                        }
+                        id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges;
+                        for(i=0;i<iomodel->numberoffaces;i++){
+                                if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
+                                        if(my_faces[i]){
+                                                node = new Node(id0+3*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
+                                                node = new Node(id0+3*i+2,counter+2,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+4,counter+4,lid++,0,iomodel,analysis,approximation);
+                                                node = new Node(id0+3*i+3,counter+3,lid++,0,iomodel,analysis,approximation);
                                                 nodes->AddObject(node);
                                         }
-                                        counter=counter+4;
+                                        counter=counter+3;
                                 }
-                                else if(iomodel->edges[i*3+2]==3){/*Top edges*/
-                                        if(my_edges[i]){
-                                                node = new Node(id0+iomodel->numberofvertices+4*i+1,counter+1,lid++,0,iomodel,analysis,approximation);
-                                                nodes->AddObject(node);
-                                        }
-                                        counter=counter+1;
+                                else if(iomodel->faces[i*6+5]==1){/*Triangular base/top*/
+                                        /*Nothing*/
                                 }
                                 else{
                                         _error_("not supported");
@@ -290 +299 @@
                 case OneLayerP4zEnum:
                         _assert_(approximation==FSApproximationEnum);
                         EdgesPartitioning(&my_edges,iomodel);
+                        FacesPartitioning(&my_faces,iomodel);
                         /*P2xP4 velocity*/
                         for(i=0;i<iomodel->numberofvertices;i++){
                                 if(iomodel->my_vertices[i]){
@@ -359 +369 @@
         }
 
         /*Clean up*/
+        xDelete<bool>(my_faces);
         xDelete<bool>(my_edges);
         xDelete<bool>(my_nodes);
 }

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

@@ -38 +38 @@
         IssmDouble *values           = NULL;
         bool        spcpresent       = false;
 
-        /*P2 finite elements*/
-        int   v1,v2;
+        /*Higher-order finite elements*/
+        int   v1,v2,v3,v4;
         bool *my_edges = NULL;
+        bool *my_faces = NULL;
 
         switch(finite_element){
                 case P1Enum: 
@@ -71 +72 @@
                         break;
                 case P2xP4Enum:
                         EdgesPartitioning(&my_edges,iomodel);
+                        FacesPartitioning(&my_faces,iomodel);
                         break;
                 default:
                         _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
@@ -125 +127 @@
                                                         v1 = iomodel->edges[3*i+0]-1;
                                                         v2 = iomodel->edges[3*i+1]-1;
                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+4*i+1,
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*i+1,
                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+4*i+2,
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+3*i+2,
                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+4*i+3,
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+3*i+3,
                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
                                                                 count=count+3;
                                                         }
                                                 }
                                         }
-                                        if(iomodel->edges[i*3+2]==1){/*Basal edges*/
+                                        if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
                                                 if(my_edges[i]){
                                                         v1 = iomodel->edges[3*i+0]-1;
                                                         v2 = iomodel->edges[3*i+1]-1;
                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+4*i+1,
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*i+1,
                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+4*i+2,
-                                                                                                dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+4*i+3,
-                                                                                                dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+4,iomodel->nodecounter+iomodel->numberofvertices+4*i+4,
-                                                                                                dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                count=count+4;
+                                                                count=count+1;
                                                         }
                                                 }
                                         }
-                                        if(iomodel->edges[i*3+2]==3){/*Top edges*/
-                                                if(my_edges[i]){
-                                                        v1 = iomodel->edges[3*i+0]-1;
-                                                        v2 = iomodel->edges[3*i+1]-1;
+                                }
+                                for(i=0;i<iomodel->numberoffaces;i++){
+                                        if(iomodel->faces[i*6+5]==2){/*Vertical quads*/
+                                                if(my_faces[i]){
+                                                        v1 = iomodel->faces[6*i+0]-1;
+                                                        v2 = iomodel->faces[6*i+1]-1;
                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
-                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+4*i+1,
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1,
                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
-                                                                count=count+1;
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+2,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2,
+                                                                                                dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+                                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+3,iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3,
+                                                                                                dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
+                                                                count=count+3;
                                                         }
                                                 }
                                         }

../trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -2833 +2833 @@
                 case P2xP4Enum:
                         numnodes         = 30;
                         penta_node_ids   = xNew<int>(numnodes);
-                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
-                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
-                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
-                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
-                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
-                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
-                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+0]+1;
-                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+1]+1;
-                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+2]+1;
-                        penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+1;
-                        penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+1;
-                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+1;
-                        penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+6]+1;
-                        penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+7]+1;
-                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+8]+1;
-                        penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+0]+2;
-                        penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+1]+2;
-                        penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+2]+2;
-                        penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+0]+3;
-                        penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+1]+3;
-                        penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+2]+3;
-                        penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+2;
-                        penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+2;
-                        penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+2;
-                        penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+3;
-                        penta_node_ids[25]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+3;
-                        penta_node_ids[26]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+3;
-                        penta_node_ids[27]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+3]+4;
-                        penta_node_ids[28]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+4]+4;
-                        penta_node_ids[29]=iomodel->nodecounter+iomodel->numberofvertices+4*iomodel->elementtoedgeconnectivity[9*index+5]+4;
+                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; /*Vertex 1*/
+                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1]; /*Vertex 2*/
+                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2]; /*Vertex 3*/
+                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3]; /*Vertex 4*/
+                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4]; /*Vertex 5*/
+                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5]; /*Vertex 6*/
+                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
+                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
+                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
+                        penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
+                        penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
+                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
+                        penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
+                        penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
+                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
+                        penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/
+                        penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/
+                        penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/
+                        penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 2/4 vertical edge 1*/
+                        penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 2/4 vertical edge 2*/
+                        penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 2/4 vertical edge 3*/
+                        penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+1; /* 1/4 vertical face 1*/
+                        penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+1; /* 1/4 vertical face 2*/
+                        penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+1; /* 1/4 vertical face 3*/
+                        penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+2; /* 2/4 vertical face 1*/
+                        penta_node_ids[25]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+2; /* 2/4 vertical face 2*/
+                        penta_node_ids[26]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+2; /* 2/4 vertical face 3*/
+                        penta_node_ids[27]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+3; /* 3/4 vertical face 1*/
+                        penta_node_ids[28]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+3; /* 3/4 vertical face 2*/
+                        penta_node_ids[29]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+3; /* 3/4 vertical face 3*/
                         break;
                 case P2Enum:
                         numnodes         = 15;