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

Timestamp:

01/22/15 14:51:25 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: testing new calving law

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (1 diff)

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

@@ -326 +326 @@
         delete gauss;
 
+}
+/*}}}*/
+void       Tria::CalvingRateDev(){/*{{{*/
+
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+        IssmDouble  calvingratex[NUMVERTICES];
+        IssmDouble  calvingratey[NUMVERTICES];
+        IssmDouble  calvingrate[NUMVERTICES];
+        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+        /*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);
+
+        /* Start looping on the number of vertices: */
+        GaussTria* gauss=new GaussTria();
+        for(int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Get velocity components and thickness*/
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                vel=sqrt(vx*vx+vy*vy)+1.e-14;
+
+                /*Compute strain rate and viscosity: */
+                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+
+                /*Get Eigen values*/
+                Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
+
+                /*Process Eigen values*/
+                lambda1>0? lambda1 = pow(lambda1,.3) : lambda1=0.;
+                lambda2>0? lambda2 = pow(lambda2,.3) : lambda2=0.;
+                lambda1 = lambda1*5.e-2;
+                lambda2 = lambda2*5.e-2;
+                _assert_(!xIsNan<IssmDouble>(lambda1));
+                _assert_(!xIsNan<IssmDouble>(lambda2));
+
+                /*Assign values*/
+                //calvingratex[iv]=ex*lambda1 - ey*lambda2;
+                //calvingratey[iv]=ey*lambda1 + ex*lambda2;
+                calvingratex[iv]=vx/vel*(lambda1 + lambda2);
+                calvingratey[iv]=vy/vel*(lambda1 + lambda2);
+                calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+        }
+
+        /*Add input*/
+        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+        /*Clean up and return*/
+        delete gauss;
 }
 /*}}}*/