Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp (revision 17134) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateFaces.cpp (revision 17135) @@ -159,9 +159,9 @@ 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"); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp (revision 17134) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp (revision 17135) @@ -51,9 +51,9 @@ 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"); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17134) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17135) @@ -13,6 +13,7 @@ /*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,6 +144,7 @@ break; case P2xP4Enum: EdgesPartitioning(&my_edges,iomodel); + FacesPartitioning(&my_faces,iomodel); for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,approximation)); @@ -152,34 +154,41 @@ for(i=0;inumberofedges;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;inumberoffaces;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,6 +299,7 @@ case OneLayerP4zEnum: _assert_(approximation==FSApproximationEnum); EdgesPartitioning(&my_edges,iomodel); + FacesPartitioning(&my_faces,iomodel); /*P2xP4 velocity*/ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ @@ -359,6 +369,7 @@ } /*Clean up*/ + xDelete(my_faces); xDelete(my_edges); xDelete(my_nodes); } Index: ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 17134) +++ ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 17135) @@ -38,9 +38,10 @@ 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,6 +72,7 @@ break; case P2xP4Enum: EdgesPartitioning(&my_edges,iomodel); + FacesPartitioning(&my_faces,iomodel); break; default: _error_("Finite element "<edges[3*i+0]-1; v2 = iomodel->edges[3*i+1]-1; if(!xIsNan(spcdata[v1]) && !xIsNan(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(spcdata[v1]) && !xIsNan(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;inumberoffaces;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(spcdata[v1]) && !xIsNan(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; } } } Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17134) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17135) @@ -2833,36 +2833,36 @@ case P2xP4Enum: numnodes = 30; penta_node_ids = xNew(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;