Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16891) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16892) @@ -820,4 +820,6 @@ case SSAApproximationEnum: return CreateKMatrixSSA(element); + case L1L2ApproximationEnum: + return CreateKMatrixL1L2(element); case HOApproximationEnum: return CreateKMatrixHO(element); @@ -1441,4 +1443,86 @@ /*L1L2*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2(Element* element){/*{{{*/ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixL1L2Viscous(element); + ElementMatrix* Ke2=CreateKMatrixL1L2Friction(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2Friction(Element* element){/*{{{*/ + + if(!element->IsOnBed() || !element->IsFloating()) return NULL; + Element* basalelement = element->SpawnBasalElement(); + ElementMatrix* Ke = CreateKMatrixSSAFriction(basalelement); + + /*clean-up and return*/ + basalelement->DeleteMaterials(); delete basalelement; + return Ke; +}/*}}}*/ +ElementMatrix* StressbalanceAnalysis::CreateKMatrixL1L2Viscous(Element* element){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity,Jdet; + IssmDouble *xyz_list = NULL; + + /*Get element on base*/ + Element* basalelement = element->GetBasalElement()->SpawnBasalElement(); + + /*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); + + /*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(); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussBase(gauss_base); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + this->GetBSSA(B,basalelement,xyz_list,gauss_base); + this->GetBSSAprime(Bprime,basalelement,xyz_list,gauss_base); + + element->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::CreatePVectorL1L2(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16891) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16892) @@ -39,4 +39,7 @@ void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element); /*L1L2*/ + ElementMatrix* CreateKMatrixL1L2(Element* element); + ElementMatrix* CreateKMatrixL1L2Viscous(Element* element); + ElementMatrix* CreateKMatrixL1L2Friction(Element* element); ElementVector* CreatePVectorL1L2(Element* element); ElementVector* CreatePVectorL1L2DrivingStress(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16892) @@ -170,4 +170,5 @@ virtual void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; virtual void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; + virtual void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0; virtual void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16891) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16892) @@ -3540,4 +3540,72 @@ /*Assign output pointer*/ *pviscosity=viscosity; +} +/*}}}*/ +/*FUNCTION Penta::ViscosityL1L2{{{*/ +void Penta::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){ + /*Compute the L1L2 viscosity + * + * 1 + * mu = - A^-1 (sigma'_e)^(1-n) + * 2 + * + * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18) + * + * L1L2 assumptions: + * + * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//| + * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2 + * + * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|) + * + * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0 */ + + int i; + IssmDouble z,s,viscosity,p,q,delta; + IssmDouble tau_perp,tau_par,eps_b,A; + IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/ + IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/ + IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ + IssmDouble z_list[NUMVERTICES]; + IssmDouble slope[3]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); + + /*Get tau_perp*/ + for(i=0;iGetInputValue(&s,gauss); + surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); + PentaRef::GetInputValue(&z,&z_list[0],gauss); + tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + + /* Get eps_b*/ + vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; + eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]); + if(eps_b==0.){ + *pviscosity = 2.5e+17; + return; + } + + /*Get A*/ + _assert_(material->GetN()==3.0); + A=material->GetA(); + + /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/ + p = tau_perp *tau_perp; + q = - eps_b/A; + delta = q *q + p*p*p*4./27.; + _assert_(delta>0); + tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.); + + /*Viscosity*/ + viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp)); + _assert_(!xIsNan(viscosity)); + _assert_(viscosity > 0.); + + /*Assign output pointer*/ + *pviscosity = viscosity; } /*}}}*/ @@ -6849,5 +6917,5 @@ gauss->GaussPoint(ig); - gauss->SynchronizeGaussTria(gauss_tria); + gauss->SynchronizeGaussBase(gauss_tria); GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); @@ -7054,5 +7122,5 @@ gauss->GaussPoint(ig); - gauss->SynchronizeGaussTria(gauss_tria); + gauss->SynchronizeGaussBase(gauss_tria); GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); @@ -7470,5 +7538,5 @@ gauss->GaussPoint(ig); - gauss->SynchronizeGaussTria(gauss_tria); + gauss->SynchronizeGaussBase(gauss_tria); GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); @@ -7618,5 +7686,5 @@ gauss->GaussPoint(ig); - gauss->SynchronizeGaussTria(gauss_tria); + gauss->SynchronizeGaussBase(gauss_tria); GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); @@ -7625,5 +7693,5 @@ /*Get viscosity for L1L2 model*/ - GetL1L2Viscosity(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input); + ViscosityL1L2(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input); for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet; @@ -9576,73 +9644,4 @@ delete gauss; return Ke; -} -/*}}}*/ -/*FUNCTION Penta::GetL1L2Viscosity{{{*/ -void Penta::GetL1L2Viscosity(IssmDouble* pviscosity,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input,Input* surface_input){ - /*Compute the L1L2 viscosity - * - * 1 - * mu = - A^-1 (sigma'_e)^(1-n) - * 2 - * - * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18) - * - * L1L2 assumptions: - * - * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//| - * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2 - * - * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|) - * - * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0 */ - - int i; - IssmDouble z,s,viscosity,p,q,delta; - IssmDouble tau_perp,tau_par,eps_b,A; - IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/ - IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/ - IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ - IssmDouble z_list[NUMVERTICES]; - IssmDouble slope[3]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); - - /*Get tau_perp*/ - for(i=0;iGetInputValue(&s,gauss); - surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); - PentaRef::GetInputValue(&z,&z_list[0],gauss); - tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); - - /* Get eps_b*/ - vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss); - vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss); - for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]); - if(eps_b==0.){ - *pviscosity = 2.5e+17; - return; - } - - /*Get A*/ - _assert_(material->GetN()==3.0); - A=material->GetA(); - - /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/ - p = tau_perp *tau_perp; - q = - eps_b/A; - delta = q *q + p*p*p*4./27.; - _assert_(delta>0); - tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.); - - /*Viscosity*/ - viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp)); - _assert_(!xIsNan(viscosity)); - _assert_(viscosity > 0.); - - /*Assign output pointer*/ - *pviscosity = viscosity; - return; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16892) @@ -297,4 +297,5 @@ void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); @@ -366,5 +367,4 @@ ElementVector* CreatePVectorStressbalanceVertVolume(void); ElementVector* CreatePVectorStressbalanceVertBase(void); - void GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*); #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16892) @@ -172,4 +172,5 @@ void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");}; void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; + void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not implemented yet");}; void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");}; void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16892) @@ -329,4 +329,5 @@ void ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");}; + void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not implemented yet");}; void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/gauss/Gauss.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/Gauss.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/Gauss.h (revision 16892) @@ -19,4 +19,5 @@ virtual void GaussVertex(int iv)=0; virtual void GaussNode(int finitelement,int iv)=0; + virtual void SynchronizeGaussBase(Gauss* gauss)=0; }; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 16892) @@ -674,6 +674,9 @@ } /*}}}*/ -/*FUNCTION GaussPenta::SynchronizeGaussTria{{{*/ -void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){ +/*FUNCTION GaussPenta::SynchronizeGaussBase{{{*/ +void GaussPenta::SynchronizeGaussBase(Gauss* gauss){ + + _assert_(gauss->Enum()==GaussTriaEnum); + GaussTria* gauss_tria = dynamic_cast(gauss); gauss_tria->coord1=this->coord1; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h (revision 16892) @@ -48,5 +48,5 @@ void GaussNode(int finitelement,int iv); void GaussFaceTria(int index1, int index2, int index3, int order); - void SynchronizeGaussTria(GaussTria* gauss_tria); + void SynchronizeGaussBase(Gauss* gauss); }; #endif Index: /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 16892) @@ -123,2 +123,8 @@ } /*}}}*/ +/*FUNCTION GaussSeg::SynchronizeGaussBase{{{*/ +void GaussSeg::SynchronizeGaussBase(Gauss* gauss){ + + _error_("not supported"); +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h (revision 16892) @@ -34,4 +34,5 @@ void GaussVertex(int iv); void GaussNode(int finitelement,int iv); + void SynchronizeGaussBase(Gauss* gauss); }; #endif Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 16892) @@ -503,2 +503,8 @@ } /*}}}*/ +/*FUNCTION GaussTria::SynchronizeGaussBase{{{*/ +void GaussTria::SynchronizeGaussBase(Gauss* gauss){ + + _error_("not supported"); +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 16891) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 16892) @@ -44,4 +44,5 @@ void GaussNode(int finitelement,int iv); void GaussEdgeCenter(int index1,int index2); + void SynchronizeGaussBase(Gauss* gauss); }; #endif /* _GAUSSTRIA_H_ */