Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25609) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25610) @@ -31,5 +31,5 @@ IssmDouble rho_ice; IssmDouble FSreconditioning; - bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling; + bool isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling; bool spcpresent = false; int Mx,Nx; @@ -59,12 +59,13 @@ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); + iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); - /*Now, is the flag macayaealHO on? otherwise, do nothing: */ - if(!isSSA && !isHO && !isFS && !isL1L2) return; + /*Is this model only SIA??*/ + if(!isSSA && !isHO && !isFS && !isL1L2 && !isMLHO) return; /*Do we have coupling*/ - if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) iscoupling = true; else @@ -77,4 +78,5 @@ if(isSSA) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_SSA"); else if(isL1L2) finiteelement = P1Enum; + else if(isMLHO) finiteelement = P1Enum; else if(isHO) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO"); else if(isFS){ iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS"); @@ -452,5 +454,5 @@ int count; int penpair_ids[2]; - bool isSSA,isL1L2,isHO,isFS; + bool isSSA,isL1L2,isMLHO,isHO,isFS; int numpenalties,numrifts,numriftsegments; IssmDouble *riftinfo = NULL; @@ -460,4 +462,5 @@ /*Fetch parameters: */ iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); + iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); @@ -465,6 +468,6 @@ iomodel->FindConstant(&numrifts,"md.rifts.numrifts"); - /*Now, is the flag macayaealHO on? otherwise, do nothing: */ - if(!isSSA && !isHO && !isFS && !isL1L2) return; + /*Is this SIA only?*/ + if(!isSSA && !isHO && !isFS && !isL1L2 && !isMLHO) return; /*Initialize counter: */ @@ -511,5 +514,5 @@ /*Intermediary*/ - bool isSSA,isL1L2,isHO,isFS,iscoupling; + bool isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling; int finiteelement=-1,approximation=-1; @@ -517,12 +520,13 @@ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); + iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); /*Now, check that we have non SIA elements */ - if(!isSSA & !isL1L2 & !isHO & !isFS) return; + if(!isSSA && !isL1L2 && !isMLHO && !isHO && !isFS) return; /*Do we have coupling*/ - if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) iscoupling = true; else @@ -539,4 +543,8 @@ else if(isL1L2){ approximation = L1L2ApproximationEnum; + finiteelement = P1Enum; + } + else if(isMLHO){ + approximation = MLHOApproximationEnum; finiteelement = P1Enum; } @@ -613,5 +621,6 @@ } break; - case L1L2ApproximationEnum: numdofs =2; break; + case L1L2ApproximationEnum: numdofs = 2; break; + case MLHOApproximationEnum: numdofs = 4; break; case HOApproximationEnum: switch(domaintype){ @@ -678,5 +687,5 @@ int FrictionCoupling; int* finiteelement_list=NULL; - bool isSSA,isL1L2,isHO,isFS,iscoupling; + bool isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling; bool control_analysis; bool dakota_analysis; @@ -686,4 +695,5 @@ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); + iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); @@ -695,5 +705,5 @@ /*return if no processing required*/ - if(!isSSA & !isL1L2 & !isHO & !isFS) return; + if(!isSSA && !isL1L2 && !isMLHO && !isHO && !isFS) return; /*Fetch data needed and allocate vectors: */ @@ -702,5 +712,5 @@ /*Do we have coupling*/ - if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) iscoupling = true; else @@ -711,4 +721,5 @@ if(isSSA) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_SSA"); else if(isL1L2) finiteelement = P1Enum; + else if(isMLHO) finiteelement = P1Enum; else if(isHO) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_HO"); else if(isFS) iomodel->FindConstant(&finiteelement,"md.flowequation.fe_FS"); @@ -926,4 +937,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isSSA",FlowequationIsSSAEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isL1L2",FlowequationIsL1L2Enum)); + parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isMLHO",FlowequationIsMLHOEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isHO",FlowequationIsHOEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum)); @@ -1023,5 +1035,5 @@ /*Intermediaries*/ - bool isSSA,isL1L2,isHO,isFS; + bool isSSA,isL1L2,isMLHO,isHO,isFS; bool conserve_loads = true; int newton,domaintype,fe_FS; @@ -1030,4 +1042,5 @@ femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum); femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum); + femmodel->parameters->FindParam(&isMLHO,FlowequationIsMLHOEnum); femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); @@ -1036,5 +1049,5 @@ femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum); - if(isFS && !(isSSA || isHO || isL1L2)){ + if(isFS && !(isSSA || isHO || isL1L2 || isMLHO)){ if(VerboseSolution()) _printf0_(" computing velocities\n"); @@ -1061,5 +1074,5 @@ solutionsequence_nonlinear(femmodel,conserve_loads); } - else if(!isFS && (isSSA || isHO || isL1L2)){ + else if(!isFS && (isSSA || isHO || isL1L2 || isMLHO)){ if(VerboseSolution()) _printf0_(" computing velocities\n"); @@ -1077,5 +1090,5 @@ } } - else if ((isSSA || isL1L2 || isHO) && isFS){ + else if ((isSSA || isL1L2 || isMLHO || isHO) && isFS){ if(VerboseSolution()) _printf0_(" computing coupling between lower order models and FS\n"); solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads); @@ -1126,4 +1139,6 @@ case L1L2ApproximationEnum: return CreateKMatrixL1L2(element); + case MLHOApproximationEnum: + return CreateKMatrixMLHO(element); case HOApproximationEnum: return CreateKMatrixHO(element); @@ -1153,4 +1168,6 @@ case L1L2ApproximationEnum: return CreatePVectorL1L2(element); + case MLHOApproximationEnum: + return CreatePVectorMLHO(element); case HOApproximationEnum: return CreatePVectorHO(element); @@ -1177,8 +1194,5 @@ GetSolutionFromInputsFS(solution,element); return; - case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum: - GetSolutionFromInputsHoriz(solution,element); - return; - case L1L2ApproximationEnum: + case SSAApproximationEnum: case HOApproximationEnum: case L1L2ApproximationEnum: case MLHOApproximationEnum: case SIAApproximationEnum: GetSolutionFromInputsHoriz(solution,element); return; @@ -1261,4 +1275,7 @@ case L1L2ApproximationEnum: InputUpdateFromSolutionL1L2(solution,element); + return; + case MLHOApproximationEnum: + InputUpdateFromSolutionMLHO(solution,element); return; case SSAHOApproximationEnum: @@ -2388,11 +2405,8 @@ /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); - int numdof = numnodes*2; /*Initialize Element matrix and vectors*/ ElementMatrix* Ke = basalelement->NewElementMatrix(L1L2ApproximationEnum); - IssmDouble* B = xNew(3*numdof); - IssmDouble* Bprime = xNew(3*numdof); - IssmDouble* D = xNewZeroInit(3*3); + IssmDouble* dbasis = xNew(3*numnodes); /*Retrieve all inputs and parameters*/ @@ -2409,15 +2423,24 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); - this->GetBSSA(B,basalelement,2,xyz_list,gauss_base); - this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); - for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet; - - TripleMultiply(B,3,numdof,1, - D,3,3,0, - Bprime,3,numdof,0, - &Ke->values[0],1); + for(int i=0;ivalues[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*( + 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + ); + Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*( + 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + ); + Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*( + 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + ); + Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*( + dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + ); + } + } } @@ -2430,7 +2453,5 @@ basalelement->DeleteMaterials(); delete basalelement; xDelete(xyz_list); - xDelete(D); - xDelete(Bprime); - xDelete(B); + xDelete(dbasis); return Ke; }/*}}}*/ @@ -2573,4 +2594,328 @@ }/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/ + + int i,dim,domaintype; + IssmDouble rho_ice,g; + int* doflist=NULL; + IssmDouble* xyz_list=NULL; + Element* basalelement=NULL; + + /*Deal with pressure first*/ + int numvertices = element->GetNumberOfVertices(); + IssmDouble* pressure = xNew(numvertices); + IssmDouble* thickness = xNew(numvertices); + IssmDouble* surface = xNew(numvertices); + + element->FindParam(&dim,DomainDimensionEnum); + element->FindParam(&domaintype,DomainTypeEnum); + rho_ice =element->FindParam(MaterialsRhoIceEnum); + g =element->FindParam(ConstantsGEnum); + if(dim==2){ + element->GetInputListOnVertices(thickness,ThicknessEnum); + for(i=0;iGetVerticesCoordinates(&xyz_list); + element->GetInputListOnVertices(surface,SurfaceEnum); + for(i=0;iAddInput(PressureEnum,pressure,P1Enum); + xDelete(pressure); + xDelete(thickness); + xDelete(surface); + + /*Get basal element*/ + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()){xDelete(xyz_list); return;} + basalelement=element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + int numdof = numnodes*2; + + /*Fetch dof list and allocate solution vectors*/ + basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numdof); + IssmDouble* vx = xNew(numnodes); + IssmDouble* vy = xNew(numnodes); + IssmDouble* vz = xNew(numnodes); + IssmDouble* vel = xNew(numnodes); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iTransformSolutionCoord(&values[0],XYEnum); + basalelement->FindParam(&domaintype,DomainTypeEnum); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + for(i=0;i(vx[i])) _error_("NaN found in solution vector"); + if(xIsInf(vx[i])) _error_("Inf found in solution vector"); + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + if(xIsInf(vy[i])) _error_("Inf found in solution vector"); + } + + /*Get Vz and compute vel*/ + basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); + for(i=0;iAddBasalInput(VxEnum,vx,element->GetElementType()); + element->AddBasalInput(VyEnum,vy,element->GetElementType()); + element->AddBasalInput(VelEnum,vel,element->GetElementType()); + + /*Free ressources:*/ + xDelete(vel); + xDelete(vz); + xDelete(vy); + xDelete(vx); + xDelete(values); + xDelete(xyz_list); + xDelete(doflist); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; +}/*}}}*/ + +/*MLHO*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ + _error_("not implemented yet"); + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixMLHOViscous(element); + ElementMatrix* Ke2=CreateKMatrixMLHOFriction(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOFriction(Element* element){/*{{{*/ + _error_("not implemented yet"); + + if(!element->IsOnBase() || element->IsFloating()) return NULL; + Element* basalelement = element->SpawnBasalElement(); + ElementMatrix* Ke = CreateKMatrixSSAFriction(basalelement); + + /*clean-up and return*/ + basalelement->DeleteMaterials(); delete basalelement; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOViscous(Element* element){/*{{{*/ + _error_("not implemented yet"); + + /*Intermediaries*/ + IssmDouble viscosity,Jdet; + IssmDouble *xyz_list = NULL; + + /*Get element on base*/ + Element* basalelement = element->GetBasalElement()->SpawnBasalElement(true); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + int numdof = numnodes*2; + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum); + IssmDouble* B = xNew(3*numdof); + IssmDouble* Bprime = xNew(3*numdof); + IssmDouble* D = xNewZeroInit(3*3); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + Gauss* gauss_base = basalelement->NewGauss(); + while(gauss->next()){ + gauss->SynchronizeGaussBase(gauss_base); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + this->GetBSSA(B,basalelement,2,xyz_list,gauss_base); + this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base); + + element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); + + for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet; + + TripleMultiply(B,3,numdof,1, + D,3,3,0, + Bprime,3,numdof,0, + &Ke->values[0],1); + } + + /*Transform Coordinate System*/ + basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum); + + /*Clean up and return*/ + delete gauss; + delete gauss_base; + basalelement->DeleteMaterials(); delete basalelement; + xDelete(xyz_list); + xDelete(D); + xDelete(Bprime); + xDelete(B); + return Ke; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorMLHO(Element* element){/*{{{*/ + _error_("not implemented yet"); + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Intermediaries*/ + int domaintype; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: case Domain2DverticalEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<DeleteMaterials(); delete basalelement;}; + delete pe1; + delete pe2; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorMLHODrivingStress(Element* element){/*{{{*/ + _error_("not implemented yet"); + + /*Intermediaries */ + IssmDouble thickness,Jdet,slope[2]; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and vectors*/ + ElementVector* pe = element->NewElementVector(MLHOApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input); + IssmDouble rhog = element->FindParam(MaterialsRhoIceEnum)*element->FindParam(ConstantsGEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + while(gauss->next()){ + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis, gauss); + + thickness_input->GetInputValue(&thickness,gauss); + surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); + + for(int i=0;ivalues[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,XYEnum); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete gauss; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorMLHOFront(Element* element){/*{{{*/ + _error_("not implemented yet"); + + /*If no front, return NULL*/ + if(!element->IsIcefront()) return NULL; + + /*Intermediaries*/ + IssmDouble Jdet,thickness,bed,sealevel,water_pressure,ice_pressure; + IssmDouble surface_under_water,base_under_water,pressure; + IssmDouble *xyz_list = NULL; + IssmDouble *xyz_list_front = NULL; + IssmDouble normal[2]; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(MLHOApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); + IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + IssmDouble gravity = element->FindParam(ConstantsGEnum); + element->GetVerticesCoordinates(&xyz_list); + element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + element->NormalSection(&normal[0],xyz_list_front); + + /*Start looping on Gaussian points*/ + Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3); + while(gauss->next()){ + thickness_input->GetInputValue(&thickness,gauss); + base_input->GetInputValue(&bed,gauss); + sealevel_input->GetInputValue(&sealevel,gauss); + element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); + element->NodalFunctions(basis,gauss); + + surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level + base_under_water = min(0.,bed-sealevel); // 0 if the bottom of the glacier is above water level + water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water); + ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness; + pressure = ice_pressure + water_pressure; + + for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; + pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,XYEnum); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(xyz_list_front); + xDelete(basis); + delete gauss; + return pe; +}/*}}}*/ +void StressbalanceAnalysis::InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element){/*{{{*/ + + _error_("not implemented yet"); int i,dim,domaintype; Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 25609) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 25610) @@ -56,4 +56,12 @@ ElementVector* CreatePVectorL1L2DrivingStress(Element* element); void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element); + /*MLHO*/ + ElementMatrix* CreateKMatrixMLHO(Element* element); + ElementMatrix* CreateKMatrixMLHOFriction(Element* element); + ElementMatrix* CreateKMatrixMLHOViscous(Element* element); + ElementVector* CreatePVectorMLHO(Element* element); + ElementVector* CreatePVectorMLHOFront(Element* element); + ElementVector* CreatePVectorMLHODrivingStress(Element* element); + void InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element); /*HO*/ ElementMatrix* CreateJacobianMatrixHO(Element* element); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 25609) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 25610) @@ -10,5 +10,5 @@ /*Intermediaries*/ - bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling; + bool isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling; /*Fetch parameters: */ @@ -16,4 +16,5 @@ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); + iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); @@ -23,5 +24,5 @@ /*Do we have coupling*/ - if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) iscoupling = true; else Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 25609) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 25610) @@ -10,5 +10,5 @@ /*Intermediary*/ - bool isSIA,isSSA,isL1L2,isHO,isFS,iscoupling; + bool isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,iscoupling; int Mz,Nz; IssmDouble *spcvz = NULL; @@ -21,9 +21,10 @@ iomodel->FindConstant(&isSSA,"md.flowequation.isSSA"); iomodel->FindConstant(&isL1L2,"md.flowequation.isL1L2"); + iomodel->FindConstant(&isMLHO,"md.flowequation.isMLHO"); iomodel->FindConstant(&isHO,"md.flowequation.isHO"); iomodel->FindConstant(&isFS,"md.flowequation.isFS"); /*Do we have coupling*/ - if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isMLHO?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) iscoupling = true; else Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25609) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25610) @@ -3192,4 +3192,5 @@ this->parameters->FindParam(&temp,FlowequationIsSSAEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isSSA")); this->parameters->FindParam(&temp,FlowequationIsL1L2Enum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isL1L2")); + this->parameters->FindParam(&temp,FlowequationIsMLHOEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isMLHO")); this->parameters->FindParam(&temp,FlowequationIsHOEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isHO")); this->parameters->FindParam(&temp,FlowequationIsFSEnum); iomodel->AddConstant(new IoConstant(temp,"md.flowequation.isFS")); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 25609) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 25610) @@ -117,4 +117,7 @@ } if(in_approximation==L1L2ApproximationEnum && !reCast(iomodel->Data("md.mesh.vertexonbase")[io_index])){ + this->HardDeactivate(); + } + if(in_approximation==MLHOApproximationEnum && !reCast(iomodel->Data("md.mesh.vertexonbase")[io_index])){ this->HardDeactivate(); }