Changeset 16940

Timestamp:

11/25/13 20:06:23 (11 years ago)

Author:

seroussi

Message:

CHG: Added Pvector coupling SSAFS viscous and friction

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (3 diffs)
• StressbalanceAnalysis.h (modified) (1 diff)

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

@@ -855 +855 @@
                 case HOFSApproximationEnum: 
                         return CreatePVectorHOFS(element);
+                case SSAFSApproximationEnum: 
+                        return CreatePVectorSSAFS(element);
                 case NoneApproximationEnum:
                         return NULL;
@@ -3356 +3358 @@
         return pe;
 }/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAFS(Element* element){/*{{{*/
+
+        /*compute all load vectors for this element*/
+        int init = element->FiniteElement();
+        element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+        ElementVector* pe1=CreatePVectorSSA(element);
+        element->SetTemporaryElementType(init); // P1 needed for HO
+        ElementVector* pe2=CreatePVectorFS(element);
+        int indices[3]={18,19,20};
+        element->SetTemporaryElementType(MINIcondensedEnum); // P1 needed for HO
+        ElementMatrix* Ke = CreateKMatrixFS(element);
+        element->SetTemporaryElementType(init); // P1 needed for HO
+        pe2->StaticCondensation(Ke,3,&indices[0]);
+        delete Ke;
+        ElementVector* pe3=CreatePVectorCouplingSSAFS(element);
+        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+
+        /*clean-up and return*/
+        delete pe1;
+        delete pe2;
+        delete pe3;
+        return pe;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
 
@@ -3543 +3568 @@
         return pe;
 }/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFS(Element* element){/*{{{*/
+
+        /*compute all load vectors for this element*/
+        ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(element);
+        ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(element);
+        ElementVector* pe =new ElementVector(pe1,pe2);
+
+        /*clean-up and return*/
+        delete pe1;
+        delete pe2;
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSFriction(Element* element){/*{{{*/
+
+        /*Intermediaries*/
+        int         i,j,approximation;
+        int         dim=3;
+        IssmDouble  Jdet,Jdet2d,FSreconditioning;
+        IssmDouble      bed_normal[3];
+        IssmDouble  viscosity, w, alpha2_gauss;
+        IssmDouble  dw[3];
+        IssmDouble  basis[6]; //for the six nodes of the penta
+        IssmDouble      *xyz_list_tria = NULL;
+        IssmDouble  *xyz_list      = NULL;
+
+        /*Initialize Element vector and return if necessary*/
+        if(!element->IsOnBed() || element->IsFloating()) return NULL;
+        element->GetInputValue(&approximation,ApproximationEnum);
+        if(approximation!=SSAFSApproximationEnum) return NULL;
+        int vnumnodes = element->GetNumberOfNodesVelocity();
+        int pnumnodes = element->GetNumberOfNodesPressure();
+
+        /*Prepare coordinate system list*/
+        int* cs_list     = xNew<int>(vnumnodes+pnumnodes);
+        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
+        for(i=0;i<vnumnodes;i++){
+                cs_list[i]             = XYZEnum;
+                node_list[i]           = element->GetNode(i);
+        }
+        for(i=0;i<pnumnodes;i++){
+                cs_list[vnumnodes+i]   = PressureEnum;
+                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
+        }
+        ElementVector* pe=element->NewElementVector(FSvelocityEnum);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GetVerticesCoordinatesBase(&xyz_list_tria);
+        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        Input* vx_input=   element->GetInput(VxEnum);    _assert_(vx_input);
+        Input* vy_input=   element->GetInput(VyEnum);    _assert_(vy_input);
+        Input* vz_input=   element->GetInput(VzEnum);    _assert_(vz_input);
+        Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input);
+
+        /*build friction object, used later on: */
+        Friction* friction=new Friction(element,3);
+
+        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
+        Gauss* gauss=element->NewGaussBase(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
+                element->NodalFunctionsP1(basis, gauss);
+
+                vzSSA_input->GetInputValue(&w, gauss);
+                vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
+
+                element->NormalBase(&bed_normal[0],xyz_list_tria);
+                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+                friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
+
+                for(i=0;i<3;i++){
+                        pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
+                        pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
+                        pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
+                }
+        }
+
+        /*Transform coordinate system*/
+        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
+
+        /*Clean up and return*/
+        xDelete<int>(cs_list);
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(xyz_list_tria);
+        xDelete<Node*>(node_list);
+        delete gauss;
+        delete friction;
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSViscous(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        int         i,approximation;
+        IssmDouble  viscosity,Jdet,FSreconditioning;
+        IssmDouble  dw[3];
+        IssmDouble  *xyz_list = NULL;
+        IssmDouble  basis[6]; //for the six nodes of the penta
+        IssmDouble  dbasis[3][6]; //for the six nodes of the penta
+
+        /*Initialize Element vector and return if necessary*/
+        element->GetInputValue(&approximation,ApproximationEnum);
+        if(approximation!=SSAFSApproximationEnum) return NULL;
+        int vnumnodes = element->GetNumberOfNodesVelocity();
+        int pnumnodes = element->GetNumberOfNodesPressure();
+
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+        Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
+        for(i=0;i<vnumnodes;i++){
+                cs_list[i]             = XYZEnum;
+                node_list[i]           = element->GetNode(i);
+        }
+        for(i=0;i<pnumnodes;i++){
+                cs_list[vnumnodes+i]   = PressureEnum;
+                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
+        }
+        ElementVector* pe=element->NewElementVector(FSvelocityEnum);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        Input* vx_input   =element->GetInput(VxEnum);      _assert_(vx_input);
+        Input* vy_input   =element->GetInput(VyEnum);      _assert_(vy_input);
+        Input* vz_input   =element->GetInput(VzEnum);      _assert_(vz_input);
+        Input* vzSSA_input=element->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(5);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctionsP1(&basis[0], gauss);
+                element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+
+                vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
+                element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
+
+                for(i=0;i<6;i++){
+                        pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
+                        pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
+                        pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
+                        pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
+                }
+        }
+
+        /*Transform coordinate system*/
+        element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
+
+        /*Clean up and return*/
+        xDelete<int>(cs_list);
+        delete gauss;
+        return pe;
+}/*}}}*/
 void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. 

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

@@ -78 +78 @@
                 ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
                 ElementVector* CreatePVectorSSAHO(Element* element);
+                ElementVector* CreatePVectorSSAFS(Element* element);
+                ElementVector* CreatePVectorCouplingSSAFS(Element* element);
+                ElementVector* CreatePVectorCouplingSSAFSFriction(Element* element);
+                ElementVector* CreatePVectorCouplingSSAFSViscous(Element* element);
                 ElementVector* CreatePVectorHOFS(Element* element);
                 ElementVector* CreatePVectorCouplingHOFS(Element* element);