Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 18178) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 18179) @@ -354,8 +354,9 @@ ./analyses/EnumToAnalysis.cpp\ ./analyses/Analysis.h\ + ./solutionsequences/solutionsequence_la.cpp\ + ./solutionsequences/solutionsequence_la_theta.cpp\ ./solutionsequences/solutionsequence_linear.cpp\ ./solutionsequences/solutionsequence_nonlinear.cpp\ ./solutionsequences/solutionsequence_newton.cpp\ - ./solutionsequences/solutionsequence_la_theta.cpp\ ./solutionsequences/convergence.cpp\ ./classes/Options/Options.h\ @@ -538,4 +539,7 @@ issm_sources += ./analyses/StressbalanceAnalysis.cpp endif +if UZAWAPRESSURE +issm_sources += ./analyses/UzawaPressureAnalysis.cpp +endif if STRESSBALANCESIA issm_sources += ./analyses/StressbalanceSIAAnalysis.cpp Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 18179) @@ -95,4 +95,7 @@ case ThermalAnalysisEnum : return new ThermalAnalysis(); #endif + #ifdef _HAVE_UZAWAPRESSURE_ + case UzawaPressureAnalysisEnum : return new UzawaPressureAnalysis(); + #endif #ifdef _HAVE_GIA_ case GiaAnalysisEnum : return new GiaAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18179) @@ -867,4 +867,6 @@ if (fe_FS==XTaylorHoodEnum) solutionsequence_la_theta(femmodel); + else if (fe_FS==LATaylorHoodEnum) + solutionsequence_la(femmodel); else if(newton>0) solutionsequence_newton(femmodel); @@ -2852,4 +2854,6 @@ if(fe_FS==XTaylorHoodEnum) Ke1=CreateKMatrixFSViscousXTH(element); + else if(fe_FS==LATaylorHoodEnum) + Ke1=CreateKMatrixFSViscousLATH(element); else Ke1=CreateKMatrixFSViscous(element); @@ -2863,4 +2867,95 @@ delete Ke2; delete Ke3; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLATH(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,dim,epssize; + IssmDouble r,Jdet,viscosity,DU; + IssmDouble *xyz_list = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + element->FindParam(&r,AugmentedLagrangianREnum); + if(dim==2) epssize = 3; + else epssize = 6; + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes; + if(dim==2) pnumnodes=3; + else pnumnodes=6; + //int pnumnodes = element->NumberofNodes(P1Enum); + int numdof = vnumnodes*dim; + int pnumdof = pnumnodes; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes); + if(dim==2) for(i=0;iNewElementMatrix(FSvelocityEnum); + IssmDouble* B = xNew(epssize*numdof); + IssmDouble* Bprime = xNew(epssize*numdof); + IssmDouble* BtBUzawa = xNew(numdof*pnumdof); + IssmDouble* BU = xNew(pnumdof); + IssmDouble* BprimeU = xNew(numdof); + IssmDouble* D = xNewZeroInit(epssize*epssize); + + /*Retrieve all 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);} + + /* 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); + + element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + for(i=0;iweight*Jdet; + + TripleMultiply(B,epssize,numdof,1, + D,epssize,epssize,0, + Bprime,epssize,numdof,0, + &Ke->values[0],1); + + this->GetBFSUzawa(BU,element,dim,xyz_list,gauss); + this->GetBFSprimeUzawa(BprimeU,element,dim,xyz_list,gauss); + + DU = gauss->weight*Jdet; + + TripleMultiply(BU,1,pnumdof,1, + &DU,1,1,0, + BprimeU,1,numdof,0, + BtBUzawa,1); + + } + + MatrixMultiply(BtBUzawa,pnumdof,numdof,1, + BtBUzawa,pnumdof,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(BprimeU); + xDelete(BU); + xDelete(BtBUzawa); + xDelete(cs_list); return Ke; }/*}}}*/ @@ -3305,4 +3400,17 @@ ElementVector* petemp =new ElementVector(pe1,pe2,pe3); ElementVector* pe4=CreatePVectorFSViscousXTH(element); + pe = new ElementVector(petemp,pe4); + delete pe1; + delete pe2; + delete pe3; + delete petemp; + delete pe4; + } + else if(fe_FS==LATaylorHoodEnum){ + ElementVector* pe1=CreatePVectorFSViscous(element); + ElementVector* pe2=CreatePVectorFSShelf(element); + ElementVector* pe3=CreatePVectorFSFront(element); + ElementVector* petemp =new ElementVector(pe1,pe2,pe3); + ElementVector* pe4=CreatePVectorFSViscousLATH(element); pe = new ElementVector(petemp,pe4); delete pe1; @@ -3565,4 +3673,60 @@ xDelete(vdbasis); xDelete(tbasis); + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLATH(Element* element){/*{{{*/ + + int i,dim; + IssmDouble Jdet,r,pressure; + IssmDouble *xyz_list = NULL; + Gauss* gauss = NULL; + + /*Get problem dimension*/ + element->FindParam(&dim,DomainDimensionEnum); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(numnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* dbasis = xNew(3*numnodes); + + /*Retrieve all inputs and parameters*/ + element->FindParam(&r,AugmentedLagrangianREnum); + element->GetVerticesCoordinates(&xyz_list); + + /*Get d and tau*/ + Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); + + gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + + pressure_input->GetInputValue(&pressure, gauss); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + for(i=0;ivalues[i*dim+0] += r*pressure*gauss->weight*Jdet*dbasis[0*numnodes+i]; + pe->values[i*dim+1] += r*pressure*gauss->weight*Jdet*dbasis[1*numnodes+i]; + if(dim==3){ + pe->values[i*dim+2]+=r*pressure*gauss->weight*Jdet*dbasis[2*numnodes+i]; + } + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(xyz_list); + xDelete(dbasis); return pe; }/*}}}*/ @@ -3955,4 +4119,62 @@ xDelete(pbasis); }/*}}}*/ +void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. + */ + + /*Fetch number of nodes for this finite element*/ + int pnumnodes; + if(dim==2) pnumnodes=3; + else pnumnodes=6; + //int pnumnodes = element->NumberofNodes(P1Enum); + + /*Get nodal functions derivatives*/ + IssmDouble* basis =xNew(pnumnodes); + element->NodalFunctionsP1(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bvi' = [ dphi/dx dphi/dy ] + * + * In 3d + * Bvi=[ dh/dx dh/dy dh/dz ] + * where phi is the finiteelement function for node i. + */ + + /*Fetch number of nodes for this finite element*/ + int vnumnodes = element->NumberofNodesVelocity(); + + /*Get nodal functions derivatives*/ + IssmDouble* vdbasis=xNew(dim*vnumnodes); + element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss); + + /*Build B_prime: */ + if(dim==2){ + for(int i=0;i(vdbasis); +}/*}}}*/ void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18178) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18179) @@ -68,4 +68,5 @@ ElementMatrix* CreateJacobianMatrixFS(Element* element); ElementMatrix* CreateKMatrixFS(Element* element); + ElementMatrix* CreateKMatrixFSViscousLATH(Element* element); ElementMatrix* CreateKMatrixFSViscousXTH(Element* element); ElementMatrix* CreateKMatrixFSViscous(Element* element); @@ -74,4 +75,5 @@ ElementVector* CreatePVectorFS(Element* element); ElementVector* CreatePVectorFSViscous(Element* element); + ElementVector* CreatePVectorFSViscousLATH(Element* element); ElementVector* CreatePVectorFSViscousXTH(Element* element); ElementVector* CreatePVectorFSShelf(Element* element); @@ -80,4 +82,6 @@ void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); void GetSolutionFromInputsFS(Vector* solution,Element* element); void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp (revision 18179) +++ /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp (revision 18179) @@ -0,0 +1,189 @@ +#include "./UzawaPressureAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +int UzawaPressureAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 1; +}/*}}}*/ +void UzawaPressureAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum)); +}/*}}}*/ +void UzawaPressureAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Update elements: */ + int finiteelement = P1Enum; + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchDataToInput(elements,VzEnum); + iomodel->FetchDataToInput(elements,PressureEnum); +}/*}}}*/ +void UzawaPressureAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ + + int finiteelement = P1Enum; + ::CreateNodes(nodes,iomodel,UzawaPressureAnalysisEnum,finiteelement); +}/*}}}*/ +void UzawaPressureAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + return; +}/*}}}*/ +void UzawaPressureAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + return; +}/*}}}*/ + +/*Finite Element Analysis*/ +void UzawaPressureAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* UzawaPressureAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* UzawaPressureAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ +_error_("Not implemented"); +}/*}}}*/ +ElementMatrix* UzawaPressureAnalysis::CreateKMatrix(Element* element){/*{{{*/ + + /*Intermediaries*/ + IssmDouble D_scalar,Jdet; + IssmDouble *xyz_list = NULL; + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* M = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + Gauss* gauss = element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + this->GetM(M,element,gauss); + D_scalar=gauss->weight*Jdet; + TripleMultiply(M,1,numnodes,1, + &D_scalar,1,1,0, + M,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(M); + return Ke; +}/*}}}*/ +ElementVector* UzawaPressureAnalysis::CreatePVector(Element* element){/*{{{*/ + + /*Intermediaries*/ + int dim; + IssmDouble Jdet,r,divu; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble *xyz_list = NULL; + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element matrix and vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->FindParam(&dim,DomainDimensionEnum); + element->FindParam(&r,AugmentedLagrangianREnum); + 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);} + + Gauss* gauss = element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis, gauss); + 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); + } + + for(int i=0;ivalues[i]+=divu*r*Jdet*gauss->weight*basis[i]; + } + } + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); + xDelete(basis); + return pe; +}/*}}}*/ +void UzawaPressureAnalysis::GetM(IssmDouble* M,Element* element,Gauss* gauss){/*{{{*/ + /*Compute B matrix. M=[M1 M2 M3] */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void UzawaPressureAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void UzawaPressureAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ + _error_("Not implemented yet"); +}/*}}}*/ +void UzawaPressureAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + + int *doflist = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Fetch dof list and allocate solution vector*/ + element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numnodes); + IssmDouble* pressure = xNew(numnodes); + element->GetInputListOnNodes(&pressure[0],PressureEnum); + + for(int i=0;iAddInput(PressureEnum,values,element->GetElementType()); + + /*Free ressources:*/ + xDelete(values); + xDelete(pressure); + xDelete(doflist); + +}/*}}}*/ +void UzawaPressureAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + /*Default, do nothing*/ + return; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.h (revision 18179) +++ /issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.h (revision 18179) @@ -0,0 +1,34 @@ +/*! \file UzawaPressureAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _UzawaPressureAnalysis_ +#define _UzawaPressureAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class UzawaPressureAnalysis: public Analysis{ + + public: + /*Model processing*/ + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); + void CreateNodes(Nodes* nodes,IoModel* iomodel); + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetM(IssmDouble* M,Element* element,Gauss* gauss); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); +}; +#endif Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 18178) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 18179) @@ -37,4 +37,5 @@ #include "./StressbalanceSIAAnalysis.h" #include "./StressbalanceVerticalAnalysis.h" +#include "./UzawaPressureAnalysis.h" #include "./L2ProjectionBaseAnalysis.h" #include "./L2ProjectionEPLAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 18179) @@ -945,4 +945,5 @@ switch(finiteelement){ + case NoneEnum: return 0; case P0Enum: return NUMNODESP0; case P1Enum: return NUMNODESP1; @@ -962,4 +963,5 @@ case MINIEnum: return NUMNODESP1b+NUMNODESP1; case TaylorHoodEnum: return NUMNODESP2+NUMNODESP1; + case LATaylorHoodEnum: return NUMNODESP2; case OneLayerP4zEnum: return NUMNODESP2xP4+NUMNODESP1; case CrouzeixRaviartEnum: return NUMNODESP2b+NUMNODESP1; @@ -978,4 +980,5 @@ case MINIEnum: return P1bubbleEnum; case TaylorHoodEnum: return P2Enum; + case LATaylorHoodEnum: return P2Enum; case OneLayerP4zEnum: return P2xP4Enum; case CrouzeixRaviartEnum:return P2bubbleEnum; @@ -994,4 +997,5 @@ case MINIEnum: return P1Enum; case TaylorHoodEnum: return P1Enum; + case LATaylorHoodEnum: return P1Enum; case OneLayerP4zEnum: return P1Enum; case CrouzeixRaviartEnum:return P1DGEnum; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 18179) @@ -863,4 +863,18 @@ tetra_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[4*index+3]; break; + case LATaylorHoodEnum: + numnodes = 10; + tetra_node_ids = xNew(numnodes); + tetra_node_ids[0]=iomodel->nodecounter+iomodel->elements[4*index+0]; + tetra_node_ids[1]=iomodel->nodecounter+iomodel->elements[4*index+1]; + tetra_node_ids[2]=iomodel->nodecounter+iomodel->elements[4*index+2]; + tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3]; + tetra_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+0]+1; + tetra_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+1]+1; + tetra_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+2]+1; + tetra_node_ids[7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+3]+1; + tetra_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+4]+1; + tetra_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+5]+1; + break; default: _error_("Finite element "<nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+2]; break; + case LATaylorHoodEnum: + numnodes = 6; + 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; + break; case CrouzeixRaviartEnum: numnodes = 10; Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18179) @@ -406,4 +406,5 @@ switch(finiteelement){ + case NoneEnum: return 0; case P0Enum: return NUMNODESP0; case P1Enum: return NUMNODESP1; @@ -419,4 +420,5 @@ case MINIEnum: return NUMNODESP1b+NUMNODESP1; case TaylorHoodEnum: return NUMNODESP2+NUMNODESP1; + case LATaylorHoodEnum: return NUMNODESP2; case XTaylorHoodEnum: return NUMNODESP2+NUMNODESP1; case CrouzeixRaviartEnum: return NUMNODESP2b+NUMNODESP1; @@ -435,4 +437,5 @@ case MINIEnum: return P1bubbleEnum; case TaylorHoodEnum: return P2Enum; + case LATaylorHoodEnum: return P2Enum; case XTaylorHoodEnum: return P2Enum; case CrouzeixRaviartEnum:return P2bubbleEnum; @@ -451,4 +454,5 @@ case MINIEnum: return P1Enum; case TaylorHoodEnum: return P1Enum; + case LATaylorHoodEnum: return NoneEnum; case XTaylorHoodEnum: return P1Enum; case CrouzeixRaviartEnum: return P1DGEnum; Index: /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 18179) @@ -26,5 +26,5 @@ case StressbalanceSolutionEnum: - numanalyses=6; + numanalyses=7; analyses=xNew(numanalyses); analyses[0]=StressbalanceAnalysisEnum; @@ -34,4 +34,5 @@ analyses[4]=ExtrudeFromBaseAnalysisEnum; analyses[5]=DepthAverageAnalysisEnum; + analyses[6]=UzawaPressureAnalysisEnum; break; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 18179) @@ -378,4 +378,39 @@ } break; + case LATaylorHoodEnum: + _assert_(approximation==FSApproximationEnum); + /*P2 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"); + } + } + } + + /*No pressure*/ + break; case OneLayerP4zEnum: _assert_(approximation==FSApproximationEnum); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18178) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18179) @@ -390,4 +390,5 @@ ThermalSolutionEnum, TransientSolutionEnum, + UzawaPressureAnalysisEnum, GiaSolutionEnum, GiaAnalysisEnum, @@ -613,4 +614,5 @@ MINIcondensedEnum, TaylorHoodEnum, + LATaylorHoodEnum, XTaylorHoodEnum, OneLayerP4zEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18179) @@ -393,4 +393,5 @@ case ThermalSolutionEnum : return "ThermalSolution"; case TransientSolutionEnum : return "TransientSolution"; + case UzawaPressureAnalysisEnum : return "UzawaPressureAnalysis"; case GiaSolutionEnum : return "GiaSolution"; case GiaAnalysisEnum : return "GiaAnalysis"; @@ -602,4 +603,5 @@ case MINIcondensedEnum : return "MINIcondensed"; case TaylorHoodEnum : return "TaylorHood"; + case LATaylorHoodEnum : return "LATaylorHood"; case XTaylorHoodEnum : return "XTaylorHood"; case OneLayerP4zEnum : return "OneLayerP4z"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18178) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18179) @@ -402,4 +402,5 @@ else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; + else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum; else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; else if (strcmp(name,"Friction")==0) return FrictionEnum; - else if (strcmp(name,"Internal")==0) return InternalEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MassFlux")==0) return MassFluxEnum; + if (strcmp(name,"Internal")==0) return InternalEnum; + else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; else if (strcmp(name,"Misfit")==0) return MisfitEnum; @@ -614,4 +615,5 @@ else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; + else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"Time")==0) return TimeEnum; else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; - else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; - else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; + if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; + else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; + else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum; else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp (revision 18179) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la.cpp (revision 18179) @@ -0,0 +1,110 @@ +/*!\file: solutionsequence_la.cpp + * \brief: numerical core of la solutions + */ + +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +void solutionsequence_la(FemModel* femmodel){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug_old = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + Vector* pug = NULL; + Vector* pug_old = NULL; + IssmDouble eps_rel,r,theta; // 0 .15SetCurrentConfiguration(StressbalanceAnalysisEnum); + + /*Recover parameters: */ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum); + femmodel->parameters->FindParam(&r,AugmentedLagrangianREnum); + + /*Update constraints*/ + femmodel->UpdateConstraintsx(); + + /*Convergence criterion*/ + int count = 0; + GetSolutionFromInputsx(&ug,femmodel); + Vector* vx = NULL; + Vector* vx_old = NULL; + GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum); + + while(true){ + count++; + + /*save pointer to old velocity*/ + delete ug_old;ug_old=ug; + delete vx_old;vx_old=vx; + + /*Solve KU=F*/ + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + 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; + + /*Update solution*/ + InputUpdateFromSolutionx(femmodel,ug); + GetVectorFromInputsx(&vx,femmodel,VxEnum,VertexEnum); + + femmodel->SetCurrentConfiguration(UzawaPressureAnalysisEnum); + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + 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(&pug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; + + /*Update solution*/ + InputUpdateFromSolutionx(femmodel,pug); + + /*Check for convergence*/ + //Vector* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0); + //IssmDouble ndu=dug->Norm(NORM_TWO); delete dug; + //IssmDouble nu =ug_old->Norm(NORM_TWO); + Vector* dvx=vx_old->Duplicate(); vx_old->Copy(dvx); dvx->AYPX(vx,-1.0); + IssmDouble ndu=dvx->Norm(NORM_TWO); delete dvx; + IssmDouble nu =vx_old->Norm(NORM_TWO); + if (xIsNan(ndu) || xIsNan(nu)) _error_("convergence criterion is NaN!"); + if((ndu/nu) " << eps_rel*100 << " %\n"); + } + + if(count>=max_nonlinear_iterations){ + _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); + break; + } + } + + if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); + + delete ug; + delete ug_old; + delete pug; + delete pug_old; + delete vx; + delete vx_old; + delete stressanalysis; + delete pressureanalysis; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 18178) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 18179) @@ -18,4 +18,5 @@ void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads); void solutionsequence_linear(FemModel* femmodel); +void solutionsequence_la(FemModel* femmodel); void solutionsequence_la_theta(FemModel* femmodel); void solutionsequence_adjoint_linear(FemModel* femmodel); Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18178) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18179) @@ -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','CrouzeixRaviart'}); + md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LaTaylorHood','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 18178) +++ /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18179) @@ -47,5 +47,5 @@ string="%s\n%s"%(string,fielddisplay(self,'fe_SSA',"Finite Element for SSA: 'P1', 'P1bubble' 'P1bubblecondensed' 'P2'")) string="%s\n%s"%(string,fielddisplay(self,'fe_HO' ,"Finite Element for HO: 'P1' 'P1bubble' 'P1bubblecondensed' 'P1xP2' 'P2xP1' 'P2'")) - string="%s\n%s"%(string,fielddisplay(self,'fe_FS' ,"Finite Element for FS: 'P1P1' (debugging only) 'P1P1GLS' 'MINIcondensed' 'MINI' 'TaylorHood' 'XTaylorHood'")) + string="%s\n%s"%(string,fielddisplay(self,'fe_FS' ,"Finite Element for FS: 'P1P1' (debugging only) 'P1P1GLS' 'MINIcondensed' 'MINI' 'TaylorHood' 'LATaylorHood' 'XTaylorHood'")) string="%s\n%s"%(string,fielddisplay(self,'vertex_equation',"flow equation for each vertex")) string="%s\n%s"%(string,fielddisplay(self,'element_equation',"flow equation for each element")) Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18178) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18179) @@ -385,4 +385,5 @@ def ThermalSolutionEnum(): return StringToEnum("ThermalSolution")[0] def TransientSolutionEnum(): return StringToEnum("TransientSolution")[0] +def UzawaPressureAnalysisEnum(): return StringToEnum("UzawaPressureAnalysis")[0] def GiaSolutionEnum(): return StringToEnum("GiaSolution")[0] def GiaAnalysisEnum(): return StringToEnum("GiaAnalysis")[0] @@ -594,6 +595,8 @@ def MINIcondensedEnum(): return StringToEnum("MINIcondensed")[0] def TaylorHoodEnum(): return StringToEnum("TaylorHood")[0] +def LATaylorHoodEnum(): return StringToEnum("LATaylorHood")[0] def XTaylorHoodEnum(): return StringToEnum("XTaylorHood")[0] def OneLayerP4zEnum(): return StringToEnum("OneLayerP4z")[0] +def CrouzeixRaviartEnum(): return StringToEnum("CrouzeixRaviart")[0] def SaveResultsEnum(): return StringToEnum("SaveResults")[0] def BoolExternalResultEnum(): return StringToEnum("BoolExternalResult")[0] Index: /issm/trunk-jpl/src/m/enum/LATaylorHoodEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/LATaylorHoodEnum.m (revision 18179) +++ /issm/trunk-jpl/src/m/enum/LATaylorHoodEnum.m (revision 18179) @@ -0,0 +1,11 @@ +function macro=LATaylorHoodEnum() +%LATAYLORHOODENUM - Enum of LATaylorHood +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=LATaylorHoodEnum() + +macro=StringToEnum('LATaylorHood'); Index: /issm/trunk-jpl/src/m/enum/UzawaPressureAnalysisEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/UzawaPressureAnalysisEnum.m (revision 18179) +++ /issm/trunk-jpl/src/m/enum/UzawaPressureAnalysisEnum.m (revision 18179) @@ -0,0 +1,11 @@ +function macro=UzawaPressureAnalysisEnum() +%UZAWAPRESSUREANALYSISENUM - Enum of UzawaPressureAnalysis +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=UzawaPressureAnalysisEnum() + +macro=StringToEnum('UzawaPressureAnalysis');