Changeset 16874

Timestamp:

11/21/13 19:01:27 (11 years ago)

Author:

seroussi

Message:

NEW: KMatrix friction SSA

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)

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

@@ -975 +975 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
+
+        /*Return if element is inactive*/
+        if(element->IsFloating()) return NULL;
+
         /*Intermediaries*/
         bool       mainlyfloating;
-        int        analysis_type,migration_style,point1;
-        IssmDouble alpha2,Jdet,fraction1,fraction2;
-        IssmDouble phi=1.0;
-        IssmDouble D_scalar,gllevelset;
-        IssmDouble *xyz_list = NULL;
-        Friction   *friction = NULL;
-        Gauss*     gauss     = NULL;
-
-        /*Return if element is inactive*/
-        if(element->IsFloating()) return NULL;
+        int         migration_style,point1;
+        IssmDouble  alpha2,Jdet,fraction1,fraction2;
+        IssmDouble  gllevelset,phi=1.;
+        IssmDouble *xyz_list  = NULL;
+        Gauss*      gauss     = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -993 +992 @@
 
         /*Initialize Element matrix and vectors*/
-        ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
-        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
-        IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
+        ElementMatrix* Ke      = element->NewElementMatrix(SSAApproximationEnum);
+        IssmDouble*    B       = xNew<IssmDouble>(2*numdof);
+        IssmDouble     D[2][2] = {0.};
 
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
         element->FindParam(&migration_style,GroundinglineMigrationEnum);
-        Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
-        Input* vx_input=element->GetInput(VxEnum);           _assert_(vx_input);
-        Input* vy_input=element->GetInput(VyEnum);           _assert_(vy_input);
-        Input* gllevelset_input=NULL;
-        element->FindParam(&analysis_type,AnalysisTypeEnum);
+        Input* surface_input    = element->GetInput(SurfaceEnum); _assert_(surface_input);
+        Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
+        Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
+        Input* gllevelset_input = NULL;
 
         /*build friction object, used later on: */
-        //friction=new Friction("2d",inputs,matpar,analysis_type);
+        Friction* friction=new Friction(element,2);
 
         /*Recover portion of element that is grounded*/
@@ -1022 +1020 @@
         /* Start  looping on the number of gaussian points: */
         for(int ig=gauss->begin();ig<gauss->end();ig++){
-
                 gauss->GaussPoint(ig);
 
-                //friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
-                alpha2=1.0;
+                friction->GetAlpha2(&alpha2, gauss,vx_input,vy_input,NULL);
                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
                 if(migration_style==SubelementMigration2Enum){
                         gllevelset_input->GetInputValue(&gllevelset, gauss);
-                        if(gllevelset<0) alpha2=0;
+                        if(gllevelset<0.) alpha2=0.;
                 }
 
                 this->GetBSSAFriction(B, element, xyz_list, gauss);
                 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
-                D_scalar=alpha2*gauss->weight*Jdet;
-                for(int i=0;i<2;i++) D[i*2+i]=D_scalar;
+                for(int i=0;i<2;i++) D[i][i]=alpha2*gauss->weight*Jdet;
 
                 TripleMultiply(B,2,numdof,1,
-                                        D,2,2,0,
+                                        &D[0][0],2,2,0,
                                         B,2,numdof,0,
                                         &Ke->values[0],1);
@@ -1051 +1046 @@
         delete friction;
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(D);
         xDelete<IssmDouble>(B);
         return Ke;
@@ -1252 +1246 @@
         delete gauss;
         return pe;
-        return NULL;
 }/*}}}*/
 void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/