Changeset 17812

Timestamp:

04/22/14 12:19:15 (11 years ago)

Author:

jbondzio

Message:

ADD: Compute Strainrates for Penta and Tria

Location:

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

Files:

• Penta.cpp (modified) (1 diff)
• Tria.cpp (modified) (1 diff)

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

@@ -279 +279 @@
 void  Penta::ComputeStrainRate(Vector<IssmDouble>* eps){
 
-        _error_("Not implemented yet");
-
+        IssmDouble      xyz_list[NUMVERTICES][3];
+        IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+        IssmDouble      eps_xx[NUMVERTICES];
+        IssmDouble               eps_yy[NUMVERTICES];
+        IssmDouble               eps_zz[NUMVERTICES];
+        IssmDouble      eps_xy[NUMVERTICES];
+        IssmDouble               eps_xz[NUMVERTICES];
+        IssmDouble               eps_yz[NUMVERTICES];
+        GaussPenta*     gauss=NULL;
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+        /*Retrieve all inputs we will be needing: */
+        Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);             _assert_(vz_input);
+
+        /* Start looping on the number of vertices: */
+        gauss=new GaussPenta();
+        for (int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Compute strain rate viscosity and pressure: */
+                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+
+                eps_xx[iv]=epsilon[0]; 
+                eps_yy[iv]=epsilon[1];
+                eps_zz[iv]=epsilon[2];
+                eps_xy[iv]=epsilon[3];
+                eps_xz[iv]=epsilon[4];
+                eps_yz[iv]=epsilon[5];
+        }
+
+        /*Add Stress tensor components into inputs*/
+        this->inputs->AddInput(new PentaInput(StrainRatexxEnum,&eps_xx[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(StrainRatexyEnum,&eps_xy[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(StrainRatexzEnum,&eps_xz[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(StrainRateyyEnum,&eps_yy[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(StrainRateyzEnum,&eps_yz[0],P1Enum));
+        this->inputs->AddInput(new PentaInput(StrainRatezzEnum,&eps_zz[0],P1Enum));
+
+        /*Clean up and return*/
+        delete gauss;
 }
 /*}}}*/

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

@@ -204 +204 @@
 /*FUNCTION Tria::ComputeStrainRate {{{*/
 void  Tria::ComputeStrainRate(Vector<IssmDouble>* eps){
-        _error_("Not Implemented yet");
+
+        IssmDouble      xyz_list[NUMVERTICES][3];
+        IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+        IssmDouble      eps_xx[NUMVERTICES];
+        IssmDouble               eps_yy[NUMVERTICES];
+        IssmDouble               eps_zz[NUMVERTICES]={0,0,0};
+        IssmDouble      eps_xy[NUMVERTICES];
+        IssmDouble               eps_xz[NUMVERTICES]={0,0,0};
+        IssmDouble               eps_yz[NUMVERTICES]={0,0,0};
+        GaussPenta*     gauss=NULL;
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+        /*Retrieve all inputs we will be needing: */
+        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: */
+        gauss=new GaussPenta();
+        for (int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Compute strain rate viscosity and pressure: */
+                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+
+                eps_xx[iv]=epsilon[0]; 
+                eps_yy[iv]=epsilon[1];
+                eps_xy[iv]=epsilon[2];
+        }
+
+        /*Add Stress tensor components into inputs*/
+        this->inputs->AddInput(new TriaInput(StrainRatexxEnum,&eps_xx[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(StrainRatexyEnum,&eps_xy[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(StrainRatexzEnum,&eps_xz[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(StrainRateyyEnum,&eps_yy[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(StrainRateyzEnum,&eps_yz[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(StrainRatezzEnum,&eps_zz[0],P1Enum));
+
+        /*Clean up and return*/
+        delete gauss;
+
 }
 /*}}}*/