source: issm/oecreview/Archive/16554-17801/ISSM-16858-16859.diff@ 17802

Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16858)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16859)
@@ -37,6 +37,9 @@
                void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
                /*L1L2*/
+               ElementVector* CreatePVectorL1L2(Element* element);
+               ElementVector* CreatePVectorL1L2DrivingStress(Element* element);
+               ElementVector* CreatePVectorL1L2Front(Element* element);
                void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
                /*HO*/
                ElementMatrix* CreateKMatrixHO(Element* element);
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16858)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16859)
@@ -838,6 +838,8 @@
                        return NULL;
                case SSAApproximationEnum: 
                        return CreatePVectorSSA(element);
+               case L1L2ApproximationEnum: 
+                       return CreatePVectorL1L2(element);
                case HOApproximationEnum: 
                        return CreatePVectorHO(element);
                case FSApproximationEnum: 
@@ -1338,6 +1340,144 @@
 }/*}}}*/
 
 /*L1L2*/
+ElementVector* StressbalanceAnalysis::CreatePVectorL1L2(Element* element){/*{{{*/
+
+       /*Intermediaries*/
+       int      meshtype;
+       Element* basalelement;
+
+       /*Get basal element*/
+       element->FindParam(&meshtype,MeshTypeEnum);
+       switch(meshtype){
+               case Mesh2DhorizontalEnum:
+                       basalelement = element;
+                       break;
+               case Mesh3DEnum:
+                       if(!element->IsOnBed()) return NULL;
+                       basalelement = element->SpawnBasalElement();
+                       break;
+               default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+       }
+
+       /*compute all load vectors for this element*/
+       ElementVector* pe1=CreatePVectorL1L2DrivingStress(basalelement);
+       ElementVector* pe2=CreatePVectorL1L2Front(basalelement);
+       ElementVector* pe =new ElementVector(pe1,pe2);
+
+       /*clean-up and return*/
+       if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+       delete pe1;
+       delete pe2;
+       return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorL1L2DrivingStress(Element* element){/*{{{*/
+
+       /*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(L1L2ApproximationEnum);
+       IssmDouble*    basis = xNew<IssmDouble>(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->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+
+       /* Start  looping on the number of gaussian points: */
+       Gauss* gauss=element->NewGauss(2);
+       for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+               gauss->GaussPoint(ig);
+
+               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;i<numnodes;i++){
+                       pe->values[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<IssmDouble>(xyz_list);
+       xDelete<IssmDouble>(basis);
+       delete gauss;
+       return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorL1L2Front(Element* element){/*{{{*/
+
+       /*If no front, return NULL*/
+       if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+
+       /*Intermediaries*/
+       IssmDouble  Jdet,thickness,bed,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(L1L2ApproximationEnum);
+       IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+       /*Retrieve all inputs and parameters*/
+       Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+       Input* bed_input       = element->GetInput(BedEnum);       _assert_(bed_input);
+       IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+       IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
+       IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
+       element->GetVerticesCoordinates(&xyz_list);
+       element->ZeroLevelsetCoordinates(&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);
+       for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+               gauss->GaussPoint(ig);
+               thickness_input->GetInputValue(&thickness,gauss);
+               bed_input->GetInputValue(&bed,gauss);
+               element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
+               element->NodalFunctions(basis,gauss);
+
+               surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
+               base_under_water    = min(0.,bed);           // 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;i<numnodes;i++){
+                       pe->values[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<IssmDouble>(xyz_list);
+       xDelete<IssmDouble>(xyz_list_front);
+       xDelete<IssmDouble>(basis);
+       delete gauss;
+       return pe;
+       return NULL;
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
 
        int         i,meshtype;