Changeset 10277

Timestamp:

10/24/11 16:04:31 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added SU

File:

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

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

@@ -501 +501 @@
         int        stabilization;
         int        i,j,ig,dim;
-        double     Jdettria,DL_scalar,dt;
-        double     vx,vy,dvxdx,dvydy;
+        double     Jdettria,DL_scalar,dt,h;
+        double     vel,vx,vy,dvxdx,dvydy;
         double     dvx[2],dvy[2];
         double     v_gauss[2]={0.0};
@@ -533 +533 @@
                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
         }
-
-        //Create Artificial diffusivity once for all if requested
-        if(stabilization==1){
-                gauss=new GaussTria();
-                gauss->GaussCenter();
-                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
-                delete gauss;
-
-                vxaverage_input->GetInputAverage(&v_gauss[0]);
-                vyaverage_input->GetInputAverage(&v_gauss[1]);
-
-                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
-                K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
-        }
+        h=sqrt(2*this->GetArea());
 
         /* Start  looping on the number of gaussian points: */
@@ -591 +578 @@
                                         &Ke->values[0],1);
 
-                if(stabilization==1){
+                if(stabilization==2){
+                        /*Streamline upwinding*/
+                        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;
+                }
+                else if(stabilization==1){
+                        /*MacAyeal*/
+                        vxaverage_input->GetInputAverage(&vx);
+                        vyaverage_input->GetInputAverage(&vy);
+                        K[0][0]=h/2.0*fabs(vx) /sqrt(2*sqrt(3)); // the second part should not be there
+                        K[0][1]=0.;
+                        K[1][0]=0.;
+                        K[1][1]=h/2.0*fabs(vy) /sqrt(2*sqrt(3)); // the second part should not be there
+                }
+                if(stabilization==1 || stabilization==2){
                         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];
-
                         TripleMultiply( &Bprime[0][0],2,numdof,1,
                                                 &KDL[0][0],2,2,0,
@@ -1078 +1083 @@
         x3=xyz_list[2][0]; y3=xyz_list[2][1];
  
+        _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
         return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
 }