Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23576) @@ -545,45 +545,36 @@ else{ /*Coupling: we are going to create P1 Elements only*/ - - Node* node = NULL; - int lid=0; - if(!nodes) nodes = new Nodes(); - iomodel->FetchData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS", "md.flowequation.vertex_equation","md.stressbalance.referential"); if(isFS){ - /*P1+ velocity*/ + int* approximations = xNew(2*iomodel->numberofvertices+iomodel->numberofelements); for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - approximation=IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i])); - if(approximation==FSApproximationEnum) approximation=FSvelocityEnum; - nodes->AddObject(new Node(i+1,i,lid++,0,i,false,iomodel,StressbalanceAnalysisEnum,approximation)); - } + approximation = IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i])); + if(approximation==FSApproximationEnum) approximation=FSvelocityEnum; + approximations[i] = approximation; } - for(int i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - node = new Node(iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum); - node->Deactivate(); - nodes->AddObject(node); - } - } - /*P1 pressure*/ - for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - approximation=IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i])); - node = new Node(iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,0,i,false,iomodel,StressbalanceAnalysisEnum,FSpressureEnum); + for(int i=0;inumberofelements;i++) approximations[iomodel->numberofvertices+i] = FSvelocityEnum; + for(int i=0;inumberofvertices;i++) approximations[iomodel->numberofvertices+iomodel->numberofelements+i] = FSpressureEnum; + ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,MINIcondensedEnum,0,approximations); + xDelete(approximations); + + for(int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + int sid = node->Sid(); + if(sid>=iomodel->numberofvertices+iomodel->numberofelements){ + /*Constrain pressure if not FS*/ + int id = sid - (iomodel->numberofvertices+iomodel->numberofelements); + approximation=IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[id])); if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){ node->Deactivate(); } - nodes->AddObject(node); } } } else{ - for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(i+1,i,lid++,0,i,false,iomodel,StressbalanceAnalysisEnum,IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i])))); - } - } + int* approximations = xNew(iomodel->numberofvertices); + for(int i=0;inumberofvertices;i++) approximations[i] = IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i])); + ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,P1Enum,0,approximations); + xDelete(approximations); } iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS", Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 23576) @@ -80,5 +80,4 @@ /*Intermediaries*/ bool isSIA; - Node* node = NULL; /*Fetch parameters: */ @@ -95,20 +94,11 @@ } - for(int i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - - /*Create new node if is in this processor's partition*/ - node = new Node(i+1,i,lid++,0,i,false,iomodel,StressbalanceSIAAnalysisEnum,IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i]))); - - /*Deactivate node if not SIA*/ - if(IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[i]))!=SIAApproximationEnum){ - node->Deactivate(); - } - - /*Add to Nodes dataset*/ - nodes->AddObject(node); - } - } - + ::CreateNodes(nodes,iomodel,StressbalanceSIAAnalysisEnum,P1Enum,SIAApproximationEnum); + for(int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + int sid = node->Sid(); + int approximation=IoCodeToEnumVertexEquation(reCast(iomodel->Data("md.flowequation.vertex_equation")[sid])); + if(approximation!=SIAApproximationEnum) node->Deactivate(); + } iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS","md.flowequation.vertex_equation","md.stressbalance.referential"); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23576) @@ -2918,7 +2918,7 @@ penta_node_ids[ 4]=iomodel->elements[6*index+4]; penta_node_ids[ 5]=iomodel->elements[6*index+5]; - penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; - penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; - penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; + penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; + penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; + penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; break; case P1xP3Enum: @@ -2931,10 +2931,10 @@ penta_node_ids[ 4]=iomodel->elements[6*index+4]; penta_node_ids[ 5]=iomodel->elements[6*index+5]; - penta_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+0]+1; - penta_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+1]+1; - penta_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+2]+1; - penta_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+0]+2; - penta_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+1]+2; - penta_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+2]+2; + penta_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; + penta_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; + penta_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; + penta_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; + penta_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; + penta_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; break; case P2xP1Enum: @@ -2947,10 +2947,10 @@ penta_node_ids[ 4]=iomodel->elements[6*index+4]; penta_node_ids[ 5]=iomodel->elements[6*index+5]; - penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; - penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; - penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; - penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; - penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; - penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; + penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+0]+1; + penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+1]+1; + penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+2]+1; + penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+3]+1; + penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+4]+1; + penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+5]+1; break; case P1xP4Enum: @@ -2963,13 +2963,13 @@ penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/ penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/ - penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/ - penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/ - penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/ - penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/ - penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/ - penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/ - penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 3/4 vertical edge 1*/ - penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 3/4 vertical edge 2*/ - penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 3/4 vertical edge 3*/ + penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /*mid vertical edge 1*/ + penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /*mid vertical edge 2*/ + penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /*mid vertical edge 3*/ + penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 1/4 vertical edge 1*/ + penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 1/4 vertical edge 2*/ + penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 1/4 vertical edge 3*/ + penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+3; /* 3/4 vertical edge 1*/ + penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+3; /* 3/4 vertical edge 2*/ + penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+3; /* 3/4 vertical edge 3*/ break; case P2xP4Enum: @@ -2982,28 +2982,28 @@ penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/ penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/ - penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/ - penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/ - penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/ - penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/ - penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/ - penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/ - penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/ - penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/ - penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/ - penta_node_ids[15]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/ - penta_node_ids[16]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/ - penta_node_ids[17]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/ - penta_node_ids[18]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 3/4 vertical edge 1*/ - penta_node_ids[19]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 3/4 vertical edge 2*/ - penta_node_ids[20]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 3/4 vertical edge 3*/ - penta_node_ids[21]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+1; /* 1/4 vertical face 1*/ - penta_node_ids[22]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+1; /* 1/4 vertical face 2*/ - penta_node_ids[23]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+1; /* 1/4 vertical face 3*/ - penta_node_ids[24]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+2; /* 2/4 vertical face 1*/ - penta_node_ids[25]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+2; /* 2/4 vertical face 2*/ - penta_node_ids[26]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+2; /* 2/4 vertical face 3*/ - penta_node_ids[27]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+3; /* 3/4 vertical face 1*/ - penta_node_ids[28]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+3; /* 3/4 vertical face 2*/ - penta_node_ids[29]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+3; /* 3/4 vertical face 3*/ + penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/ + penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/ + penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/ + penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/ + penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/ + penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/ + penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/ + penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/ + penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/ + penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /* 1/4 vertical edge 1*/ + penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /* 1/4 vertical edge 2*/ + penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /* 1/4 vertical edge 3*/ + penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 3/4 vertical edge 1*/ + penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 3/4 vertical edge 2*/ + penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 3/4 vertical edge 3*/ + penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; /* 1/4 vertical face 1*/ + penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; /* 1/4 vertical face 2*/ + penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; /* 1/4 vertical face 3*/ + penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+2; /* 2/4 vertical face 1*/ + penta_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+2; /* 2/4 vertical face 2*/ + penta_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+2; /* 2/4 vertical face 3*/ + penta_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+3; /* 3/4 vertical face 1*/ + penta_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+3; /* 3/4 vertical face 2*/ + penta_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+3; /* 3/4 vertical face 3*/ break; case P2Enum: @@ -3025,7 +3025,7 @@ penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; - penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1; - penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1; - penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1; + penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; + penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; + penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; break; case P2bubbleEnum: case P2bubblecondensedEnum: @@ -3105,14 +3105,14 @@ penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; - penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1; - penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1; - penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1; - - penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+0]; - penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+1]; - penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+2]; - penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+3]; - penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+4]; - penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+5]; + penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; + penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; + penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; + + penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+0]; + penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+1]; + penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+2]; + penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+3]; + penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+4]; + penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+5]; break; case LATaylorHoodEnum: @@ -3147,35 +3147,35 @@ penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/ penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/ - penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/ - penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/ - penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/ - penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/ - penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/ - penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/ - penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/ - penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/ - penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/ - penta_node_ids[15]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/ - penta_node_ids[16]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/ - penta_node_ids[17]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/ - penta_node_ids[18]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 2/4 vertical edge 1*/ - penta_node_ids[19]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 2/4 vertical edge 2*/ - penta_node_ids[20]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 2/4 vertical edge 3*/ - penta_node_ids[21]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+1; /* 1/4 vertical face 1*/ - penta_node_ids[22]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+1; /* 1/4 vertical face 2*/ - penta_node_ids[23]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+1; /* 1/4 vertical face 3*/ - penta_node_ids[24]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+2; /* 2/4 vertical face 1*/ - penta_node_ids[25]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+2; /* 2/4 vertical face 2*/ - penta_node_ids[26]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+2; /* 2/4 vertical face 3*/ - penta_node_ids[27]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+3; /* 3/4 vertical face 1*/ - penta_node_ids[28]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+3; /* 3/4 vertical face 2*/ - penta_node_ids[29]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+3; /* 3/4 vertical face 3*/ - - penta_node_ids[30]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+0]; - penta_node_ids[31]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+1]; - penta_node_ids[32]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+2]; - penta_node_ids[33]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+3]; - penta_node_ids[34]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+4]; - penta_node_ids[35]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+5]; + penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/ + penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/ + penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/ + penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/ + penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/ + penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/ + penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/ + penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/ + penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/ + penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /* 1/4 vertical edge 1*/ + penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /* 1/4 vertical edge 2*/ + penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /* 1/4 vertical edge 3*/ + penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 3/4 vertical edge 1*/ + penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 3/4 vertical edge 2*/ + penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 3/4 vertical edge 3*/ + penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; /* 1/4 vertical face 1*/ + penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; /* 1/4 vertical face 2*/ + penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; /* 1/4 vertical face 3*/ + penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+2; /* 2/4 vertical face 1*/ + penta_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+2; /* 2/4 vertical face 2*/ + penta_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+2; /* 2/4 vertical face 3*/ + penta_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+3; /* 3/4 vertical face 1*/ + penta_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+3; /* 3/4 vertical face 2*/ + penta_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+3; /* 3/4 vertical face 3*/ + + penta_node_ids[30]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+0]; + penta_node_ids[31]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+1]; + penta_node_ids[32]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+2]; + penta_node_ids[33]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+3]; + penta_node_ids[34]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+4]; + penta_node_ids[35]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+5]; break; case CrouzeixRaviartEnum: Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23576) @@ -2707,4 +2707,5 @@ } SpcNodesx(new_nodes_list[i],new_constraints_list[i],this->parameters,analysis_type); + new_nodes_list[i]->FlagClones(analysis_type);/*FIXME: should be removed !*/ NodesDofx(new_nodes_list[i],this->parameters,analysis_type); } Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 23576) @@ -113,5 +113,5 @@ } - int count=0; + int count =0; if(M==iomodel->numberofvertices){ switch(finite_element){ @@ -127,4 +127,37 @@ break; case P2Enum: + for(i=0;inumberofvertices;i++){ + if((iomodel->my_vertices[i])){ + if (!xIsNan(spcdata[i])){ + constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type)); + count++; + } + } + } + for(i=0;inumberofedges;i++){ + if(iomodel->my_edges[i] && boundaryedge[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(count+1,iomodel->numberofvertices+i+1, + dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + count++; + } + } + } + if(iomodel->meshelementtype==PentaEnum){ + for(i=0;inumberofverticalfaces;i++){ + if(iomodel->my_vfaces[i]){ + value=0.; + for(j=0;j<4;j++) value += spcdata[iomodel->verticalfaces[i*4+j] -1]/4.; + if(!xIsNan(value)){ + constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+iomodel->numberofedges+i+1,dof,value,analysis_type)); + count++; + } + } + } + } + break; + case P2bubbleEnum: for(i=0;inumberofvertices;i++){ if((iomodel->my_vertices[i])){ @@ -166,45 +199,4 @@ } } - break; - case P2bubbleEnum: - for(i=0;inumberofvertices;i++){ - if((iomodel->my_vertices[i])){ - if (!xIsNan(spcdata[i])){ - constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type)); - count++; - } - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i] && boundaryedge[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(count+1,iomodel->numberofvertices+i+1, - dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); - count++; - } - } - } - if(iomodel->meshelementtype==PentaEnum){ - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - numfacevertices = iomodel->faces[i*iomodel->facescols+3]; - value=0.; - for(j=0;jfaces[i*iomodel->facescols+4+j] -1]; - } - value = value/reCast(numfacevertices); - if(!xIsNan(value)){ - constraints->AddObject(new SpcStatic(count+1, - iomodel->numberofvertices+iomodel->numberofedges+i+1, - dof,value,analysis_type)); - count++; - } - } - } - } - } for(i=0;inumberofelements;i++){ if(iomodel->my_elements[i]){ @@ -233,52 +225,43 @@ } for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){/*Vertical edges*/ - if(iomodel->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(count+1,iomodel->numberofvertices+3*i+1, - dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); - constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2, - dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); - constraints->AddObject(new SpcStatic(count+3,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){/*Horizontal edges*/ - if(iomodel->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(count+1,iomodel->numberofvertices+3*i+1, - dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); - count=count+1; - } - } - } - } - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - numfacevertices = iomodel->faces[i*iomodel->facescols+3]; - value=0.; - for(j=0;jfaces[i*iomodel->facescols+4+j] -1]; - } - value = value/reCast(numfacevertices); - if(!xIsNan(value)){ - constraints->AddObject(new SpcStatic(count+1, - iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1, - dof,value,analysis_type)); - constraints->AddObject(new SpcStatic(count+2, - iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2, - dof,value,analysis_type)); - constraints->AddObject(new SpcStatic(count+3, - iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3, - dof,value,analysis_type)); - count=count+3; - } + if(iomodel->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(count+1,iomodel->numberofvertices+i+1, + dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); + count++; + } + } + } + for(i=0;inumberofverticaledges;i++){ + if(iomodel->my_vedges[i]){ + v1 = iomodel->verticaledges[2*i+0]-1; + v2 = iomodel->verticaledges[2*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + /*FIXME: should not be 1/2 but 1/4 and 3/4!*/ + constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+iomodel->numberofedges+2*i+1, + dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); + constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+iomodel->numberofedges+2*i+2, + dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); + count=count+2; + } + } + } + for(i=0;inumberofverticalfaces;i++){ + if(iomodel->my_vfaces[i]){ + value=0.; + for(j=0;j<4;j++) value += spcdata[iomodel->verticalfaces[i*4+j]-1]/4.; + if(!xIsNan(value)){ + constraints->AddObject(new SpcStatic(count+1, + iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+1, + dof,value,analysis_type)); + constraints->AddObject(new SpcStatic(count+2, + iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+2, + dof,value,analysis_type)); + constraints->AddObject(new SpcStatic(count+3, + iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+3, + dof,value,analysis_type)); + count=count+3; } } @@ -326,14 +309,11 @@ } } - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){ - if(iomodel->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(count+1,iomodel->numberofvertices+i+1, - dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); - count++; - } + for(i=0;inumberofverticaledges;i++){ + if(iomodel->my_vedges[i]){ + v1 = iomodel->verticaledges[2*i+0]-1; + v2 = iomodel->verticaledges[2*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + count++; } } @@ -349,16 +329,14 @@ } } - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){ - if(iomodel->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(count+1,iomodel->numberofvertices+2*i+1, - dof,2./3.*spcdata[v1]+1./3.*spcdata[v2],analysis_type)); - constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+2*i+2, - dof,1./3.*spcdata[v1]+2./3.*spcdata[v2],analysis_type)); - count=count+2; - } + for(i=0;inumberofverticaledges;i++){ + if(iomodel->my_vedges[i]){ + v1 = iomodel->verticaledges[2*i+0]-1; + v2 = iomodel->verticaledges[2*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+2*i+1, + dof,2./3.*spcdata[v1]+1./3.*spcdata[v2],analysis_type)); + constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+2*i+2, + dof,1./3.*spcdata[v1]+2./3.*spcdata[v2],analysis_type)); + count=count+2; } } @@ -374,18 +352,16 @@ } } - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){/*Vertical edges*/ - if(iomodel->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(count+1,iomodel->numberofvertices+3*i+1, - dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); - constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2, - dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); - constraints->AddObject(new SpcStatic(count+3,iomodel->numberofvertices+3*i+3, - dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); - count=count+3; - } + for(i=0;inumberofverticaledges;i++){ + if(iomodel->my_vedges[i]){ + v1 = iomodel->verticaledges[2*i+0]-1; + v2 = iomodel->verticaledges[2*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+3*i+1, + dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); + constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2, + dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); + constraints->AddObject(new SpcStatic(count+3,iomodel->numberofvertices+3*i+3, + dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type)); + count=count+3; } } @@ -401,14 +377,11 @@ } } - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]!=2){ - if(iomodel->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(count+1,iomodel->numberofvertices+i+1, - dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); - count++; - } + for(i=0;inumberofhorizontaledges;i++){ + if(iomodel->my_hedges[i]){ + v1 = iomodel->horizontaledges[2*i+0]-1; + v2 = iomodel->horizontaledges[2*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + count++; } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 23576) @@ -9,595 +9,680 @@ #include "./ModelProcessorx.h" -void CreateNodes(Nodes* nodes, IoModel* iomodel,int analysis,int finite_element,int approximation){ +void CreateNodes(Nodes* nodes, IoModel* iomodel,int analysis,int finite_element,int approximation,int* approximations){ /*Intermediaries*/ - int i,j,counter,vnodes,lid=0; - int numberoffaces,elementnbv; - bool *my_nodes = NULL; - Node *node = NULL; - - int id0 = 0; - switch(finite_element){ - case P1Enum: - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - break; - - case P1DGEnum: - DiscontinuousGalerkinNodesPartitioning(&my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel); - for(i=0;inumberofelements;i++){ - for(j=0;j<3;j++){ - if(my_nodes[3*i+j]){ - nodes->AddObject(new Node(id0+3*i+j+1,id0+3*i+j,lid++,0,iomodel->elements[3*i+j]-1,false,iomodel,analysis,approximation)); - - } - } - } - break; - - case P1bubbleEnum: - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation)); - } - } - break; - - case P1bubblecondensedEnum: - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - node = new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation); - node->HardDeactivate(); - nodes->AddObject(node); - } - } - break; - - case P1xP2Enum: - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - - counter = iomodel->numberofvertices; - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter++; - } - } - break; - - case P1xP3Enum: - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - - counter = iomodel->numberofvertices; - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+2*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+2*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter=counter+2; - } - } - break; - case P1xP4Enum: - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - counter = iomodel->numberofvertices; - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){/*Vertical edges*/ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter=counter+3; - } - } - break; - - case P2xP1Enum: - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - - counter = iomodel->numberofvertices; - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]!=2){ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter++; - } - } - break; - - case P2Enum: - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation)); - } - } - id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; - if(iomodel->meshelementtype==PentaEnum){ + int numnodes; + const int MAXCONNECTIVITY = 5; + int element_numnodes; + int element_node_ids[40] = {0}; + + /*Get partitioning variables*/ + int my_rank = IssmComm::GetRank(); + int num_procs = IssmComm::GetSize(); + int* epart = iomodel->epart; + + + /*Determine how many nodes we have in total*/ + /*Define nodes sids for each element*/ + if(iomodel->meshelementtype==TriaEnum){ + switch(finite_element){ + case P1Enum: + numnodes = iomodel->numberofvertices; + break; + case P1DGEnum: + numnodes = 3*iomodel->numberofelements; + break; + case P1bubbleEnum: case P1bubblecondensedEnum: + numnodes = iomodel->numberofvertices+iomodel->numberofelements; + break; + case P2Enum: + EdgesPartitioning(iomodel); + numnodes = iomodel->numberofvertices+iomodel->numberofedges; + break; + case MINIEnum: case MINIcondensedEnum: + /*P1+ velocity (bubble statically condensed), P1 pressure*/ + numnodes = (iomodel->numberofvertices+iomodel->numberofelements) + (iomodel->numberofvertices); + break; + case TaylorHoodEnum: case XTaylorHoodEnum: + /*P2 velocity, P1 pressure*/ + EdgesPartitioning(iomodel); + numnodes = (iomodel->numberofvertices+iomodel->numberofedges) + (iomodel->numberofvertices); + break; + case LATaylorHoodEnum: + /*P2 velocity*/ + EdgesPartitioning(iomodel); + numnodes = (iomodel->numberofvertices+iomodel->numberofedges); + break; + case CrouzeixRaviartEnum: + /*P2b velocity, P1 DG pressure*/ + EdgesPartitioning(iomodel); + numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements) + (3*iomodel->numberofelements); + break; + case LACrouzeixRaviartEnum: + /*P2b velocity*/ + EdgesPartitioning(iomodel); + numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements); + break; + default: + _error_("Finite element "<meshelementtype==PentaEnum){ + switch(finite_element){ + case P1Enum: + numnodes = iomodel->numberofvertices; + break; + case P1bubbleEnum: case P1bubblecondensedEnum: + numnodes = iomodel->numberofvertices+iomodel->numberofelements; + break; + case P1xP2Enum: + EdgesPartitioning(iomodel); + numnodes = iomodel->numberofvertices+iomodel->numberofverticaledges; + break; + case P1xP3Enum: + EdgesPartitioning(iomodel); + numnodes = iomodel->numberofvertices+2*iomodel->numberofverticaledges; + break; + case P1xP4Enum: + EdgesPartitioning(iomodel); + numnodes = iomodel->numberofvertices+3*iomodel->numberofverticaledges; + break; + case P2xP1Enum: + EdgesPartitioning(iomodel); + numnodes = iomodel->numberofvertices+iomodel->numberofhorizontaledges; + break; + case P2Enum: + EdgesPartitioning(iomodel); FacesPartitioning(iomodel); - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); + numnodes = iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces; + break; + case P2xP4Enum: + EdgesPartitioning(iomodel); + FacesPartitioning(iomodel); + numnodes = iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces; + break; + case MINIEnum: case MINIcondensedEnum: + /*P1+ velocity (bubble statically condensed), P1 pressure*/ + numnodes = (iomodel->numberofvertices+iomodel->numberofelements) + (iomodel->numberofvertices); + break; + case TaylorHoodEnum: case XTaylorHoodEnum: + /*P2 velocity, P1 pressure*/ + EdgesPartitioning(iomodel); + FacesPartitioning(iomodel); + numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces) + (iomodel->numberofvertices); + break; + case OneLayerP4zEnum: + /*P2xP4 velocity, P1 pressure*/ + EdgesPartitioning(iomodel); + FacesPartitioning(iomodel); + numnodes = (iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces) + (iomodel->numberofvertices); + break; + default: + _error_("Finite element "<meshelementtype) <<" not supported yet"); + } + + /*create matrix that keeps track of all ranks that have node i, and initialize as -1 (Common to all CPUs)*/ + int* nodes_ranks = xNew(MAXCONNECTIVITY*numnodes); + for(int i=0;i(numnodes); + + /*Create vector of size total numnodes, initialized with -1, that will keep track of local ids*/ + int* nodes_lids = xNew(numnodes); + for(int i=0;i(numnodes); + if(approximations){ + for(int i=0;inumberofelements;i++){ + + /*Define nodes sids for each element*/ + if(iomodel->meshelementtype==TriaEnum){ + switch(finite_element){ + case P1Enum: + element_numnodes=3; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + break; + case P1bubbleEnum: case P1bubblecondensedEnum: + element_numnodes=4; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+i; + break; + case P1DGEnum: + element_numnodes=3; + element_node_ids[0]=3*i+0; + element_node_ids[1]=3*i+1; + element_node_ids[2]=3*i+2; + break; + case P2Enum: + element_numnodes = 6; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0]; + element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1]; + element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2]; + break; + case MINIEnum: case MINIcondensedEnum: + element_numnodes = 4+3; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+i; + for(int n=0;n<4;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[4]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+0]-1; + element_node_ids[5]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+1]-1; + element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+2]-1; + for(int n=4;n<7;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; + break; + case TaylorHoodEnum: case XTaylorHoodEnum: + element_numnodes = 6+3; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0]; + element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1]; + element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2]; + for(int n=0;n<6;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+0]-1; + element_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+1]-1; + element_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+2]-1; + for(int n=6;n<9;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; + break; + case LATaylorHoodEnum: + element_numnodes = 6; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0]; + element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1]; + element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2]; + for(int n=0;n<6;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + break; + case CrouzeixRaviartEnum: + element_numnodes = 7+3; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0]; + element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1]; + element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2]; + element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+i; + for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+0; + element_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+1; + element_node_ids[9]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+2; + for(int n=7;n<10;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; + break; + case LACrouzeixRaviartEnum: + element_numnodes = 7; + element_node_ids[0]=iomodel->elements[3*i+0]-1; + element_node_ids[1]=iomodel->elements[3*i+1]-1; + element_node_ids[2]=iomodel->elements[3*i+2]-1; + element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0]; + element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1]; + element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2]; + element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+i; + for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + break; + default: + _error_("Finite element "<meshelementtype==PentaEnum){ + switch(finite_element){ + case P1Enum: + element_numnodes=6; + element_node_ids[0]=iomodel->elements[6*i+0]-1; + element_node_ids[1]=iomodel->elements[6*i+1]-1; + element_node_ids[2]=iomodel->elements[6*i+2]-1; + element_node_ids[3]=iomodel->elements[6*i+3]-1; + element_node_ids[4]=iomodel->elements[6*i+4]-1; + element_node_ids[5]=iomodel->elements[6*i+5]-1; + break; + case P1bubbleEnum: case P1bubblecondensedEnum: + element_numnodes=7; + element_node_ids[0]=iomodel->elements[6*i+0]-1; + element_node_ids[1]=iomodel->elements[6*i+1]-1; + element_node_ids[2]=iomodel->elements[6*i+2]-1; + element_node_ids[3]=iomodel->elements[6*i+3]-1; + element_node_ids[4]=iomodel->elements[6*i+4]-1; + element_node_ids[5]=iomodel->elements[6*i+5]-1; + element_node_ids[6]=iomodel->numberofvertices+i; + break; + case P1xP2Enum: + element_numnodes=9; + element_node_ids[0]=iomodel->elements[6*i+0]-1; + element_node_ids[1]=iomodel->elements[6*i+1]-1; + element_node_ids[2]=iomodel->elements[6*i+2]-1; + element_node_ids[3]=iomodel->elements[6*i+3]-1; + element_node_ids[4]=iomodel->elements[6*i+4]-1; + element_node_ids[5]=iomodel->elements[6*i+5]-1; + element_node_ids[6]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+0]; + element_node_ids[7]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+1]; + element_node_ids[8]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+2]; + break; + case P1xP3Enum: + element_numnodes=12; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; + element_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+0; + element_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+0; + element_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+0; + element_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; + element_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; + element_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; + break; + case P1xP4Enum: + element_numnodes=15; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/ + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/ + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/ + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/ + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/ + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/ + element_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+0; /*mid vertical edge 1*/ + element_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+0; /*mid vertical edge 2*/ + element_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+0; /*mid vertical edge 3*/ + element_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 1/4 vertical edge 1*/ + element_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 1/4 vertical edge 2*/ + element_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 1/4 vertical edge 3*/ + element_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+2; /* 3/4 vertical edge 1*/ + element_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+2; /* 3/4 vertical edge 2*/ + element_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+2; /* 3/4 vertical edge 3*/ + break; + case P2xP1Enum: + element_numnodes=12; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; + element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+0]; + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+1]; + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+2]; + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+3]; + element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+4]; + element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+5]; + break; + case P2Enum: + element_numnodes = 18; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; + element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; + element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; + element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; + element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; + element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; + element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; + element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+0]; + element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+1]; + element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+2]; + break; + case P2xP4Enum: + element_numnodes = 30; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/ + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/ + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/ + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/ + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/ + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/ + element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; /*mid vertical edge 1*/ + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; /*mid vertical edge 2*/ + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; /*mid vertical edge 3*/ + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; /*mid basal edge 1*/ + element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; /*mid basal edge 2*/ + element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; /*mid basal edge 3*/ + element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; /*mid top edge 1*/ + element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; /*mid top edge 2*/ + element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; /*mid top edge 3*/ + element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]; /* 1/4 vertical edge 1*/ + element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]; /* 1/4 vertical edge 2*/ + element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]; /* 1/4 vertical edge 3*/ + element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 3/4 vertical edge 1*/ + element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 3/4 vertical edge 2*/ + element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 3/4 vertical edge 3*/ + element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+0; /* 1/4 vertical face 1*/ + element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+0; /* 1/4 vertical face 2*/ + element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+0; /* 1/4 vertical face 3*/ + element_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+1; /* 2/4 vertical face 1*/ + element_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+1; /* 2/4 vertical face 2*/ + element_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+1; /* 2/4 vertical face 3*/ + element_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+2; /* 3/4 vertical face 1*/ + element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/ + element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/ + break; + case MINIEnum: case MINIcondensedEnum: + element_numnodes = 7+6; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; + element_node_ids[ 6]=iomodel->numberofvertices+i; + if(!approximations) for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+0]-1; + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+1]-1; + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+2]-1; + element_node_ids[10]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+3]-1; + element_node_ids[11]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+4]-1; + element_node_ids[12]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+5]-1; + if(!approximations) for(int n=7;n<13;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; + break; + case TaylorHoodEnum: case XTaylorHoodEnum: + element_numnodes = 18+6; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; + element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; + element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; + element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; + element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; + element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; + element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; + element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+0]; + element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+1]; + element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+2]; + for(int n=0;n<18;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+0]-1; + element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+1]-1; + element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+2]-1; + element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+3]-1; + element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+4]-1; + element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+5]-1; + for(int n=18;n<24;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; + break; + case OneLayerP4zEnum: + element_numnodes = 30+6; + element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/ + element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/ + element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/ + element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/ + element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/ + element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/ + element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; /*mid vertical edge 1*/ + element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; /*mid vertical edge 2*/ + element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; /*mid vertical edge 3*/ + element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; /*mid basal edge 1*/ + element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; /*mid basal edge 2*/ + element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; /*mid basal edge 3*/ + element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; /*mid top edge 1*/ + element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; /*mid top edge 2*/ + element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; /*mid top edge 3*/ + element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]; /* 1/4 vertical edge 1*/ + element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]; /* 1/4 vertical edge 2*/ + element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]; /* 1/4 vertical edge 3*/ + element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 3/4 vertical edge 1*/ + element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 3/4 vertical edge 2*/ + element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 3/4 vertical edge 3*/ + element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+0; /* 1/4 vertical face 1*/ + element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+0; /* 1/4 vertical face 2*/ + element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+0; /* 1/4 vertical face 3*/ + element_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+1; /* 2/4 vertical face 1*/ + element_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+1; /* 2/4 vertical face 2*/ + element_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+1; /* 2/4 vertical face 3*/ + element_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+2; /* 3/4 vertical face 1*/ + element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/ + element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/ + for(int n=0;n<30;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum; + element_node_ids[30]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+0]-1; + element_node_ids[31]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+1]-1; + element_node_ids[32]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+2]-1; + element_node_ids[33]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+3]-1; + element_node_ids[34]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+4]-1; + element_node_ids[35]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+5]-1; + for(int n=30;n<36;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum; + break; + default: + _error_("Finite element "<meshelementtype) <<" not supported yet"); + } + + for(int j=0;jMAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)"); + } + } + } + if(finite_element==P1DGEnum){/* Special case for DG...{{{*/ + int node_list[4]; + if(!iomodel->domaintype==Domain2DhorizontalEnum) _error_("not implemented yet"); + CreateEdges(iomodel); + CreateFaces(iomodel); + for(int i=0;inumberoffaces;i++){ + int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2] + if(e2!=-2){ + if(epart[e1]!=epart[e2]){ + int i1=iomodel->faces[4*i+0]; + int i2=iomodel->faces[4*i+1]; + int pos=-1; + for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j; + if( pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;} + else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;} + else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;} + else _error_("not supposed to happen"); + pos=-1; + for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j; + if( pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;} + else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;} + else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;} + else _error_("not supposed to happen"); + for(int j=0;j<4;j++){ + int nid = node_list[j]; + bool found = false; + for(int k=0;kMAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)"); } } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - } - break; - case P2bubbleEnum: - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation)); - } - } - id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; - if(iomodel->meshelementtype==PentaEnum){ - FacesPartitioning(iomodel); - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); + for(int j=0;j<4;j++){ + int nid = node_list[j]; + bool found = false; + for(int k=0;kMAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)"); } } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ + } + } + } + }/*}}}*/ + //if(my_rank==0) printarray(nodes_ranks,numnodes,MAXCONNECTIVITY); + + /*Now, Count how many clones we have with other partitions*/ + int* common_send = xNew(num_procs); + int* common_recv = xNew(num_procs); + int** common_send_ids = xNew(num_procs); + int** common_recv_ids = xNew(num_procs); + + /*First step: allocate, Step 2: populate*/ + for(int step=0;step<2;step++){ + if(step==1){ + /*Allocate send and receive arrays of ids now*/ + for(int i=0;i=0 && common_recv[i]>=0); + common_send_ids[i] = xNew(common_send[i]); + common_recv_ids[i] = xNew(common_recv[i]); + } + } + /*Re/Initialize counters to 0*/ + for(int i=0;i=0); + int rank = nodes_ranks[MAXCONNECTIVITY*i+j]; + if(step==1){ + common_send_ids[rank][common_send[rank]] = nodes_lids[i]; + } + common_send[rank]++; } - else{ - _error_("not supported"); - } - } - id0 = id0+iomodel->numberoffaces; - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - nodes->AddObject(new Node(id0+i+1,id0+i,lid++,0,0,false,iomodel,analysis,approximation)); - } - } - break; - case P2xP4Enum: - EdgesPartitioning(iomodel); - FacesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation)); - } - } - counter = iomodel->numberofvertices; - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){/*Vertical edges*/ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter=counter+3; - } - else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter=counter+1; - } - else{ - _error_("not supported"); - } - } - id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges; - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - node = new Node(id0+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,approximation); - nodes->AddObject(node); - } - counter=counter+3; - } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - break; - - /*Stokes elements*/ - case P1P1Enum: - _assert_(approximation==FSApproximationEnum); - /*P1 velocity*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - /*P1 pressure*/ - vnodes = id0+iomodel->numberofvertices; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum)); - } - } - break; - case P1P1GLSEnum: - _assert_(approximation==FSApproximationEnum); - /*P1 velocity*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - /*P1 pressure*/ - vnodes = id0+iomodel->numberofvertices; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum)); - } - } - break; - case MINIcondensedEnum: - _assert_(approximation==FSApproximationEnum); - /*P1+ velocity (bubble statically condensed)*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - node = new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); + } + } + else{ + /*3. It is a slave, record that we need to receive for this cpu*/ + int rank = nodes_ranks[MAXCONNECTIVITY*i+0]; + if(step==1){ + common_recv_ids[rank][common_recv[rank]] = nodes_lids[i]; + } + common_recv[rank]++; + } + } + } + xDelete(nodes_proc_count); + + /*Final step: prepare pids (parallel ids), first count number of masters for each proc*/ + int* num_masters = xNewZeroInit(num_procs); + for(int i=0;i(num_procs); + for(int i=0;i(num_masters); + /*Now create pids*/ + int* nodes_pids = xNew(numnodes); + int* rank_counters = xNewZeroInit(numnodes); + for(int i=0;i(rank_counters); + xDelete(rank_offsets); + + /*Go ahead and create vertices now that we have all we need*/ + for(int i=0;inumberofvertices) io_index = i; + Node* node=new Node(i+1,i,nodes_lids[i],nodes_pids[i],io_index,isclone,iomodel,analysis,nodes_approx[i]); + if(finite_element==MINIcondensedEnum || finite_element==P1bubblecondensedEnum){ + /*Bubble function is collapsed, needs to constrain it, maybe this is not the best place to do this, but that's life!*/ + if(i>=iomodel->numberofvertices && inumberofvertices+iomodel->numberofelements){ node->HardDeactivate(); - nodes->AddObject(node); - } - } - /*P1 pressure*/ - vnodes = id0+iomodel->numberofvertices+iomodel->numberofelements; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum)); - } - } - break; - case MINIEnum: - _assert_(approximation==FSApproximationEnum); - /*P1+ velocity*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - /*P1 pressure*/ - vnodes = id0+iomodel->numberofvertices+iomodel->numberofelements; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum)); - } - } - break; - case TaylorHoodEnum: - case XTaylorHoodEnum: - _assert_(approximation==FSApproximationEnum); - /*P2 velocity*/ - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; - if(iomodel->meshelementtype==PentaEnum){ - FacesPartitioning(iomodel); - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - } - - /*P1 pressure*/ - if(iomodel->meshelementtype==PentaEnum){ - numberoffaces=iomodel->numberoffaces; - } - else{ - numberoffaces=0; - } - vnodes = id0+numberoffaces; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofedges+numberoffaces+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum)); - } - } - break; - case LATaylorHoodEnum: - _assert_(approximation==FSApproximationEnum); - /*P2 velocity*/ - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; - if(iomodel->meshelementtype==PentaEnum){ - FacesPartitioning(iomodel); - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - } - - /*No pressure*/ - break; - case OneLayerP4zEnum: - _assert_(approximation==FSApproximationEnum); - EdgesPartitioning(iomodel); - FacesPartitioning(iomodel); - /*P2xP4 velocity*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - counter = iomodel->numberofvertices; - for(i=0;inumberofedges;i++){ - if(iomodel->edges[i*3+2]==2){/*Vertical edges*/ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - counter=counter+3; - } - else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/ - if(iomodel->my_edges[i]){ - node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - counter=counter+1; - } - else{ - _error_("not supported"); - } - } - id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges; - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - node = new Node(id0+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - node = new Node(id0+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - counter=counter+3; - } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - - /*P1 pressure*/ - vnodes = id0+3*iomodel->numberoffaces; - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(vnodes+i+1,counter+1,lid++,0,i,false,iomodel,analysis,FSpressureEnum)); - } - counter++; - } - break; - case CrouzeixRaviartEnum: - _assert_(approximation==FSApproximationEnum); - /*P2b velocity*/ - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; - if(iomodel->meshelementtype==PentaEnum){ - FacesPartitioning(iomodel); - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - id0 = id0+iomodel->numberoffaces; - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - nodes->AddObject(new Node(id0+i+1,id0+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - - /*P1 DG pressure*/ - vnodes = id0+iomodel->numberofelements; - switch(iomodel->meshelementtype){ - case TriaEnum: elementnbv = 3; break; - case TetraEnum: elementnbv = 4; break; - case PentaEnum: elementnbv = 6; break; - default: _error_("mesh dimension not supported yet"); - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - for(j=0;jAddObject(new Node(vnodes+elementnbv*i+j+1,vnodes+elementnbv*i+j,lid++,0,iomodel->elements[+elementnbv*i+j]-1,false,iomodel,analysis,FSpressureEnum)); - - } - } - } - break; - case LACrouzeixRaviartEnum: - _assert_(approximation==FSApproximationEnum); - /*P2b velocity*/ - EdgesPartitioning(iomodel); - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum)); - } - } - for(i=0;inumberofedges;i++){ - if(iomodel->my_edges[i]){ - nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; - if(iomodel->meshelementtype==PentaEnum){ - FacesPartitioning(iomodel); - for(i=0;inumberoffaces;i++){ - if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ - if(iomodel->my_faces[i]){ - node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum); - nodes->AddObject(node); - } - } - else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/ - /*Nothing*/ - } - else{ - _error_("not supported"); - } - } - id0 = id0+iomodel->numberoffaces; - } - for(i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - nodes->AddObject(new Node(id0+i+1,id0+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum)); - } - } - - /*No pressure*/ - break; - - default: - _error_("Finite element "<(my_nodes); + } + } + nodes->AddObject(node); + } + } + xDelete(nodes_approx); + xDelete(nodes_ranks); + xDelete(nodes_lids); + xDelete(nodes_pids); + + nodes->common_send=common_send; + nodes->common_recv=common_recv; + nodes->common_send_ids=common_send_ids; + nodes->common_recv_ids=common_recv_ids; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 23575) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 23576) @@ -22,5 +22,5 @@ void UpdateElementsAndMaterialsDakota(Elements* elements,Materials* materials, IoModel* iomodel); void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel); -void CreateNodes(Nodes*nodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum); +void CreateNodes(Nodes*nodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum,int* approximations=NULL); /*partitioning: */ Index: /issm/trunk-jpl/src/c/modules/NodesDofx/NodesDofx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/NodesDofx/NodesDofx.cpp (revision 23575) +++ /issm/trunk-jpl/src/c/modules/NodesDofx/NodesDofx.cpp (revision 23576) @@ -20,5 +20,5 @@ /*Ensure that only for each cpu, the partition border nodes only will be taken into account once * across the cluster. To do so, we flag all the clone nodes: */ - nodes->FlagClones(configuration_type); + //nodes->FlagClones(configuration_type); /*Not needed anymore!*/ /*Go through all nodes, and build degree of freedom lists. Each node gets a fixed number of dofs. When