Changeset 8117

Timestamp:

05/03/11 09:08:29 (14 years ago)

Author:

Mathieu Morlighem

Message:

Added mesh velocity for thermal model. Still need to initialize it

Location:

issm/trunk/src/c

Files:

• EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• EnumDefinitions/EnumToString.cpp (modified) (1 diff)
• EnumDefinitions/StringToEnum.cpp (modified) (1 diff)
• objects/Elements/Penta.cpp (modified) (7 diffs)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -436 +436 @@
         RheologyLawEnum,
         PatersonEnum,
-        ArrheniusEnum
-        /*}}}*/
+        ArrheniusEnum,
+        /*}}}*/
+        VxMeshEnum,
+        VyMeshEnum,
+        VzMeshEnum
 };
 

issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

@@ -384 +384 @@
                 case PatersonEnum : return "Paterson";
                 case ArrheniusEnum : return "Arrhenius";
+                case VxMeshEnum : return "VxMesh";
+                case VyMeshEnum : return "VyMesh";
+                case VzMeshEnum : return "VzMesh";
                 default : return "unknown";
 

issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

@@ -382 +382 @@
         else if (strcmp(name,"Paterson")==0) return PatersonEnum;
         else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
+        else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
+        else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
+        else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
         else _error_("Enum %s not found",name);
 

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

@@ -1703 +1703 @@
         int        artdiff;
         int        i,j,ig,found=0;
-        double     Jdet,u,v,w,epsvel;
+        double     Jdet,u,v,w,um,vm,wm,epsvel;
         double     gravity,rho_ice,rho_water;
         double     heatcapacity,thermalconductivity,dt;
@@ -1740 +1740 @@
         this->parameters->FindParam(&artdiff,ArtDiffEnum);
         this->parameters->FindParam(&epsvel,EpsVelEnum);
-        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);
+        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);
+        Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input);
+        Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input);
+        Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input);
         if (artdiff==2) diameter=MinEdgeLength(xyz_list);
 
@@ -1777 +1780 @@
                 vy_input->GetParameterValue(&v, gauss);
                 vz_input->GetParameterValue(&w, gauss);
+                vxm_input->GetParameterValue(&um,gauss);
+                vym_input->GetParameterValue(&vm,gauss);
+                vzm_input->GetParameterValue(&wm,gauss);
 
                 D_scalar_advec=gauss->weight*Jdet;
                 if(dt) D_scalar_advec=D_scalar_advec*dt;
 
-                D[0][0]=D_scalar_advec*u;D[0][1]=0;         D[0][2]=0;
-                D[1][0]=0;         D[1][1]=D_scalar_advec*v;D[1][2]=0;
-                D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar_advec*w;
+                D[0][0]=D_scalar_advec*(um-u);D[0][1]=0;                    D[0][2]=0;
+                D[1][0]=0;                    D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0;
+                D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*(w-wm);
 
                 TripleMultiply(&B_advec[0][0],3,numdof,1,
@@ -1810 +1816 @@
                 if(artdiff==1){
                         /*Build K: */
-                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
+                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
                         if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
                         K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
@@ -1826 +1832 @@
                         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
 
-                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
 
                         for(i=0;i<numdof;i++){
                                 for(j=0;j<numdof;j++){
-                                        Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]);
+                                        Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
+                                          ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
                                 }
                         }
@@ -1836 +1843 @@
                                 for(i=0;i<numdof;i++){
                                         for(j=0;j<numdof;j++){
-                                                Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
+                                                Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
                                         }
                                 }
@@ -6238 +6245 @@
                         break;
 
+                case ThermalAnalysisEnum:
+                        /*Initialize mesh velocity*/
+                        for(i=0;i<6;i++)nodeinputs[i]=0;
+                        this->inputs->AddInput(new PentaVertexInput(VxMeshEnum,nodeinputs));
+                        this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs));
+                        this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs));
+                        break;
+
                 default:
                         /*No update for other solution types*/