Changeset 26830 for issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

Timestamp:

01/27/22 18:17:58 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added sigmavm for Julia

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (6 diffs)

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

@@ -271 +271 @@
 
         if(!this->IsOnBase()) return;
-
-        IssmDouble  xyz_list[NUMVERTICES][3];
-        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+        this->ComputeSigmaVM();
+
         IssmDouble  calvingratex[NUMVERTICES];
         IssmDouble  calvingratey[NUMVERTICES];
         IssmDouble  calvingrate[NUMVERTICES];
-        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-        IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
-        IssmDouble  epse_2,groundedice,bed,sealevel;
-        IssmDouble  sigma_vm[NUMVERTICES];
-
-        /* Get node coordinates and dof list: */
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        IssmDouble  vx,vy;
+        IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
+        IssmDouble  groundedice,bed,sealevel;
 
         /*Depth average B for stress calculation*/
-        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
@@ -295 +289 @@
         Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
         Input* bs_input = this->GetInput(BaseEnum);                    _assert_(bs_input);
-        Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
-        Input* n_input  = this->GetInput(MaterialsRheologyNEnum);   _assert_(n_input);
         Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
         Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
         Input* sl_input  = this->GetInput(SealevelEnum); _assert_(sl_input);
+        Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input);
 
         /* Start looping on the number of vertices: */
@@ -307 +300 @@
 
                 /*Get velocity components and thickness*/
-                B_input->GetInputValue(&B,&gauss);
-                n_input->GetInputValue(&n,&gauss);
                 vx_input->GetInputValue(&vx,&gauss);
                 vy_input->GetInputValue(&vy,&gauss);
+                sigma_vm_input->GetInputValue(&sigma_vm,&gauss);
                 gr_input->GetInputValue(&groundedice,&gauss);
                 bs_input->GetInputValue(&bed,&gauss);
                 smax_fl_input->GetInputValue(&sigma_max_floating,&gauss);
                 smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss);
-                vel=sqrt(vx*vx+vy*vy)+1.e-14;
                 sl_input->GetInputValue(&sealevel,&gauss);
-
-                /*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]);
-                _assert_(!xIsNan<IssmDouble>(lambda1));
-                _assert_(!xIsNan<IssmDouble>(lambda2));
-
-                /*Process Eigen values (only account for extension)*/
-                lambda1 = max(lambda1,0.);
-                lambda2 = max(lambda2,0.);
-
-                /*Calculate sigma_vm*/
-                epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
-                sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n));
 
                 /*Tensile stress threshold*/
@@ -346 +321 @@
                 }
                 else{
-                        calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
-                        calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
+                        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]);
@@ -356 +331 @@
         this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
         this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
-        this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
 
         this->InputExtrude(CalvingratexEnum,-1);
         this->InputExtrude(CalvingrateyEnum,-1);
         this->InputExtrude(CalvingCalvingrateEnum,-1);
-        this->InputExtrude(SigmaVMEnum,-1);
 }
 /*}}}*/
@@ -847 +820 @@
         this->AddInput(StressTensoryzEnum,&sigma_yz[0],P1DGEnum);
         this->AddInput(StressTensorzzEnum,&sigma_zz[0],P1DGEnum);
+}
+/*}}}*/
+void       Penta::ComputeSigmaVM(){/*{{{*/
+
+        if(!this->IsOnBase()) return;
+
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
+        IssmDouble  B,n;
+        IssmDouble  sigma_vm[NUMVERTICES];
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+        /*Depth average B for stress calculation*/
+        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+        /*Retrieve all inputs and parameters we will need*/
+        Input* vx_input = this->GetInput(VxAverageEnum); _assert_(vx_input);
+        Input* vy_input = this->GetInput(VyAverageEnum); _assert_(vy_input);
+        Input* B_input  = this->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
+        Input* n_input  = this->GetInput(MaterialsRheologyNEnum);   _assert_(n_input);
+
+        /* Start looping on the number of vertices: */
+        GaussPenta gauss;
+        for(int iv=0;iv<3;iv++){
+                gauss.GaussVertex(iv);
+
+                /*Get velocity components and thickness*/
+                B_input->GetInputValue(&B,&gauss);
+                n_input->GetInputValue(&n,&gauss);
+                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]);
+                _assert_(!xIsNan<IssmDouble>(lambda1));
+                _assert_(!xIsNan<IssmDouble>(lambda2));
+
+                /*Process Eigen values (only account for extension)*/
+                lambda1 = max(lambda1,0.);
+                lambda2 = max(lambda2,0.);
+
+                /*Calculate sigma_vm*/
+                IssmDouble epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
+                sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n));
+        }
+
+        /*Add input*/
+        this->AddBasalInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum);
+        this->InputExtrude(SigmaVMEnum,-1);
 }
 /*}}}*/