Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18226) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18227) @@ -112,4 +112,5 @@ if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum){ parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum)); + parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum)); parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum)); } @@ -2879,6 +2880,7 @@ /*Intermediaries*/ int i,dim,epssize; - IssmDouble r,Jdet,viscosity,DU; + IssmDouble r,rl,Jdet,viscosity,DU,DUl; IssmDouble *xyz_list = NULL; + IssmDouble *xyz_list_base = NULL; /*Get problem dimension*/ @@ -2893,7 +2895,11 @@ if(dim==2) pnumnodes=3; else pnumnodes=6; + int lnumnodes; + if(dim==2) lnumnodes=2; + else lnumnodes=3; //int pnumnodes = element->NumberofNodes(P1Enum); int numdof = vnumnodes*dim; int pnumdof = pnumnodes; + int lnumdof = lnumnodes; /*Prepare coordinate system list*/ @@ -2910,4 +2916,7 @@ IssmDouble* BprimeU = xNew(numdof); IssmDouble* D = xNewZeroInit(epssize*epssize); + IssmDouble* CtCUzawa = xNewZeroInit(numdof*pnumdof); + IssmDouble* C = xNew(lnumdof); + IssmDouble* Cprime = xNew(numdof); /*Retrieve all inputs and parameters*/ @@ -2947,4 +2956,24 @@ } + if(element->IsOnBase()){ + element->FindParam(&rl,AugmentedLagrangianRlambdaEnum); + element->GetVerticesCoordinatesBase(&xyz_list_base); + + delete gauss; + Gauss* gauss=element->NewGaussBase(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + this->GetCFS(C,element,dim,xyz_list_base,gauss); + this->GetCFSprime(Cprime,element,dim,xyz_list_base,gauss); + DUl = gauss->weight*Jdet*sqrt(rl); + TripleMultiply(C,1,lnumdof,1, + &DU,1,1,0, + Cprime,1,numdof,0, + CtCUzawa,1); + } + } + /*Transform Coordinate System*/ element->TransformStiffnessMatrixCoord(Ke,cs_list); @@ -2953,4 +2982,9 @@ MatrixMultiply(BtBUzawa,pnumdof,numdof,1, BtBUzawa,pnumdof,numdof,0, + &Ke->values[0],1); + + /*The sigma naugmentation should not be transformed*/ + MatrixMultiply(CtCUzawa,lnumdof,numdof,1, + CtCUzawa,lnumdof,numdof,0, &Ke->values[0],1); @@ -2964,4 +2998,7 @@ xDelete(BU); xDelete(BtBUzawa); + xDelete(Cprime); + xDelete(C); + xDelete(CtCUzawa); xDelete(cs_list); return Ke; @@ -4372,4 +4409,60 @@ xDelete(vbasis); }/*}}}*/ +void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute C matrix. C=[Cp1 Cp2 ...] where Cpi=phi_pi. + */ + + /*Fetch number of nodes for this finite element*/ + int lnumnodes; + if(dim==2) lnumnodes=3; + else lnumnodes=6; + //int pnumnodes = element->NumberofNodes(P1Enum); + + /*Get nodal functions derivatives*/ + IssmDouble* basis =xNew(lnumnodes); + element->NodalFunctionsP1(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /* Compute C' matrix. C'=[C1' C2' C3'] + * Ci' = [ dphi/dx dphi/dy ] + * + * In 3d + * Ci=[ dh/dx dh/dy dh/dz ] + * where phi is the finiteelement function for node i. + */ + + /*Fetch number of nodes for this finite element*/ + int lnumnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(dim*lnumnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build C_prime: */ + if(dim==2){ + for(int i=0;i(dbasis); +}/*}}}*/ void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18226) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18227) @@ -86,4 +86,6 @@ 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 GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss); + void GetCFSprime(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/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18226) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18227) @@ -669,4 +669,6 @@ AugmentedLagrangianREnum, AugmentedLagrangianRhopEnum, + AugmentedLagrangianRlambdaEnum, + AugmentedLagrangianRholambdaEnum, AugmentedLagrangianThetaEnum, /*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18226) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18227) @@ -650,4 +650,6 @@ case AugmentedLagrangianREnum : return "AugmentedLagrangianR"; case AugmentedLagrangianRhopEnum : return "AugmentedLagrangianRhop"; + case AugmentedLagrangianRlambdaEnum : return "AugmentedLagrangianRlambda"; + case AugmentedLagrangianRholambdaEnum : return "AugmentedLagrangianRholambda"; case AugmentedLagrangianThetaEnum : return "AugmentedLagrangianTheta"; case NoneEnum : return "None"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18226) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18227) @@ -665,4 +665,6 @@ else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum; else if (strcmp(name,"AugmentedLagrangianRhop")==0) return AugmentedLagrangianRhopEnum; + else if (strcmp(name,"AugmentedLagrangianRlambda")==0) return AugmentedLagrangianRlambdaEnum; + else if (strcmp(name,"AugmentedLagrangianRholambda")==0) return AugmentedLagrangianRholambdaEnum; else if (strcmp(name,"AugmentedLagrangianTheta")==0) return AugmentedLagrangianThetaEnum; else if (strcmp(name,"None")==0) return NoneEnum; Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18226) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 18227) @@ -6,20 +6,22 @@ classdef flowequation properties (SetAccess=public) - isSIA = 0; - isSSA = 0; - isL1L2 = 0; - isHO = 0; - isFS = 0; - fe_SSA = ''; - fe_HO = ''; - fe_FS = ''; - augmented_lagrangian_r = 1.; - augmented_lagrangian_rhop = 1.; - XTH_theta = 0.; - vertex_equation = NaN; - element_equation = NaN; - borderSSA = NaN; - borderHO = NaN; - borderFS = NaN; + isSIA = 0; + isSSA = 0; + isL1L2 = 0; + isHO = 0; + isFS = 0; + fe_SSA = ''; + fe_HO = ''; + fe_FS = ''; + augmented_lagrangian_r = 1.; + augmented_lagrangian_rhop = 1.; + augmented_lagrangian_rlambda = 1.; + augmented_lagrangian_rholambda = 1.; + XTH_theta = 0.; + vertex_equation = NaN; + element_equation = NaN; + borderSSA = NaN; + borderHO = NaN; + borderFS = NaN; end methods (Static) @@ -137,5 +139,7 @@ md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'}); md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>',0.); + md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rlambda','numel',[1],'>',0.); md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rhop','numel',[1],'>',0.); + md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rholambda','numel',[1],'>',0.); md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5); md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]); @@ -196,4 +200,6 @@ WriteData(fid,'enum',AugmentedLagrangianREnum(),'data',obj.augmented_lagrangian_r ,'format','Double'); WriteData(fid,'enum',AugmentedLagrangianRhopEnum(),'data',obj.augmented_lagrangian_rhop ,'format','Double'); + WriteData(fid,'enum',AugmentedLagrangianRlambdaEnum(),'data',obj.augmented_lagrangian_rlambda ,'format','Double'); + WriteData(fid,'enum',AugmentedLagrangianRholambdaEnum(),'data',obj.augmented_lagrangian_rholambda ,'format','Double'); WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',obj.XTH_theta ,'format','Double'); WriteData(fid,'object',obj,'fieldname','borderSSA','format','DoubleMat','mattype',1); Index: /issm/trunk-jpl/src/m/classes/flowequation.py =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18226) +++ /issm/trunk-jpl/src/m/classes/flowequation.py (revision 18227) @@ -17,20 +17,22 @@ def __init__(self): # {{{ - self.isSIA = 0 - self.isSSA = 0 - self.isL1L2 = 0 - self.isHO = 0 - self.isFS = 0 - self.fe_SSA = '' - self.fe_HO = '' - self.fe_FS = '' - self.augmented_lagrangian_r = 1. - self.augmented_lagrangian_rhop = 1. - self.XTH_theta = 0. - self.vertex_equation = float('NaN') - self.element_equation = float('NaN') - self.borderSSA = float('NaN') - self.borderHO = float('NaN') - self.borderFS = float('NaN') + self.isSIA = 0 + self.isSSA = 0 + self.isL1L2 = 0 + self.isHO = 0 + self.isFS = 0 + self.fe_SSA = '' + self.fe_HO = '' + self.fe_FS = '' + self.augmented_lagrangian_r = 1. + self.augmented_lagrangian_rhop = 1. + self.augmented_lagrangian_rlambda = 1. + self.augmented_lagrangian_rholambda = 1. + self.XTH_theta = 0. + self.vertex_equation = float('NaN') + self.element_equation = float('NaN') + self.borderSSA = float('NaN') + self.borderHO = float('NaN') + self.borderFS = float('NaN') #set defaults @@ -85,4 +87,6 @@ md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>',0.) md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rhop','numel',[1],'>',0.) + md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rlambda','numel',[1],'>',0.) + md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_rholambda','numel',[1],'>',0.) md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',.5) if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'): @@ -115,4 +119,6 @@ WriteData(fid,'enum',AugmentedLagrangianREnum(),'data',self.augmented_lagrangian_r ,'format','Double') WriteData(fid,'enum',AugmentedLagrangianRhopEnum(),'data',self.augmented_lagrangian_rhop ,'format','Double') + WriteData(fid,'enum',AugmentedLagrangianRlambdaEnum(),'data',self.augmented_lagrangian_rlambda ,'format','Double') + WriteData(fid,'enum',AugmentedLagrangianRholambdaEnum(),'data',self.augmented_lagrangian_rholambda ,'format','Double') WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',self.XTH_theta ,'format','Double') WriteData(fid,'object',self,'fieldname','borderSSA','format','DoubleMat','mattype',1) Index: /issm/trunk-jpl/src/m/enum/AugmentedLagrangianRholambdaEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/AugmentedLagrangianRholambdaEnum.m (revision 18227) +++ /issm/trunk-jpl/src/m/enum/AugmentedLagrangianRholambdaEnum.m (revision 18227) @@ -0,0 +1,11 @@ +function macro=AugmentedLagrangianRholambdaEnum() +%AUGMENTEDLAGRANGIANRHOLAMBDAENUM - Enum of AugmentedLagrangianRholambda +% +% 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=AugmentedLagrangianRholambdaEnum() + +macro=StringToEnum('AugmentedLagrangianRholambda'); Index: /issm/trunk-jpl/src/m/enum/AugmentedLagrangianRlambdaEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/AugmentedLagrangianRlambdaEnum.m (revision 18227) +++ /issm/trunk-jpl/src/m/enum/AugmentedLagrangianRlambdaEnum.m (revision 18227) @@ -0,0 +1,11 @@ +function macro=AugmentedLagrangianRlambdaEnum() +%AUGMENTEDLAGRANGIANRLAMBDAENUM - Enum of AugmentedLagrangianRlambda +% +% 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=AugmentedLagrangianRlambdaEnum() + +macro=StringToEnum('AugmentedLagrangianRlambda'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18226) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18227) @@ -642,4 +642,6 @@ def AugmentedLagrangianREnum(): return StringToEnum("AugmentedLagrangianR")[0] def AugmentedLagrangianRhopEnum(): return StringToEnum("AugmentedLagrangianRhop")[0] +def AugmentedLagrangianRlambdaEnum(): return StringToEnum("AugmentedLagrangianRlambda")[0] +def AugmentedLagrangianRholambdaEnum(): return StringToEnum("AugmentedLagrangianRholambda")[0] def AugmentedLagrangianThetaEnum(): return StringToEnum("AugmentedLagrangianTheta")[0] def NoneEnum(): return StringToEnum("None")[0]