Changeset 17411

Timestamp:

03/11/14 16:45:19 (11 years ago)

Author:

seroussi

Message:

NEW: added ice shelf dampening for FS

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (3 diffs)
• StressbalanceAnalysis.h (modified) (1 diff)

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

@@ -7 +7 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 4
+//#define FSANALYTICAL 10
 
 /*Model processing*/
@@ -2810 +2810 @@
         ElementMatrix* Ke1=CreateKMatrixFSViscous(element);
         ElementMatrix* Ke2=CreateKMatrixFSFriction(element);
-        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+        ElementMatrix* Ke3=CreateKMatrixFSShelf(element);
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
         /*clean-up and return*/
         delete Ke1;
         delete Ke2;
+        delete Ke3;
         return Ke;
 }/*}}}*/
@@ -2979 +2981 @@
         xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(D);
+        return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
+
+        if(!element->IsFloating() || !element->IsOnBed()) return NULL;
+
+        /*If on not water or not FS, skip stiffness: */
+        int approximation,shelf_dampening;
+        element->GetInputValue(&approximation,ApproximationEnum);
+        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+        element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
+        if(shelf_dampening==0) return NULL;
+
+        /*Intermediaries*/
+        bool        mainlyfloating;
+        int         j,i,meshtype,dim;
+        IssmDouble  Jdet,slope2,scalar,dt;
+        IssmDouble  slope[2];
+        IssmDouble *xyz_list_base = NULL;
+        Gauss*      gauss         = NULL;
+
+        /*Get problem dimension*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2;break;
+                case Mesh3DEnum:         dim = 3;break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = element->GetNumberOfNodesVelocity();
+        int pnumnodes = element->GetNumberOfNodesPressure();
+        int numdof    = vnumnodes*dim + pnumnodes;
+
+        /*Initialize Element matrix and vectors*/
+        ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
+        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        if(dt==0) dt=pow(10,5);
+        IssmDouble  rho_water =element->GetMaterialParameter(MaterialsRhoWaterEnum);
+        IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
+        Input*      surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=element->NewGaussBase(3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                surface_input->GetInputDerivativeValue(&slope[0],xyz_list_base,gauss);
+                element->NodalFunctionsVelocity(vbasis,gauss);
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                if(dim==2) slope2=slope[0]*slope[0];
+                else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
+                scalar  = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; 
+                for(i=0;i<vnumnodes;i++){
+                        for(j=0;j<vnumnodes;j++){
+                                Ke->values[numdof*(i*dim+1)+j*dim+1] += scalar*vbasis[i]*vbasis[j];
+                        }
+                }
+        }
+
+        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+        /*Clean up and return*/
+        delete gauss;
+        xDelete<IssmDouble>(xyz_list_base);
+        xDelete<IssmDouble>(vbasis);
         return Ke;
 }/*}}}*/

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

@@ -69 +69 @@
                 ElementMatrix* CreateKMatrixFSViscous(Element* element);
                 ElementMatrix* CreateKMatrixFSFriction(Element* element);
+                ElementMatrix* CreateKMatrixFSShelf(Element* element);
                 ElementVector* CreatePVectorFS(Element* element);
                 ElementVector* CreatePVectorFSViscous(Element* element);