Changeset 16872 for issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

Timestamp:

11/21/13 16:46:58 (11 years ago)

Author:

seroussi

Message:

CHG: added SSA friction KMatrix in analysis

File:

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

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

@@ -974 +974 @@
         return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
+        /*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;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+        int numdof   = numnodes*2;
+
+        /*Initialize Element matrix and vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
+        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
+        IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
+
+        /*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);
+
+        /*build friction object, used later on: */
+        //friction=new Friction("2d",inputs,matpar,analysis_type);
+
+        /*Recover portion of element that is grounded*/
+        if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list);
+        if(migration_style==SubelementMigration2Enum){
+                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+        }
+        else{
+                gauss = element->NewGauss(2);
+        }
+
+        /* 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;
+                if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
+                if(migration_style==SubelementMigration2Enum){
+                        gllevelset_input->GetInputValue(&gllevelset, gauss);
+                        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;
+
+                TripleMultiply(B,2,numdof,1,
+                                        D,2,2,0,
+                                        B,2,numdof,0,
+                                        &Ke->values[0],1);
+        }
+
+        /*Transform Coordinate System*/
+        element->TransformStiffnessMatrixCoord(Ke,XYEnum);
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(D);
+        xDelete<IssmDouble>(B);
+        return Ke;
+}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
 
@@ -1034 +1115 @@
         xDelete<IssmDouble>(B);
         return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
-        return NULL;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
@@ -1207 +1285 @@
         /*Clean-up*/
         xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
+         * For node i, Bi can be expressed in the actual coordinate system
+         * by: 
+         *                 Bi=[ N   0 ]
+         *                    [ 0   N ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* basis=xNew<IssmDouble>(numnodes);
+        element->NodalFunctions(basis,gauss);
+
+        /*Build L: */
+        for(int i=0;i<numnodes;i++){
+                B[2*numnodes*0+2*i+0] = basis[i];
+                B[2*numnodes*0+2*i+1] = 0.;
+                B[2*numnodes*1+2*i+0] = 0.;
+                B[2*numnodes*1+2*i+1] = basis[i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(basis);
 }/*}}}*/
 void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/