Changeset 70

Timestamp:

04/27/09 16:50:16 (16 years ago)

Author:

seroussi

Message:

added overshooting in tri

Location:

issm/trunk/src/c/objects

Files:

• Tria.cpp (modified) (13 diffs)
• Tria.h (modified) (2 diffs)

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

@@ -33 +33 @@
 
 Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],
-                                int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel){
+                                int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,double tria_viscosity_overshoot){
         
         int i;
@@ -60 +60 @@
         epsvel=tria_epsvel;
         onbed=1;
+        viscosity_overshoot=tria_viscosity_overshoot;
         
         return;
@@ -88 +89 @@
         printf("   epsvel: %g\n",epsvel);
         printf("   onbed: %i\n",onbed);
+        printf("   viscosity_overshoot=%g\n",viscosity_overshoot);
         printf("   nodes: \n");
         if(nodes[0])nodes[0]->Echo();
@@ -134 +136 @@
         memcpy(marshalled_dataset,&meanvel,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
         memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
+        memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
         
         *pmarshalled_dataset=marshalled_dataset;
@@ -161 +164 @@
                 +sizeof(meanvel)
                 +sizeof(epsvel)
+                +sizeof(viscosity_overshoot)
                 +sizeof(int); //sizeof(int) for enum type
 }
@@ -200 +204 @@
         memcpy(&meanvel,marshalled_dataset,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
         memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
+        memcpy(&viscosity_overshoot,marshalled_dataset,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
 
         /*nodes and materials are not pointing to correct objects anymore:*/
@@ -292 +297 @@
         /* material data: */
         double viscosity; //viscosity
+        double newviscosity; //viscosity
+        double oldviscosity; //viscosity
         
         /* strain rate: */
         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+        double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/
 
         /* matrices: */
@@ -318 +326 @@
         /*input parameters for structural analysis (diagnostic): */
         double* velocity_param=NULL;
+        double* oldvelocity_param=NULL;
         double  vxvy_list[numgrids][2];
+        double  oldvxvy_list[numgrids][2];
         double  thickness;
 
@@ -331 +341 @@
         /*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: */
@@ -338 +349 @@
         /* Set Ke_gg to 0: */
         for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
-
 
         /*Initialize vxvy_list: */
@@ -349 +359 @@
                         vxvy_list[i][0]=0;
                         vxvy_list[i][1]=0;
+                }
+
+                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;
                 }
         }
@@ -443 +462 @@
                 /*Get strain rate from velocity: */
                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
+                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
                 
                 /*Get viscosity: */
                 matice->GetViscosity2d(&viscosity, &epsilon[0]);
+                matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
                 
                 /* Get Jacobian determinant: */
@@ -452 +473 @@
                 /* 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*thickness*gauss_weight*Jdet;
+                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+                D_scalar=newviscosity*thickness*gauss_weight*Jdet;
+
                 for (i=0;i<3;i++){
                         D[i][i]=D_scalar;

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

@@ -43 +43 @@
                 double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
                 int    onbed;
+                double viscosity_overshoot;
         
         public:
@@ -48 +49 @@
                 Tria();
                 Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],
-                                int friction_type,double p,double q,int shelf,double meanvel,double epsvel);
+                                int friction_type,double p,double q,int shelf,double meanvel,double epsvel,double viscosity_overshoot);
                 ~Tria();