Changeset 171

Timestamp:

05/01/09 14:44:07 (16 years ago)

Author:

seroussi

Message:

fixed Pattyn to be similar in cielo and ice

Location:

issm/trunk/src/c/objects

Files:

• Penta.cpp (modified) (9 diffs)
• Tria.cpp (modified) (3 diffs)

issm/trunk/src/c/objects/Penta.cpp

@@ -341 +341 @@
         /* material data: */
         double viscosity; //viscosity
+        double oldviscosity; //viscosity
+        double newviscosity; //viscosity
 
         /* strain rate: */
         double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
 
         /* matrices: */
@@ -367 +370 @@
         double* velocity_param=NULL;
         double  vxvy_list[numgrids][2];
+        double* oldvelocity_param=NULL;
+        double  oldvxvy_list[numgrids][2];
         double  thickness;
 
@@ -403 +408 @@
                 /*recover extra inputs from users, at current convergence iteration: */
                 velocity_param=ParameterInputsRecover(inputs,"velocity"); 
+                oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity"); 
 
                 /* Get node coordinates and dof list: */
@@ -427 +433 @@
                 }
                         
+                /*Initialize oldvxvy_list: */
+                for(i=0;i<numgrids;i++){
+                        if(oldvelocity_param){
+                                oldvxvy_list[i][0]=oldvelocity_param[doflist[i*numberofdofspernode+0]];
+                                oldvxvy_list[i][1]=oldvelocity_param[doflist[i*numberofdofspernode+1]];
+                        }
+                        else{
+                                oldvxvy_list[i][0]=0;
+                                oldvxvy_list[i][1]=0;
+                        }
+                }
+
                 #ifdef _DEBUGELEMENTS_
                 if(my_rank==RANK && id==ELID){ 
@@ -486 +504 @@
                                 /*Get strain rate from velocity: */
                                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
+                                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
                         
                                 /*Get viscosity: */
                                 matice->GetViscosity3d(&viscosity, &epsilon[0]);
+                                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
 
                                 /*Get B and Bprime matrices: */
@@ -499 +519 @@
                                 /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
                                   onto this scalar matrix, so that we win some computational time: */
-                                D_scalar=viscosity*gauss_weight*Jdet;
+                                
+                                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+                                D_scalar=newviscosity*gauss_weight*Jdet;
                                 for (i=0;i<5;i++){
                                         D[i][i]=D_scalar;
@@ -510 +532 @@
                                                 &Ke_gg_gaussian[0][0],0);
 
-                                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+                                /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
                                 for( i=0; i<numdofs; i++){
                                         for(j=0;j<numdofs;j++){
@@ -594 +616 @@
         /*If this element is on the bed rock, spawn a Tria element out of the bed triangle, and use it 
          * to create the stifness matrix for the base velocity: */
-        if(onbed){
-                tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 for the base
-                tria->CreateKMatrixDiagnosticBaseVert(Kgg,inputs, analysis_type);
-                delete tria;
-        }
+        //if(onbed){
+//              tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 for the base
+//              tria->CreateKMatrixDiagnosticBaseVert(Kgg,inputs, analysis_type);
+//              delete tria;
+//      }
 
         /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness 
@@ -1678 +1700 @@
         /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the 
          *diagnostic base vertical stifness: */
-        tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
-        tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
-        delete tria;
+        if(onbed){
+                tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
+                tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
+                delete tria;
+        }
 
         /*Now, handle the standard penta element: */

issm/trunk/src/c/objects/Tria.cpp

@@ -663 +663 @@
                 GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
 
-                if (velocity_param){
-                        DL_scalar=alpha2*gauss_weight*Jdet;
-                }
-                else{
-                        DL_scalar=0;
-                }
-                        
+                DL_scalar=alpha2*gauss_weight*Jdet;
                 for (i=0;i<2;i++){
                         DL[i][i]=DL_scalar;
@@ -2126 +2120 @@
                                         &Ke_gg_gaussian[0][0],0);
 
+                
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
                 for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
@@ -2355 +2350 @@
                 /*Get bed slope: */
                 GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
-                dbdx=slope[0];
-                dbdy=slope[1];
+                dbdx=-slope[0];
+                dbdy=-slope[1];
 
                 /* Get Jacobian determinant: */