Changeset 16018

Timestamp:

08/29/13 13:56:40 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added stabilization

Location:

issm/trunk-jpl/src/c

Files:

• analyses/balancevelocity_core.cpp (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp

@@ -19 +19 @@
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
-        if(VerboseSolution()) _printf0_("call computational core:\n");
-        //femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
-        //solutionsequence_linear(femmodel);
-        //femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
-        //solutionsequence_linear(femmodel);
-        surfaceslope_core(femmodel);
+        if(VerboseSolution()) _printf0_("computing smoothed slopes:\n");
+        femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
+        solutionsequence_linear(femmodel);
+        femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
+        solutionsequence_linear(femmodel);
+        //surfaceslope_core(femmodel);
 
         if(VerboseSolution()) _printf0_("call computational core:\n");

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -7231 +7231 @@
         IssmDouble h,gamma,thickness;
         IssmDouble hnx,hny,dhnx[2],dhny[2];
-        int        i,j;
-        IssmDouble D_scalar;
+        IssmDouble vel,vx,vy,dvxdx,dvydy;
+        IssmDouble D[2][2];
+        int        stabilization=1;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -7239 +7240 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
         IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
@@ -7295 +7297 @@
                                                         );
                         }
+                }
+
+                vx=hnx;
+                vy=hny;
+                if(stabilization==2){
+                        /*Streamline upwinding*/
+                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                        D[0][0]=h/(2*vel)*vx*vx;
+                        D[1][0]=h/(2*vel)*vy*vx;
+                        D[0][1]=h/(2*vel)*vx*vy;
+                        D[1][1]=h/(2*vel)*vy*vy;
+                }
+                else if(stabilization==1){
+                        /*MacAyeal*/
+                        D[0][0]=h/2.0*fabs(vx);
+                        D[0][1]=0.;
+                        D[1][0]=0.;
+                        D[1][1]=h/2.0*fabs(vy);
+                }
+                if(stabilization==1 || stabilization==2){
+                        GetBprimeMasstransport(B,&xyz_list[0][0],gauss);
+                        D[0][0]=gauss->weight*Jdettria*D[0][0];
+                        D[1][0]=gauss->weight*Jdettria*D[1][0];
+                        D[0][1]=gauss->weight*Jdettria*D[0][1];
+                        D[1][1]=gauss->weight*Jdettria*D[1][1];
+                        TripleMultiply(B,2,numnodes,1,
+                                                &D[0][0],2,2,0,
+                                                B,2,numnodes,0,
+                                                &Ke->values[0],1);
                 }
         }
@@ -7306 +7337 @@
         xDelete<IssmDouble>(HNx);
         xDelete<IssmDouble>(HNy);
+        xDelete<IssmDouble>(B);
         delete gauss;
         return Ke;
@@ -7317 +7349 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble D[2][2];
-        IssmDouble l=12.;
+        IssmDouble l=15.;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -7339 +7371 @@
                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
+                if(thickness<50.) thickness=50.;
 
                 D_scalar=gauss->weight*Jdet*(l*thickness)*(l*thickness);