Changeset 16804

Timestamp:

11/16/13 13:48:18 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: some clean up

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (5 diffs)
• StressbalanceAnalysis.h (modified) (1 diff)
• StressbalanceVerticalAnalysis.cpp (modified) (5 diffs)

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

@@ -812 +812 @@
                 case SSAApproximationEnum: 
                         return CreatePVectorSSA(element);
+                case HOApproximationEnum: 
+                        return CreatePVectorHO(element);
                 default:
                         _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
         }
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
+
+        /*compute all load vectors for this element*/
+        ElementVector* pe1=CreatePVectorHODrivingStress(element);
+        ElementVector* pe2=CreatePVectorHOFront(element);
+        ElementVector* pe =new ElementVector(pe1,pe2);
+
+        /*clean-up and return*/
+        delete pe1;
+        delete pe2;
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble  Jdet,slope[3];
+        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(HOApproximationEnum);
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        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(3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis, gauss);
+                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
+                        pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
+                }
+        }
+
+        /*Transform coordinate system*/
+        element->TransformLoadVectorCoord(pe,XYEnum);
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(xyz_list);
+        delete gauss;
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
+
+        return NULL;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
@@ -849 +909 @@
 
         /*Intermediaries */
-        int         i;
-        IssmDouble  driving_stress_baseline,thickness;
-        IssmDouble  Jdet;
-        IssmDouble  slope[2];
-        IssmDouble* xyz_list     = NULL;
+        IssmDouble  thickness,Jdet,slope[2];
+        IssmDouble* xyz_list = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = element->GetNumberOfNodes();
-        IssmDouble* icefrontlevel= xNew<IssmDouble>(numnodes);
 
         /*Initialize Element vector and vectors*/
@@ -865 +921 @@
         /*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);
-        Input* drag_input     =element->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
-        element->GetInputListOnVertices(icefrontlevel,BedEnum);
+        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: */
@@ -881 +936 @@
                 thickness_input->GetInputValue(&thickness,gauss);
                 surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-                driving_stress_baseline=element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum)*thickness;
-
-                /*Build load vector: */
-                for(i=0;i<numnodes;i++){
-                        pe->values[i*2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
-                        pe->values[i*2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
+
+                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];
                 }
         }
@@ -894 +947 @@
 
         /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
         delete gauss;

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

@@ -23 +23 @@
                 ElementMatrix* CreateKMatrix(Element* element);
                 ElementVector* CreatePVector(Element* element);
+                ElementVector* CreatePVectorHO(Element* element);
+                ElementVector* CreatePVectorHODrivingStress(Element* element);
+                ElementVector* CreatePVectorHOFront(Element* element);
                 ElementVector* CreatePVectorSSA(Element* element);
                 ElementVector* CreatePVectorSSADrivingStress(Element* element);

issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

@@ -115 +115 @@
         /*Intermediaries*/
         int         approximation;
-        IssmDouble  Jdet;
-        IssmDouble  dudx,dvdy,dwdz;
+        IssmDouble  Jdet,dudx,dvdy,dwdz;
         IssmDouble  du[3],dv[3],dw[3];
         IssmDouble* xyz_list = NULL;
@@ -167 +166 @@
 
         /*Intermediaries */
-        int         i,j,approximation;
+        int         approximation;
         IssmDouble *xyz_list      = NULL;
         IssmDouble *xyz_list_base = NULL;
-        IssmDouble  Jdet;
-        IssmDouble  vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
-        IssmDouble  slope[3];
+        IssmDouble  Jdet,slope[3];
+        IssmDouble  vx,vy,vz=0.,dbdx,dbdy,basalmeltingvalue;
 
         if(!element->IsOnBed()) return NULL;
@@ -208 +206 @@
                         vzFS_input->GetInputValue(&vz,gauss);
                 }
-                else vz=0.;
-
                 dbdx=slope[0];
                 dbdy=slope[1];
@@ -216 +212 @@
                 element->NodalFunctions(basis,gauss);
 
-                for(i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
+                for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
         }
 
@@ -225 +221 @@
         xDelete<IssmDouble>(xyz_list_base);
         return pe;
-
 }/*}}}*/
 void StressbalanceVerticalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/