Changeset 9512

Timestamp:

08/31/11 11:00:59 (14 years ago)

Author:

Mathieu Morlighem

Message:

better artificial diffusivity for bal H

File:

• issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (3 diffs)

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

@@ -516 +516 @@
         int        artdiff;
         int        i,j,ig,dim;
-        double     Jdettria ,vx,vy,dvxdx,dvydy;
+        double     Jdettria,vx,vy,dvxdx,dvydy,vel,h;
         double     dvx[2],dvy[2];
         double     xyz_list[NUMVERTICES][3];
@@ -546 +546 @@
                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
         }
-
-        /*Create Artificial diffusivity once for all if requested*/
-        if(artdiff){
-                gauss=new GaussTria();
-                gauss->GaussCenter();
-                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
-                delete gauss;
-
-                vxaverage_input->GetParameterAverage(&vx);
-                vyaverage_input->GetParameterAverage(&vy);
-                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
-                K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
-        }
+        h=sqrt(2*this->GetArea());
+
+        ///*Create Artificial diffusivity once for all if requested*/
+        //if(artdiff){
+        //      gauss=new GaussTria();
+        //      gauss->GaussCenter();
+        //      GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
+        //      delete gauss;
+
+        //      vxaverage_input->GetParameterAverage(&vx);
+        //      vyaverage_input->GetParameterAverage(&vy);
+        //      K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
+        //      K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
+        //}
 
         /*Start looping on the number of gaussian points:*/
@@ -596 +597 @@
 
                 if(artdiff){
+                        vel=sqrt(pow(vx,2.)+pow(vy,2.));
+                        K[0][0]=h/(2*vel)*vx*vx;
+                        K[1][0]=h/(2*vel)*vy*vx;
+                        K[0][1]=h/(2*vel)*vx*vy;
+                        K[1][1]=h/(2*vel)*vy*vy;
                         KDL[0][0]=DL_scalar*K[0][0];
+                        KDL[1][0]=DL_scalar*K[1][0];
+                        KDL[0][1]=DL_scalar*K[0][1];
                         KDL[1][1]=DL_scalar*K[1][1];
+
+                        //KDL[0][0]=DL_scalar*K[0][0];
+                        //KDL[1][1]=DL_scalar*K[1][1];
 
                         TripleMultiply( &Bprime[0][0],2,numdof,1,