Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17524) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17525) @@ -355,4 +355,5 @@ ./solutionsequences/solutionsequence_nonlinear.cpp\ ./solutionsequences/solutionsequence_newton.cpp\ + ./solutionsequences/solutionsequence_la_theta.cpp\ ./solutionsequences/convergence.cpp\ ./classes/Options/Options.h\ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17524) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17525) @@ -407,6 +407,7 @@ case MINIEnum : finiteelement = P1bubbleEnum; break; case TaylorHoodEnum : finiteelement = P2Enum; break; + case XTaylorHoodEnum : finiteelement = P2Enum; break; case OneLayerP4zEnum : finiteelement = P2xP4Enum; break; - default: _error_("finite element "<parameters->FindParam(&isHO,FlowequationIsHOEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&fe_FS,FlowequationFeFSEnum); femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum); - if((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor + if(isFS){ + if(VerboseSolution()) _printf0_(" computing velocities\n"); + + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + if (fe_FS==XTaylorHoodEnum) + solutionsequence_la_theta(femmodel); + else if(newton>0) + solutionsequence_newton(femmodel); + else + solutionsequence_nonlinear(femmodel,conserve_loads); + } + else if(isSSA || isHO || isL1L2){ if(VerboseSolution()) _printf0_(" computing velocities\n"); @@ -851,4 +864,7 @@ extrudefrombase_core(femmodel); } + } + else{ + _error_("not supported"); } @@ -2857,6 +2873,15 @@ ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/ + /*Get type of algorithm*/ + int fe_FS; + element->FindParam(&fe_FS,FlowequationFeFSEnum); + /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixFSViscous(element); + ElementMatrix* Ke1=NULL; + if(fe_FS==XTaylorHoodEnum) + Ke1=CreateKMatrixFSViscousXTH(element); + else + Ke1=CreateKMatrixFSViscous(element); + ElementMatrix* Ke2=CreateKMatrixFSFriction(element); ElementMatrix* Ke3=CreateKMatrixFSShelf(element); @@ -2867,4 +2892,79 @@ delete Ke2; delete Ke3; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousXTH(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,meshtype,dim,epssize; + IssmDouble r,FSreconditioning,Jdet; + IssmDouble *xyz_list = NULL; + + /*FIXME this should change*/ + r = 1.; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; epssize = 3; break; + case Mesh3DEnum: dim = 3; epssize = 6; break; + case Mesh3DtetrasEnum: dim = 3; epssize = 6; break; + default: _error_("mesh "<NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int numdof = vnumnodes*dim + pnumnodes; + int bsize = epssize + 2; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementMatrix(FSvelocityEnum); + IssmDouble* B = xNew(bsize*numdof); + IssmDouble* Bprime = xNew(bsize*numdof); + IssmDouble* D = xNewZeroInit(bsize*bsize); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum); + 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);} + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + this->GetBFS(B,element,dim,xyz_list,gauss); + this->GetBFSprime(Bprime,element,dim,xyz_list,gauss); + + for(i=0;iweight*Jdet; + for(i=epssize;iweight*Jdet; + + TripleMultiply(B,bsize,numdof,1, + D,bsize,bsize,0, + Bprime,bsize,numdof,0, + &Ke->values[0],1); + } + + /*Transform Coordinate System*/ + element->TransformStiffnessMatrixCoord(Ke,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(D); + xDelete(Bprime); + xDelete(B); + xDelete(cs_list); return Ke; }/*}}}*/ @@ -3179,19 +3279,113 @@ ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/ - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorFSViscous(element); - ElementVector* pe2=CreatePVectorFSShelf(element); - ElementVector* pe3=CreatePVectorFSFront(element); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); + ElementVector* pe = NULL; + + int fe_FS; + element->FindParam(&fe_FS,FlowequationFeFSEnum); + + if(fe_FS==XTaylorHoodEnum){ + ElementVector* pe1=CreatePVectorFSViscous(element); + ElementVector* pe2=CreatePVectorFSShelf(element); + ElementVector* pe3=CreatePVectorFSFront(element); + ElementVector* petemp =new ElementVector(pe1,pe2,pe3); + ElementVector* pe4=CreatePVectorFSViscousXTH(element); + ElementVector* pe = new ElementVector(petemp,pe4); + delete pe1; + delete pe2; + delete pe3; + delete petemp; + delete pe4; + } + else{ + ElementVector* pe1=CreatePVectorFSViscous(element); + ElementVector* pe2=CreatePVectorFSShelf(element); + ElementVector* pe3=CreatePVectorFSFront(element); + pe =new ElementVector(pe1,pe2,pe3); + delete pe1; + delete pe2; + delete pe3; + } /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; }/*}}}*/ #endif ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/ + + int i,meshtype,dim; + IssmDouble Jdet,forcex,forcey,forcez; + IssmDouble *xyz_list = NULL; + + /*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 "<NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); + Input* loadingforcex_input=element->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); + Input* loadingforcey_input=element->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); + Input* loadingforcez_input=NULL; + if(dim==3){ + loadingforcez_input=element->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + + loadingforcex_input->GetInputValue(&forcex,gauss); + loadingforcey_input->GetInputValue(&forcey,gauss); + if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss); + + for(i=0;ivalues[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i]; + pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i]; + if(dim==3){ + pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i]; + pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + } + else{ + pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i]; + } + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + xDelete(xyz_list); + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/ + _error_("STOP"); int i,meshtype,dim; Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17524) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17525) @@ -67,4 +67,5 @@ ElementMatrix* CreateJacobianMatrixFS(Element* element); ElementMatrix* CreateKMatrixFS(Element* element); + ElementMatrix* CreateKMatrixFSViscousXTH(Element* element); ElementMatrix* CreateKMatrixFSViscous(Element* element); ElementMatrix* CreateKMatrixFSFriction(Element* element); @@ -72,4 +73,5 @@ ElementVector* CreatePVectorFS(Element* element); ElementVector* CreatePVectorFSViscous(Element* element); + ElementVector* CreatePVectorFSViscousXTH(Element* element); ElementVector* CreatePVectorFSShelf(Element* element); ElementVector* CreatePVectorFSFront(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17524) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17525) @@ -1596,5 +1596,5 @@ void Tria::ResetFSBasalBoundaryCondition(void){ - int numnodes = this->GetNumberOfNodes(); + int numnodes = this->NumberofNodesVelocity(); int approximation; @@ -1919,4 +1919,5 @@ break; case TaylorHoodEnum: + case XTaylorHoodEnum: numnodes = 9; tria_node_ids = xNew(numnodes); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 17524) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 17525) @@ -421,4 +421,5 @@ return; case TaylorHoodEnum: + case XTaylorHoodEnum: this->element_type = P2Enum; this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss); @@ -588,4 +589,5 @@ case MINIEnum: return NUMNODESP1b+NUMNODESP1; case TaylorHoodEnum: return NUMNODESP2+NUMNODESP1; + case XTaylorHoodEnum: return NUMNODESP2+NUMNODESP1; default: _error_("Element type "<element_type)<<" not supported yet"); } @@ -603,4 +605,5 @@ case MINIEnum: return NUMNODESP1; case TaylorHoodEnum: return NUMNODESP1; + case XTaylorHoodEnum: return NUMNODESP1; default: _error_("Element type "<element_type)<<" not supported yet"); } @@ -618,4 +621,5 @@ case MINIEnum: return NUMNODESP1b; case TaylorHoodEnum: return NUMNODESP2; + case XTaylorHoodEnum: return NUMNODESP2; default: _error_("Element type "<element_type)<<" not supported yet"); } @@ -633,4 +637,5 @@ case MINIEnum: return P1bubbleEnum; case TaylorHoodEnum: return P2Enum; + case XTaylorHoodEnum: return P2Enum; default: _error_("Element type "<element_type)<<" not supported yet"); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17524) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 17525) @@ -294,4 +294,5 @@ break; case TaylorHoodEnum: + case XTaylorHoodEnum: _assert_(approximation==FSApproximationEnum); /*P2 velocity*/ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 17525) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (revision 17525) @@ -0,0 +1,34 @@ +/*!\file: solutionsequence_la_theta.cpp + * \brief: numerical core of la_theta solutions + */ + +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +void solutionsequence_la_theta(FemModel* femmodel){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->UpdateConstraintsx(); + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + _error_("STOP"); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys); delete Kfs; + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); + delete Kff; delete pf; delete df; + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; + InputUpdateFromSolutionx(femmodel,ug); + delete ug; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17524) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17525) @@ -18,4 +18,5 @@ void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads); void solutionsequence_linear(FemModel* femmodel); +void solutionsequence_la_theta(FemModel* femmodel); void solutionsequence_adjoint_linear(FemModel* femmodel);