Changeset 10278

Timestamp:

10/24/11 16:05:13 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added MacAyeal for Bal H

File:

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

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

@@ -5289 +5289 @@
         h=sqrt(2*this->GetArea());
 
-        ///*Create Artificial diffusivity once for all if requested*/
-        //if(stabilization){
-        //      gauss=new GaussTria();
-        //      gauss->GaussCenter();
-        //      GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
-        //      delete gauss;
-
-        //      vxaverage_input->GetInputAverage(&vx);
-        //      vyaverage_input->GetInputAverage(&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:*/
         gauss=new GaussTria(2);
@@ -5338 +5325 @@
 
                 if(stabilization==1){
+                        /*Streamline upwinding*/
                         vel=sqrt(pow(vx,2.)+pow(vy,2.));
                         K[0][0]=h/(2*vel)*vx*vx;
@@ -5343 +5331 @@
                         K[0][1]=h/(2*vel)*vx*vy;
                         K[1][1]=h/(2*vel)*vy*vy;
+                }
+                else if(stabilization==2){
+                        /*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];
-
-                        //KDL[0][0]=DL_scalar*K[0][0];
-                        //KDL[1][1]=DL_scalar*K[1][1];
-
                         TripleMultiply( &Bprime[0][0],2,numdof,1,
                                                 &KDL[0][0],2,2,0,