Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17539) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17540) @@ -234,5 +234,4 @@ iomodel->FetchDataToInput(elements,PressureEnum,0.); if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum); - if(finiteelement==XTaylorHoodEnum) this->InitializeXTH(elements); } if(islevelset){ @@ -821,53 +820,4 @@ } }/*}}}*/ -void StressbalanceAnalysis::InitializeXTH(Elements* elements){/*{{{*/ - - /*Intermediaries*/ - int meshtype,dim; - IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble viscosity; - IssmDouble *xyz_list = NULL; - - for(int i=0;iSize();i++){ - Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - - /*Get problem dimension*/ - element->FindParam(&meshtype,MeshTypeEnum); - switch(meshtype){ - case Mesh2DverticalEnum: dim = 2; break; - case Mesh3DEnum: dim = 3; break; - case Mesh3DtetrasEnum: dim = 3; break; - default: _error_("mesh "<GetNumberOfVertices(); //Tensors, P1 DG - - /*Get inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input; - if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} - - /*Get d and tau*/ - Gauss* gauss = element->NewGauss(); - for(int i=0;iGaussNode(P1DGEnum,i); - - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - if(dim==3){ - vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); - } - - _error_("STOP"); - element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); - } - - /*Clean up*/ - xDelete(xyz_list); - } - -}/*}}}*/ /*Finite Element Analysis*/ @@ -4109,4 +4059,88 @@ xDelete(vvalues); }/*}}}*/ +void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/ + + /*Intermediaries*/ + int meshtype,dim; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble viscosity; + IssmDouble *xyz_list = NULL; + + /*Get problem dimension*/ + parameters->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + case Mesh3DtetrasEnum: dim = 3; break; + default: _error_("mesh "<Size();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + + /*Get inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input; + if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);} + + /*Allocate new inputs*/ + int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG + IssmDouble* epsxx = xNew(tnumnodes); IssmDouble* sigmapxx = xNew(tnumnodes); + IssmDouble* epsyy = xNew(tnumnodes); IssmDouble* sigmapyy = xNew(tnumnodes); + IssmDouble* epsxy = xNew(tnumnodes); IssmDouble* sigmapxy = xNew(tnumnodes); + IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL; + IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL; + IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL; + if(dim==3){ + epszz = xNew(tnumnodes); sigmapzz = xNew(tnumnodes); + epsxz = xNew(tnumnodes); sigmapxz = xNew(tnumnodes); + epsyz = xNew(tnumnodes); sigmapyz = xNew(tnumnodes); + } + + /*Get d and tau*/ + Gauss* gauss = element->NewGauss(); + for(int i=0;iGaussNode(P1DGEnum,i); + + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + if(dim==3){ + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + } + + element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + epsxx[i] = dvx[0]; sigmapxx[i] = 2.*viscosity*epsxx[i]; + epsyy[i] = dvy[1]; sigmapyy[i] = 2.*viscosity*epsyy[i]; + epsxy[i] = 0.5*(dvx[1] + dvy[0]); sigmapxy[i] = 2.*viscosity*epsxy[i]; + if(dim==3){ + epszz[i] = dvz[2]; sigmapzz[i] = 2.*viscosity*epszz[i]; + epsxz[i] = 0.5*(dvx[2] + dvz[0]); sigmapxz[i] = 2.*viscosity*epsxz[i]; + epsyz[i] = 0.5*(dvy[2] + dvz[1]); sigmapyz[i] = 2.*viscosity*epsyz[i]; + } + } + + /*Add inputs*/ + element->AddInput(StrainRatexxEnum,epsxx,P1DGEnum); element->AddInput(DeviatoricStressxxEnum,sigmapxx,P1DGEnum); + element->AddInput(StrainRateyyEnum,epsyy,P1DGEnum); element->AddInput(DeviatoricStressyyEnum,sigmapyy,P1DGEnum); + element->AddInput(StrainRatexyEnum,epsxy,P1DGEnum); element->AddInput(DeviatoricStressxyEnum,sigmapxy,P1DGEnum); + if(dim==3){ + element->AddInput(StrainRatezzEnum,epszz,P1DGEnum); element->AddInput(DeviatoricStresszzEnum,sigmapzz,P1DGEnum); + element->AddInput(StrainRatexzEnum,epsxz,P1DGEnum); element->AddInput(DeviatoricStressxzEnum,sigmapxz,P1DGEnum); + element->AddInput(StrainRateyzEnum,epsyz,P1DGEnum); element->AddInput(DeviatoricStressyzEnum,sigmapyz,P1DGEnum); + } + + /*Clean up*/ + xDelete(xyz_list); + xDelete(epsxx); + xDelete(epsyy); + xDelete(epszz); + xDelete(epsxy); + xDelete(epsxz); + xDelete(epsyz); + } + +}/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17539) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17540) @@ -19,5 +19,4 @@ void CreateConstraints(Constraints* constraints,IoModel* iomodel); void CreateLoads(Loads* loads, IoModel* iomodel); - void InitializeXTH(Elements* elements); /*Finite element Analysis*/ @@ -82,4 +81,5 @@ void GetSolutionFromInputsFS(Vector* solution,Element* element); void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element); + void InitializeXTH(Elements* elements,Parameters* parameters); /*Coupling*/ ElementMatrix* CreateKMatrixSSA3d(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 17539) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 17540) @@ -246,4 +246,9 @@ this->element_type = TaylorHoodEnum; return; + case XTaylorHoodEnum: + this->element_type = P2Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = XTaylorHoodEnum; + return; default: _error_("Element type "<element_type)<<" not supported yet"); @@ -280,4 +285,9 @@ this->GetNodalFunctions(basis,gauss); this->element_type = TaylorHoodEnum; + return; + case XTaylorHoodEnum: + this->element_type = P1Enum; + this->GetNodalFunctions(basis,gauss); + this->element_type = XTaylorHoodEnum; return; default: @@ -439,5 +449,5 @@ this->element_type = P2Enum; this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); - this->element_type = TaylorHoodEnum; + this->element_type = XTaylorHoodEnum; return; default: Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 17539) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 17540) @@ -20,9 +20,16 @@ int configuration_type; + /*Create analysis*/ + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + /*Recover parameters: */ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*Update constraints and initialize d and tau if necessary*/ femmodel->UpdateConstraintsx(); + analysis->InitializeXTH(femmodel->elements,femmodel->parameters); + + /*Solve KU=F*/ SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); - _error_("STOP"); CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); Reduceloadx(pf, Kfs, ys); delete Kfs; @@ -30,5 +37,12 @@ delete Kff; delete pf; delete df; Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; + + /*Update solution*/ InputUpdateFromSolutionx(femmodel,ug); + + /*Update d and tau accordingly*/ + _error_("STOP"); + delete ug; + delete analysis; }