Changeset 11262

Timestamp:

01/31/12 09:15:07 (13 years ago)

Author:

seroussi

Message:

started to work on enthalpy

File:

• issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (modified) (11 diffs)

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

@@ -3207 +3207 @@
         int        i,j,ig,found=0;
         double     Jdet,u,v,w,um,vm,wm;
+        double     h,hx,hy,hz,vx,vy,vz,vel;
         double     gravity,rho_ice,rho_water;
         double     epsvel=2.220446049250313e-16;
@@ -3214 +3215 @@
         double     tau_parameter,diameter;
         double     xyz_list[NUMVERTICES][3];
-        double     B[3][numdof];
-        double     Bprime[3][numdof];
         double     B_conduct[3][numdof];
         double     B_advec[3][numdof];
-        double     B_stab[2][numdof];
         double     Bprime_advec[3][numdof];
         double     L[numdof];
@@ -3225 +3223 @@
         double     D_scalar_trans,D_scalar_stab;
         double     D[3][3];
-        double     K[2][2]={0.0};
+        double     K[3][3]={0.0};
         Tria*      tria=NULL;
         GaussPenta *gauss=NULL;
@@ -3262 +3260 @@
                 /*Conduction: */  
                 /*Need to change that depending on enthalpy value -> cold or temperate ice: */  
-
                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
@@ -3281 +3278 @@
 
                 /*Advection: */
-
                 GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); 
                 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
@@ -3291 +3287 @@
                 vym_input->GetInputValue(&vm,gauss);
                 vzm_input->GetInputValue(&wm,gauss);
+                vx=u-um; vy=v-vm; vz=w-wm;
 
                 D_scalar_advec=gauss->weight*Jdet;
                 if(dt) D_scalar_advec=D_scalar_advec*dt;
 
-                D[0][0]=D_scalar_advec*(u-um);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);
+                D[0][0]=D_scalar_advec*vx;D[0][1]=0;                D[0][2]=0;
+                D[1][0]=0;                D[1][1]=D_scalar_advec*vy;D[1][2]=0;
+                D[2][0]=0;                D[2][1]=0;                D[2][2]=D_scalar_advec*vz;
 
                 TripleMultiply(&B_advec[0][0],3,numdof,1,
@@ -3305 +3302 @@
 
                 /*Transient: */
-
                 if(dt){
                         GetNodalFunctionsP1(&L[0], gauss);
@@ -3318 +3314 @@
 
                 /*Artifficial diffusivity*/
-
                 if(stabilization==1){
                         /*Build K: */
-                        D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
+                        GetElementSizes(&hx,&hy,&hz);
+                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
+                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
+                        K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
+                        K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
+                        K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
+                        D_scalar_stab=gauss->weight*Jdet;
                         if(dt) D_scalar_stab=D_scalar_stab*dt;
-                        K[0][0]=D_scalar_stab*pow(u,2);       K[0][1]=D_scalar_stab*fabs(u)*fabs(v);
-                        K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2);
-                        _error_("TO BE RECODED");
-
-                        //GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 
-
-                        TripleMultiply(&B_stab[0][0],2,numdof,1,
-                                                &K[0][0],2,2,0,
-                                                &B_stab[0][0],2,numdof,0,
+                        for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j];
+
+                        GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
+
+                        TripleMultiply(&Bprime_advec[0][0],3,numdof,1,
+                                                &K[0][0],3,3,0,
+                                                &Bprime_advec[0][0],3,numdof,0,
                                                 &Ke->values[0],1);
                 }
                 else if(stabilization==2){
-
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
-
                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
 
@@ -3456 +3453 @@
         double     tau_parameter,diameter;
         double     xyz_list[NUMVERTICES][3];
-        double     B[3][numdof];
-        double     Bprime[3][numdof];
         double     B_conduct[3][numdof];
         double     B_advec[3][numdof];
-        double     B_stab[2][numdof];
         double     Bprime_advec[3][numdof];
         double     L[numdof];
@@ -3537 +3531 @@
 
                 /*Transient: */
-
                 if(dt){
                         GetNodalFunctionsP1(&L[0], gauss);
@@ -3550 +3543 @@
 
                 /*Artifficial diffusivity*/
-
                 if(stabilization==1){
                         /*Build K: */