Changeset 18531

Timestamp:

09/19/14 14:25:30 (11 years ago)

Author:

seroussi

Message:

NEW: return SigmaNN as P1Enums instead of P0

File:

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

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

@@ -167 +167 @@
 
         if(!IsOnBase()){
-                IssmDouble sigma_nn=0;
-                this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn,P0Enum));
+                IssmDouble sigma_nn[3]={0.};
+                this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn[0],P1Enum));
                 return;
         }
@@ -174 +174 @@
                 IssmDouble* xyz_list=NULL;
                 IssmDouble *xyz_list_base=NULL;
-                IssmDouble  pressure,viscosity,sigma_nn;
+                IssmDouble  pressure,viscosity;
+                IssmDouble  sigma_nn[3];
                 IssmDouble  sigma_xx,sigma_xy,sigma_yy;
                 IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
@@ -192 +193 @@
 
                 /* Start looping on the number of vertices: */
-                Gauss* gauss=NewGaussBase(1);
-                gauss->GaussPoint(0);
-
-                /*Compute strain rate viscosity and pressure: */
-                this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-                this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
-                pressure_input->GetInputValue(&pressure,gauss);
-
-                /*Compute Stress*/
-                sigma_xx=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
-                sigma_yy=2*viscosity*epsilon[1]-pressure;
-                sigma_xy=2*viscosity*epsilon[2];
-
-                /*Get normal vector to the bed */
-                NormalBase(&base_normal[0],xyz_list_base);
-
-                /*Compute sigma_nn*/
-                sigma_nn=sigma_xx*base_normal[0]*base_normal[0] + 2*sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1];
+                Gauss* gauss = this->NewGauss();
+                for(int i=0;i<NUMVERTICES;i++){
+                        gauss->GaussNode(P1Enum,i);
+
+                        /*Compute strain rate viscosity and pressure: */
+                        this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+                        this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
+                        pressure_input->GetInputValue(&pressure,gauss);
+
+                        /*Compute Stress*/
+                        sigma_xx=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
+                        sigma_yy=2*viscosity*epsilon[1]-pressure;
+                        sigma_xy=2*viscosity*epsilon[2];
+
+                        /*Get normal vector to the bed */
+                        NormalBase(&base_normal[0],xyz_list_base);
+
+                        /*Compute sigma_nn*/
+                        sigma_nn[i]=sigma_xx*base_normal[0]*base_normal[0] + 2*sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1];
+                }
 
                 /*Add Stress tensor components into inputs*/
-                this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn,P0Enum));
+                this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn[0],P1Enum));
 
                 /*Clean up and return*/
@@ -1745 +1748 @@
                         xz_plane[1]=sin(theta);       xz_plane[4]=0.;  
                         xz_plane[2]=0.;               xz_plane[5]=1.;          
+
                         if(groundedice>=0){
                                 this->nodes[i]->DofInSSet(1); //vy