Changeset 6212

Timestamp:

10/08/10 16:34:31 (14 years ago)

Author:

seroussi

Message:

added Franca2003 artificial diffusivity

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (11 diffs)
• Penta.h (modified) (1 diff)
• Tria.cpp (modified) (3 diffs)

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

@@ -2858 +2858 @@
 
         /*Intermediaries */
-        bool       artdiff;
+        int        artdiff;
         int        i,j,ig,found=0;
         double     Jdet,u,v,w,epsvel;
         double     gravity,rho_ice,rho_water;
         double     heatcapacity,thermalconductivity,dt;
+        double     scalar_artdiff;
+        double     tau_parameter,diameter,normu;
         double     xyz_list[NUMVERTICES][3];
         double     B[3][numdof];
@@ -2871 +2873 @@
         double     Bprime_advec[3][numdof];
         double     L[numdof];
+        double     dh1dh6[3][6];
         double     D_scalar;
         double     D[3][3];
@@ -2903 +2906 @@
         Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
         Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
+        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -2961 +2965 @@
                 /*Artifficial diffusivity*/
 
-                if(artdiff){
+                if(artdiff==1){
                         /*Build K: */
                         D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
@@ -2972 +2976 @@
                         MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
                         MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
+                }
+                else if(artdiff==2){
+
+                        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);
+                        scalar_artdiff=tau_parameter*gauss->weight*Jdet;
+
+                        for(i=0;i<numdof;i++){
+                                for(j=0;j<numdof;j++){
+                                        Ke_gaussian_artdiff[i][j]=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]);
+                                }
+                        }
                 }
                 else{
@@ -3746 +3767 @@
         /*Intermediaries*/
         int    i,j,ig,found=0;
-        int    friction_type;
+        int    friction_type,artdiff;
         double Jdet,phi,dt;
         double rho_ice,heatcapacity;
+        double thermalconductivity;
         double viscosity,temperature;
+        double tau_parameter,diameter,normu;
+        double u,v,w,scalar_artdiff;
         double scalar_def,scalar_transient;
         double temperature_list[NUMVERTICES];
         double xyz_list[NUMVERTICES][3];
         double L[numdof];
+        double dh1dh6[3][6];
         double epsilon[6];
         GaussPenta *gauss=NULL;
@@ -3763 +3788 @@
         /*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        this->parameters->FindParam(&dt,DtEnum);
         rho_ice=matpar->GetRhoIce();
         heatcapacity=matpar->GetHeatCapacity();
+        thermalconductivity=matpar->GetThermalConductivity();
+        this->parameters->FindParam(&dt,DtEnum);
+        this->parameters->FindParam(&artdiff,ArtDiffEnum);
         Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
@@ -3771 +3798 @@
         Input* temperature_input=NULL;
         if (dt) temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs);
+        if (artdiff==2) diameter=MinEdgeLength(xyz_list);
 
         /* Start  looping on the number of gaussian points: */
@@ -3795 +3823 @@
                         scalar_transient=temperature*Jdet*gauss->weight;
                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
+                }
+
+                if(artdiff==2){
+                        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+
+                        vx_input->GetParameterValue(&u, gauss);
+                        vy_input->GetParameterValue(&v, gauss);
+                        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);
+
+                        scalar_artdiff=tau_parameter*gauss->weight*Jdet*phi/(rho_ice*heatcapacity);
+                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
                 }
         }
@@ -5441 +5486 @@
 }
 /*}}}*/
+/*FUNCTION Penta::MinEdgeLength{{{1*/
+double Penta::MinEdgeLength(double xyz_list[6][3]){
+        /*Return the minimum lenght of the nine egdes of the penta*/
+
+        int    i,node0,node1;
+        int    edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges
+        double minlength=-1;
+        double length;
+        
+        for(i=0;i<9;i++){
+                /*Find the two nodes for this edge*/
+                node0=edges[i][0];
+                node1=edges[i][1];
+
+                /*Compute the length of this edge and compare it to the minimal length*/
+                length=pow(pow(xyz_list[node0][0]-xyz_list[node1][0],2.0)+pow(xyz_list[node0][1]-xyz_list[node1][1],2.0)+pow(xyz_list[node0][2]-xyz_list[node1][2],2.0),0.5);
+                if(length<minlength || minlength<0) minlength=length;
+        }
+
+        return minlength;
+}
+/*}}}*/
 /*FUNCTION Penta::ReduceMatrixStokes {{{1*/
 void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
@@ -5604 +5671 @@
         *(surface_normal+1)=normal[1]/normal_norm;
         *(surface_normal+2)=normal[2]/normal_norm;
-
-}
-/*}}}*/
+}
+/*}}}*/

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

@@ -218 +218 @@
                 bool    IsOnShelf(void); 
                 bool    IsOnWater(void); 
-                void      ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
+                double  MinEdgeLength(double xyz_list[6][3]);
+           void   ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
                 void      ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
                 void      SetClone(int* minranks);

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

@@ -2415 +2415 @@
 
         /*Intermediaries */
-        bool       artdiff;
+        int        artdiff;
         int        i,j,ig,dim;
         double     Jdettria ,vx,vy,dvxdx,dvydy;
@@ -2586 +2586 @@
 
         /*Intermediaries */
-        bool    artdiff;
+        int     artdiff;
         int     i,j,ig,dim;
         double  nx,ny,norm,Jdettria;
@@ -3157 +3157 @@
 
         /*Intermediaries */
-        bool       artdiff;
+        int        artdiff;
         int        i,j,ig,dim;
         double     Jdettria,DL_scalar,dt;