Changeset 6270

Timestamp:

10/12/10 15:32:12 (14 years ago)

Author:

seroussi

Message:

some minor changes in thermal

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (8 diffs)

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

@@ -2928 +2928 @@
         double     L[numdof];
         double     dh1dh6[3][6];
-        double     D_scalar;
+        double     D_scalar_conduct,D_scalar_advec;
+        double     D_scalar_trans,D_scalar_artdiff;
         double     D[3][3];
         double     K[2][2]={0.0};
@@ -2969 +2970 @@
                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
-                D_scalar=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
-                if(dt) D_scalar=D_scalar*dt;
-
-                D[0][0]=D_scalar; D[0][1]=0; D[0][2]=0;
-                D[1][0]=0; D[1][1]=D_scalar; D[1][2]=0;
-                D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
+                D_scalar_conduct=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
+                if(dt) D_scalar_conduct=D_scalar_conduct*dt;
+
+                D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0;
+                D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0;
+                D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct;
 
                 TripleMultiply(&B_conduct[0][0],3,numdof,1,
@@ -2990 +2991 @@
                 vz_input->GetParameterValue(&w, gauss);
 
-                D_scalar=gauss->weight*Jdet;
-                if(dt) D_scalar=D_scalar*dt;
-
-                D[0][0]=D_scalar*u;D[0][1]=0;         D[0][2]=0;
-                D[1][0]=0;         D[1][1]=D_scalar*v;D[1][2]=0;
-                D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
+                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;
 
                 TripleMultiply(&B_advec[0][0],3,numdof,1,
@@ -3006 +3007 @@
                 if(dt){
                         GetNodalFunctionsP1(&L[0], gauss);
-                        D_scalar=gauss->weight*Jdet;
-                        D_scalar=D_scalar;
+                        D_scalar_trans=gauss->weight*Jdet;
+                        D_scalar_trans=D_scalar_trans;
 
                         TripleMultiply(&L[0],numdof,1,0,
-                                                &D_scalar,1,1,0,
+                                                &D_scalar_trans,1,1,0,
                                                 &L[0],1,numdof,0,
                                                 &Ke_gaussian_transient[0][0],0);
@@ -3022 +3023 @@
                 if(artdiff==1){
                         /*Build K: */
-                        D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
-                        if(dt) D_scalar=D_scalar*dt;
-                        K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
-                        K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+                        D_scalar_artdiff=gauss->weight*Jdet/(pow(u,2)+pow(v,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);
+                        K[1][0]=D_scalar_artdiff*fabs(u)*fabs(v);K[1][1]=D_scalar_artdiff*pow(v,2);
 
                         GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
@@ -3039 +3040 @@
 
                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
-                        //if(dt) scalar_artdiff=scalar_artdiff*dt;
+                        if(dt) scalar_artdiff=scalar_artdiff*dt;
 
                         scalar_artdiff=tau_parameter*gauss->weight*Jdet;
@@ -3828 +3829 @@
         double viscosity,temperature;
         double tau_parameter,diameter;
-        double u,v,w,scalar_artdiff;
+        double u,v,w;
         double scalar_def,scalar_transient;
         double temperature_list[NUMVERTICES];
@@ -3888 +3889 @@
 
                         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]);
+
+                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
                 }
         }