Changeset 19154

Timestamp:

02/25/15 09:27:11 (10 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added simplistic reset levelset function and calving law

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Element.cpp (modified) (1 diff)
• Element.h (modified) (1 diff)
• Tria.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -1367 +1367 @@
 IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
         return this->matpar->PureIceEnthalpy(pressure);
+}/*}}}*/
+void       Element::ResetIceLevelset(void){/*{{{*/
+
+        /*Get Ice Levelset on vertices*/
+        int         numvertices = this->GetNumberOfVertices();
+        IssmDouble *values      = xNew<IssmDouble>(numvertices);
+        this->GetInputListOnVertices(values,MaskIceLevelsetEnum);
+
+        /*Project on {-1,+1}*/
+        for(int i=0;i<numvertices;i++){
+                if(values[i]>0) values[i]=+1.;
+                else            values[i]=-1.;
+        }
+        this->AddInput(MaskIceLevelsetEnum,values,P1Enum);
+
+        /*Cleanup and return*/
+        xDelete<IssmDouble>(values);
+
 }/*}}}*/
 void       Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -123 +123 @@
                 ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
                 IssmDouble         PureIceEnthalpy(IssmDouble pressure);
+                void               ResetIceLevelset(void);
                 void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum);
                 void               ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -336 +336 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
+        IssmDouble  sigma_vm,sigma_max;
 
         /* Get node coordinates and dof list: */
@@ -341 +342 @@
 
         /*Retrieve all inputs and parameters we will need*/
-        Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
+        Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
+        Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
+        IssmDouble  B   = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
+        IssmDouble  n   = this->GetMaterialParameter(MaterialsRheologyNEnum);
 
         /* Start looping on the number of vertices: */
@@ -366 +369 @@
                 lambda2 = max(lambda2,0.);
 
+                /*Calculate sigma_vm*/
+                sigma_vm  = sqrt(3./2.) * B * pow(lambda1*lambda1 + lambda2*lambda2,1./(2.*n));
+                sigma_max = 350.e+3;
+
                 /*Assign values*/
-                calvingratex[iv]=vx*pow(lambda1 + lambda2,1./3.)*3.e+2;
-                calvingratey[iv]=vy*pow(lambda1 + lambda2,1./3.)*3.e+2;
+                calvingratex[iv]=vx*sigma_vm/sigma_max;
+                calvingratey[iv]=vy*sigma_vm/sigma_max;
                 calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
         }