Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15620) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15621) @@ -3095,16 +3095,15 @@ /*Intermediaries*/ int i; - int penta_node_ids[6]; int penta_vertex_ids[6]; IssmDouble nodeinputs[6]; IssmDouble yts; - int stabilization; bool dakota_analysis; bool isFS; IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat; + int numnodes; + int* penta_node_ids = NULL; /*Fetch parameters: */ iomodel->Constant(&yts,ConstantsYtsEnum); - iomodel->Constant(&stabilization,PrognosticStabilizationEnum); iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); iomodel->Constant(&isFS,FlowequationIsFSEnum); @@ -3120,4 +3119,5 @@ /*}}}*/ + /*Recover element type*/ this->SetElementType(finiteelement_type,analysis_counter); @@ -3126,12 +3126,41 @@ /*Recover nodes ids needed to initialize the node hook.*/ - for(i=0;i<6;i++){ - //go recover node ids, needed to initialize the node hook. - //WARNING: We assume P1 elements here!!!!! - penta_node_ids[i]=iomodel->nodecounter+iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab + switch(finiteelement_type){ + case P1Enum: + numnodes = 6; + 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]; + break; + case P2Enum: + numnodes = 15; + 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; + break; + default: + _error_("Finite element "<SetHookNodes(penta_node_ids,6,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + this->SetHookNodes(penta_node_ids,numnodes,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + xDelete(penta_node_ids); /*Fill with IoModel*/ @@ -6843,23 +6872,23 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticHOViscous(void){ - /*Constants*/ - const int numdof=NDOF2*NUMVERTICES; - /*Intermediaries */ - int i,j; - int approximation; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble Jdet; - IssmDouble viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity - IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble D_scalar; - IssmDouble D[5][5]={0.0}; // material matrix, simple scalar matrix. - IssmDouble B[5][numdof]; - IssmDouble Bprime[5][numdof]; - Tria* tria=NULL; + int i,j; + int approximation; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdet; + IssmDouble viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity + IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble D_scalar; GaussPenta *gauss=NULL; - /*Initialize Element matrix*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum); + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,HOApproximationEnum); + IssmDouble* B = xNew(5*numdof); + IssmDouble* Bprime = xNew(5*numdof); + IssmDouble* D = xNewZeroInit(5*5); /*Retrieve all inputs and parameters*/ @@ -6879,6 +6908,6 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - GetBHO(&B[0][0], &xyz_list[0][0], gauss); - GetBprimeHO(&Bprime[0][0], &xyz_list[0][0], gauss); + GetBHO(&B[0], &xyz_list[0][0], gauss); + GetBprimeHO(&Bprime[0], &xyz_list[0][0], gauss); this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); @@ -6889,17 +6918,20 @@ D_scalar=2*newviscosity*gauss->weight*Jdet; - for (i=0;i<5;i++) D[i][i]=D_scalar; - - TripleMultiply( &B[0][0],5,numdof,1, - &D[0][0],5,5,0, - &Bprime[0][0],5,numdof,0, + for (i=0;i<5;i++) D[i*5+i]=D_scalar; + + TripleMultiply(B,5,numdof,1, + D,5,5,0, + Bprime,5,numdof,0, &Ke->values[0],1); } /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); + TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); /*Clean up and return*/ delete gauss; + xDelete(D); + xDelete(Bprime); + xDelete(B); return Ke; } @@ -6908,24 +6940,26 @@ ElementMatrix* Penta::CreateKMatrixDiagnosticHOFriction(void){ - /*Constants*/ - const int numdof = NDOF2*NUMVERTICES; - /*Intermediaries */ - int i,j; - int analysis_type,migration_style; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; - IssmDouble alpha2,Jdet; - IssmDouble phi=1.0; - IssmDouble L[2][numdof]; - IssmDouble DL[2][2]={{ 0,0 },{0,0}}; //for basal drag - IssmDouble DL_scalar; - Friction *friction = NULL; - GaussPenta *gauss=NULL; + int i,j; + int analysis_type,migration_style; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.}; + IssmDouble alpha2,Jdet; + IssmDouble phi=1.0; + IssmDouble DL_scalar; + Friction *friction = NULL; + GaussPenta *gauss = NULL; /*Initialize Element matrix and return if necessary*/ if(IsFloating() || !IsOnBed()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum); + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,HOApproximationEnum); + IssmDouble* B = xNew(2*numdof); + IssmDouble* D = xNewZeroInit(2*2); /*Retrieve all inputs and parameters*/ @@ -6951,6 +6985,6 @@ gauss->GaussPoint(ig); - GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss); - GetBHOFriction(&L[0][0],gauss); + GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss); + GetBHOFriction(&B[0],gauss); friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); @@ -6958,17 +6992,19 @@ DL_scalar=alpha2*gauss->weight*Jdet; - for (i=0;i<2;i++) DL[i][i]=DL_scalar; - - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, + for (i=0;i<2;i++) D[i*2+i]=DL_scalar; + + TripleMultiply(B,2,numdof,1, + D,2,2,0, + B,2,numdof,0, &Ke->values[0],1); } /*Transform Coordinate System*/ - TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum); + TransformStiffnessMatrixCoord(Ke,nodes,numnodes,XYEnum); /*Clean up and return*/ delete gauss; + xDelete(D); + xDelete(B); delete friction; return Ke; Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15620) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15621) @@ -45,6 +45,4 @@ /*FUNCTION PentaRef::SetElementType{{{*/ void PentaRef::SetElementType(int type,int type_counter){ - - _assert_(type==P1Enum || type==P1DGEnum); /*initialize element type*/ @@ -143,68 +141,71 @@ * where h is the interpolation function for node i. * - * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) - */ - - IssmDouble dbasis[3][NUMNODESP1]; - - /*Get dbasis in actual coordinate system: */ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss); + * We assume B has been allocated already, of size: 5x(NDOF2*numnodes) + */ + + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(3*numnodes); + GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build B: */ - for (int i=0;i(dbasis); } /*}}}*/ /*FUNCTION PentaRef::GetBprimeHO {{{*/ -void PentaRef::GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){ +void PentaRef::GetBprimeHO(IssmDouble* B,IssmDouble* xyz_list,GaussPenta* gauss){ /*Compute B prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. * For node i, Bi can be expressed in the actual coordinate system * by: * Bi=[ 2*dh/dx dh/dy ] - * [ dh/dx 2*dh/dy ] - * [ dh/dy dh/dx ] - * [ dh/dz 0 ] - * [ 0 dh/dz ] + * [ dh/dx 2*dh/dy ] + * [ dh/dy dh/dx ] + * [ dh/dz 0 ] + * [ 0 dh/dz ] * where h is the interpolation function for node i. * - * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1) - */ - IssmDouble dbasis[3][NUMNODESP1]; - - /*Get dbasis in actual coordinate system: */ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord); + * We assume B has been allocated already, of size: 5x(NDOF2*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(3*numnodes); + GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); /*Build BPrime: */ - for(int i=0;i(dbasis); } /*}}}*/ @@ -1183,4 +1184,7 @@ dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i]; } + + /*Clean up*/ + xDelete(dbasis_ref); } @@ -1622,4 +1626,5 @@ switch(this->element_type){ case P1Enum: return NUMNODESP1; + case P2Enum: return NUMNODESP2; case MINIEnum: return NUMNODESP1b; default: _error_("Element type "<element_type)<<" not supported yet"); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15620) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15621) @@ -2320,5 +2320,5 @@ /*}}}*/ /*FUNCTION Tria::Update{{{*/ -void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement_type){ +void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ /*Intermediaries*/ @@ -2346,5 +2346,5 @@ /*Recover element type*/ - this->SetElementType(finitelement_type,analysis_counter); + this->SetElementType(finiteelement_type,analysis_counter); /*Recover vertices ids needed to initialize inputs*/ @@ -2354,5 +2354,5 @@ /*Recover nodes ids needed to initialize the node hook.*/ - switch(finitelement_type){ + switch(finiteelement_type){ case P1Enum: numnodes = 3; @@ -2380,5 +2380,5 @@ break; default: - _error_("Finite element "<my_vertices[i])){ if (!xIsNan(spcdata[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcdata[i],analysis_type)); + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,spcdata[i],analysis_type)); count++; } @@ -61,5 +61,5 @@ if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 1,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); count++; } @@ -90,5 +90,5 @@ if(spcpresent){ - constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,analysis_type)); + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,dof,N,times,values,analysis_type)); count++; } @@ -108,5 +108,5 @@ } if(spcpresent){ - constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1, + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,dof, N,times,values,analysis_type)); count++; Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15620) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15621) @@ -8,5 +8,5 @@ /* local prototypes: */ -void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element); +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element,int dof=1); #endif /* _IOMODELTOELEMENTINPUTX_H */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp (revision 15620) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp (revision 15621) @@ -82,5 +82,5 @@ if(edgestemp[e*2+1]==v2+1){ exist = true; - element_edge_connectivity[i*elementnbv+j]=e; + element_edge_connectivity[i*elementnbe+j]=e; break; } @@ -96,5 +96,5 @@ /*Update Connectivity*/ - element_edge_connectivity[i*elementnbv+j]=nbe; + element_edge_connectivity[i*elementnbe+j]=nbe; /*Update chain*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15620) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15621) @@ -13,10 +13,9 @@ int i,j; int count; - IssmDouble yts; IssmDouble g; IssmDouble rho_ice; IssmDouble FSreconditioning; bool isSSA,isL1L2,isHO,isFS; - int fe_ssa; + int fe_ssa,fe_ho; bool spcpresent = false; int Mx,Nx; @@ -44,5 +43,4 @@ /*Fetch parameters: */ - iomodel->Constant(&yts,ConstantsYtsEnum); iomodel->Constant(&g,ConstantsGEnum); iomodel->Constant(&rho_ice,MaterialsRhoIceEnum); @@ -53,4 +51,5 @@ iomodel->Constant(&isFS,FlowequationIsFSEnum); iomodel->Constant(&fe_ssa,FlowequationFeSSAEnum); + iomodel->Constant(&fe_ho,FlowequationFeHOEnum); /*Recover pointer: */ @@ -108,9 +107,9 @@ count++; if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -123,9 +122,9 @@ count++; if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -145,9 +144,9 @@ count++; if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -160,13 +159,13 @@ count++; if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } if (!xIsNan(spcvz[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -185,9 +184,9 @@ count++; if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -200,13 +199,13 @@ count++; if (!xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - if (!xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + if (!xIsNan(spcvy[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,4,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } if (!xIsNan(spcvz[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,5,spcvz[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; } @@ -217,5 +216,5 @@ else{ if (Mx==iomodel->numberofvertices && !xIsNan(spcvx[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcvx[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; @@ -226,5 +225,5 @@ spcpresent=false; for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default } @@ -242,5 +241,5 @@ if (My==iomodel->numberofvertices && !xIsNan(spcvy[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy. + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,spcvy[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vy. count++; } @@ -250,5 +249,5 @@ spcpresent=false; for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default } @@ -266,5 +265,5 @@ if (reCast(vertices_type[i])==FSApproximationEnum || (reCast(vertices_type[i])==NoneApproximationEnum)){ if (Mz==iomodel->numberofvertices && !xIsNan(spcvz[i])){ - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i]/yts,DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,spcvz[i],DiagnosticHorizAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy count++; } @@ -274,5 +273,5 @@ spcpresent=false; for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default } @@ -313,5 +312,5 @@ /*SPC Quadratic elements*/ - if(isSSA && fe_ssa==1){ + if((isSSA && fe_ssa==1) || (isHO && fe_ho==1)){ int v1,v2; @@ -330,10 +329,10 @@ if(!xIsNan(spcvx[v1]) && !xIsNan(spcvx[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 1,(spcvx[v1]+spcvx[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); + 1,(spcvx[v1]+spcvx[v2])/(2.),DiagnosticHorizAnalysisEnum)); count++; } if(!xIsNan(spcvy[v1]) && !xIsNan(spcvy[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 2,(spcvy[v1]+spcvy[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); + 2,(spcvy[v1]+spcvy[v2])/(2.),DiagnosticHorizAnalysisEnum)); count++; } @@ -341,5 +340,5 @@ if(!xIsNan(spcvz[v1]) && !xIsNan(spcvz[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 3,(spcvz[v1]+spcvz[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); + 3,(spcvz[v1]+spcvz[v2])/(2.),DiagnosticHorizAnalysisEnum)); count++; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15620) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15621) @@ -49,5 +49,10 @@ else if(isHO){ approximation = HOApproximationEnum; - finiteelement = P1Enum; + iomodel->Constant(&temp,FlowequationFeHOEnum); + switch(temp){ + case 0 : finiteelement = P1Enum; break; + case 1 : finiteelement = P2Enum; break; + default: _error_("finite element "<Constant(&temp,FlowequationFeHOEnum); + switch(temp){ + case 0 : finiteelement = P1Enum; break; + case 1 : finiteelement = P2Enum; break; + default: _error_("finite element "<