Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18172) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18173) @@ -414,11 +414,12 @@ /*Deduce velocity interpolation from finite element*/ switch(finiteelement){ - case P1P1Enum : finiteelement = P1Enum; break; - case P1P1GLSEnum : finiteelement = P1Enum; break; - case MINIcondensedEnum : finiteelement = P1bubbleEnum; break; - case MINIEnum : finiteelement = P1bubbleEnum; break; - case TaylorHoodEnum : finiteelement = P2Enum; break; - case XTaylorHoodEnum : finiteelement = P2Enum; break; - case OneLayerP4zEnum : finiteelement = P2xP4Enum; break; + case P1P1Enum : finiteelement = P1Enum; break; + case P1P1GLSEnum : finiteelement = P1Enum; break; + case MINIcondensedEnum : finiteelement = P1bubbleEnum; break; + case MINIEnum : finiteelement = P1bubbleEnum; break; + case TaylorHoodEnum : finiteelement = P2Enum; break; + case XTaylorHoodEnum : finiteelement = P2Enum; break; + case OneLayerP4zEnum : finiteelement = P2xP4Enum; break; + case CrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break; default: _error_("finite element "<(vertices_type[i])==NoneApproximationEnum){ constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; + case P2bubbleEnum: + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if(reCast(vertices_type[i])==NoneApproximationEnum){ + constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+iomodel->numberoffaces+i+1,0,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); count++; } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18172) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18173) @@ -2362,4 +2362,34 @@ penta_node_ids[35]=iomodel->nodecounter+iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+5]; break; + case CrouzeixRaviartEnum: + numnodes = 25; + penta_node_ids = xNew(numnodes); + penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0]; + penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1]; + penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2]; + penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3]; + penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4]; + penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5]; + penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; + penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; + penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; + penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; + penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; + penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; + penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; + penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; + penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; + penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1; + penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1; + penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1; + penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1; + + penta_node_ids[19]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+1; + penta_node_ids[20]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+2; + penta_node_ids[21]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+3; + penta_node_ids[22]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+4; + penta_node_ids[23]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+5; + penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+6; + break; default: _error_("Finite element "<(numindices); + indices[0] = 0; + indices[1] = 1; + indices[2] = 2; + indices[3] = 9; + indices[4] = 10; + indices[5] = 11; + break; + case P2bubbleEnum: numindices = 6; indices = xNew(numindices); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18172) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18173) @@ -1975,4 +1975,19 @@ tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+2]; break; + case CrouzeixRaviartEnum: + numnodes = 10; + tria_node_ids = xNew(numnodes); + tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0]; + tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1]; + tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2]; + tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1; + tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1; + tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1; + tria_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+index+1; + + tria_node_ids[7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+1; + tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+2; + tria_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+3; + break; default: _error_("Finite element "<nodecounter; bool *my_faces = NULL; @@ -442,4 +442,60 @@ } break; + case CrouzeixRaviartEnum: + _assert_(approximation==FSApproximationEnum); + /*P2b velocity*/ + EdgesPartitioning(&my_edges,iomodel); + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum)); + } + } + for(i=0;inumberofedges;i++){ + if(my_edges[i]){ + nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum)); + } + } + id0 = id0+iomodel->numberofvertices+iomodel->numberofedges; + if(iomodel->meshelementtype==PentaEnum){ + FacesPartitioning(&my_faces,iomodel); + for(i=0;inumberoffaces;i++){ + if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/ + if(my_faces[i]){ + node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,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-iomodel->nodecounter+i,lid++,0,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-iomodel->nodecounter+elementnbv*i+j,lid++,iomodel->elements[+elementnbv*i+j]-1,iomodel,analysis,FSpressureEnum)); + + } + } + } + break; default: Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18172) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18173) @@ -615,4 +615,5 @@ XTaylorHoodEnum, OneLayerP4zEnum, + CrouzeixRaviartEnum, /*}}}*/ /*Results{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18172) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18173) @@ -604,4 +604,5 @@ case XTaylorHoodEnum : return "XTaylorHood"; case OneLayerP4zEnum : return "OneLayerP4z"; + case CrouzeixRaviartEnum : return "CrouzeixRaviart"; case SaveResultsEnum : return "SaveResults"; case BoolExternalResultEnum : return "BoolExternalResult"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18172) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18173) @@ -616,4 +616,5 @@ else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; + else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; - else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum; + if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; + else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum; else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum; else if (strcmp(name,"MisfitName")==0) return MisfitNameEnum; Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18172) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18173) @@ -134,5 +134,5 @@ md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'}); md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'}); - md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'}); + md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'}); md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.); md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5); Index: /issm/trunk-jpl/src/m/classes/flowequation.py =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18172) +++ /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18173) @@ -78,5 +78,5 @@ md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble']) md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4']) - md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z']) + md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart']) md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1]) md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])