Changeset 6268

Timestamp:

10/12/10 15:05:50 (14 years ago)

Author:

seroussi

Message:

added GetStabilizationParameter routine

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (5 diffs)
• Penta.h (modified) (1 diff)

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

@@ -2918 +2918 @@
         double     heatcapacity,thermalconductivity,dt;
         double     scalar_artdiff;
-        double     tau_parameter,diameter,normu;
+        double     tau_parameter,diameter;
         double     xyz_list[NUMVERTICES][3];
         double     B[3][numdof];
@@ -3038 +3038 @@
                         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
 
-                        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
-                        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
-                                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
-                        }
-                        else tau_parameter=diameter/(2*normu);
+                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+                        //if(dt) scalar_artdiff=scalar_artdiff*dt;
+
                         scalar_artdiff=tau_parameter*gauss->weight*Jdet;
-
                         for(i=0;i<numdof;i++){
                                 for(j=0;j<numdof;j++){
@@ -3830 +3827 @@
         double thermalconductivity;
         double viscosity,temperature;
-        double tau_parameter,diameter,normu;
+        double tau_parameter,diameter;
         double u,v,w,scalar_artdiff;
         double scalar_def,scalar_transient;
@@ -3890 +3887 @@
                         vz_input->GetParameterValue(&w, gauss);
 
-                        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
-                        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
-                                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
-                        }
-                        else tau_parameter=diameter/(2*normu);
-
+                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
                         scalar_artdiff=tau_parameter*gauss->weight*Jdet*phi/(rho_ice*heatcapacity);
+                        //if(dt) scalar_artdiff=scalar_artdiff*dt;
+
                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
                 }
@@ -4296 +4290 @@
         delete gauss;
         xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Penta::GetStabilizationParameter {{{1*/
+double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
+        /*Compute stabilization parameter*/
+
+        double normu;
+        double tau_parameter;
+
+        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
+        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
+                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
+        }
+        else tau_parameter=diameter/(2*normu);
+
+        return tau_parameter;
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/Penta.h

@@ -196 +196 @@
                 void      GetSolutionFromInputsDiagnosticVert(Vec solutiong);
                 void      GetSolutionFromInputsThermal(Vec solutiong);
+                double GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
                 void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
                 void    GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);