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

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

@@ -37 +37 @@
                 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);

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -838 +838 @@
                         return NULL;
                 case SSAApproximationEnum: 
                         return CreatePVectorSSA(element);
+                case L1L2ApproximationEnum: 
+                        return CreatePVectorL1L2(element);
                 case HOApproximationEnum: 
                         return CreatePVectorHO(element);
                 case FSApproximationEnum: 
@@ -1338 +1340 @@
 }/*}}}*/
 
 /*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;