Changeset 16907

Timestamp:

11/23/13 19:32:10 (11 years ago)

Author:

seroussi

Message:

CHG: moving coupling SSAHO in analysis

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)
• analyses/StressbalanceAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (2 diffs)
• classes/Elements/Penta.cpp (modified) (2 diffs)
• classes/Elements/Penta.h (modified) (2 diffs)
• classes/Elements/Seg.h (modified) (2 diffs)
• classes/Elements/Tria.h (modified) (2 diffs)

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

@@ -826 +826 @@
                 case FSApproximationEnum: 
                         return CreateKMatrixFS(element);
+                case SSAHOApproximationEnum: 
+                        return CreateKMatrixSSAHO(element);
                 case NoneApproximationEnum:
                         return NULL;
@@ -847 +849 @@
                 case FSApproximationEnum: 
                         return CreatePVectorFS(element);
+                case SSAHOApproximationEnum: 
+                        return CreatePVectorSSAHO(element);
                 case NoneApproximationEnum:
                         return NULL;
@@ -3059 +3063 @@
 
 /*Coupling (Tiling)*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=CreateKMatrixSSA(element);
+        ElementMatrix* Ke2=CreateKMatrixHO(element);
+        ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element);
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
+
+        /*clean-up and return*/
+        delete Ke1;
+        delete Ke2;
+        delete Ke3;
+        return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element);
+        ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element);
+        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+
+        /*clean-up and return*/
+        delete Ke1;
+        delete Ke2;
+        return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
+        return NULL;
+
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
+
+        /*Constants*/
+        int numnodes = 2*element->GetNumberOfNodes();
+        int numdofm     = 1 *numnodes; //*2/2
+        int numdofp     = 2 *numnodes;
+        int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA
+
+        /*Intermediaries */
+        int         i,j;
+        IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
+        IssmDouble  *xyz_list      = NULL;
+        IssmDouble* B              = xNew<IssmDouble>(3*numdofp);
+        IssmDouble* Bprime         = xNew<IssmDouble>(3*numdofm);
+        IssmDouble  D[3][3]={0.0}; // material matrix, simple scalar matrix.
+        IssmDouble  D_scalar;
+        IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
+        IssmDouble* Ke_gg_gaussian = xNew<IssmDouble>(numdofp*numdofm);
+        Node        *node_list[numnodes];
+        int         cs_list[numnodes];
+
+        /*Find penta on bed as HO must be coupled to the dofs on the bed: */
+        Element* pentabase=element->GetBasalElement();
+        Element* basaltria=pentabase->SpawnBasalElement();
+
+        /*prepare node list*/
+        for(i=0;i<numnodes;i++){
+                node_list[i+0*numnodes] = pentabase->GetNode(i);
+                node_list[i+1*numnodes] = element  ->GetNode(i);
+                cs_list[i+0*numnodes] = XYEnum;
+                cs_list[i+1*numnodes] = XYEnum;
+        }
+
+        /*Initialize Element matrix*/
+        ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum);
+        ElementMatrix* Ke2= element  ->NewElementMatrix(HOApproximationEnum);
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+        delete Ke1; delete Ke2;
+
+        /* Get node coordinates and dof list: */
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+        Input* vx_input   =element->GetInput(VxEnum);       _assert_(vx_input);
+        Input* vy_input   =element->GetInput(VyEnum);       _assert_(vy_input);
+        Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);
+        Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(5);
+        Gauss* gauss_tria=new GaussTria();
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                gauss->SynchronizeGaussBase(gauss_tria);
+
+                element->JacobianDeterminant(&Jdet, xyz_list,gauss);
+                this->GetBSSAHO(B, element,xyz_list, gauss);
+                //basaltria->GetBprimeSSA(Bprime, xyz_list, gauss_tria); /FIXME
+                this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria); 
+                element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
+                element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
+
+                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+                D_scalar=2*newviscosity*gauss->weight*Jdet;
+                for (i=0;i<3;i++) D[i][i]=D_scalar;
+
+                TripleMultiply( B,3,numdofp,1,
+                                        &D[0][0],3,3,0,
+                                        Bprime,3,numdofm,0,
+                                        Ke_gg_gaussian,0);
+
+                for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofp+j]+=Ke_gg_gaussian[i*numdofp+j];
+        } 
+        for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofp+j];
+        for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofp+i];
+
+        /*Transform Coordinate System*/
+        element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+
+        /*Clean-up and return*/
+        basaltria->DeleteMaterials(); delete basaltria;
+        delete gauss;
+        delete gauss_tria;
+        return Ke;
+
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/
+
+        /*compute all load vectors for this element*/
+        ElementVector* pe1=CreatePVectorSSA(element);
+        ElementVector* pe2=CreatePVectorHO(element);
+        ElementVector* pe =new ElementVector(pe1,pe2);
+
+        /*clean-up and return*/
+        delete pe1;
+        delete pe2;
+        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. 
+         * For node i, Bi can be expressed in the actual coordinate system
+         * by: 
+         *       Bi=[ dh/dx          0      ]
+         *          [   0           dh/dy   ]
+         *          [ 1/2*dh/dy  1/2*dh/dx  ]
+         * where h is the interpolation function for node i.
+         *
+         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+         */
+
+        int numnodes = element->GetNumberOfNodes();
+        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+
+        /*Get dbasis in actual coordinate system: */
+        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B[2*numnodes*0+2*i+0] = dbasis[0*3+i];
+                B[2*numnodes*0+2*i+1] = 0.;
+                B[2*numnodes*1+2*i+0] = 0.;
+                B[2*numnodes*1+2*i+1] = dbasis[1*3+i];
+                B[2*numnodes*2+2*i+0] = .5*dbasis[1*3+i];
+                B[2*numnodes*2+2*i+1] = .5*dbasis[0*3+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
 

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

@@ -71 +71 @@
                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
                 /*Coupling*/
+                ElementMatrix* CreateKMatrixSSAHO(Element* element);
+                ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
+                ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);
+                ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
+                ElementVector* CreatePVectorSSAHO(Element* element);
+                void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
                 void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -84 +84 @@
                 virtual void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum)=0;
                 virtual void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list)=0;
+                virtual void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list)=0;
                 virtual Element* GetBasalElement(void)=0;
                 virtual void    GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
@@ -120 +121 @@
                 virtual void   GetInputValue(IssmDouble* pvalue,int enum_type)=0;
                 virtual void   GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type)=0;
+                virtual Node*  GetNode(int node_number)=0;
                 virtual void   GetVerticesCoordinates(IssmDouble** xyz_list)=0;
                 virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;

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

@@ -1326 +1326 @@
 }
 /*}}}*/
+/*FUNCTION Penta::GetNode(int node_number) {{{*/
+Node* Penta::GetNode(int node_number){
+        return this->nodes[node_number];
+}
+/*}}}*/
 /*FUNCTION Penta::GetMaterialInput(int inputenum) {{{*/
 Input* Penta::GetMaterialInput(int inputenum){
@@ -3229 +3234 @@
 
         ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum_list);
+
+}
+/*}}}*/
+/*FUNCTION Penta::TransformStiffnessMatrixCoord(ElementMatrix* pe,Node** nodes_list,int numnodesint* transformenum_list){{{*/
+void Penta::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){
+
+        ::TransformStiffnessMatrixCoord(Ke,nodes_list,numnodes,transformenum_list);
 
 }

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -244 +244 @@
                 void           GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
                 void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
+                Node*          GetNode(int node_number);
                 void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
                 void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
@@ -295 +296 @@
                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
+                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list);
                 void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
                 void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -155 +155 @@
                 void        GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type){_error_("not implemented yet");};
                 void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
+                Node*       GetNode(int node_number){_error_("Not implemented");};
                 IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
                 IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
@@ -226 +227 @@
                 void        TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
                 void        TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
+                void        TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){_error_("not implemented yet");};
                 void UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};
                 void UpdateConstraintsExtrudeFromTop(){_error_("not implemented");};

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

@@ -286 +286 @@
                 void           GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type);
                 void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
+                Node*          GetNode(int node_number){_error_("not implemented yet");};
                 void             InputChangeName(int enum_type,int enum_type_old);
                 void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
@@ -326 +327 @@
                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
+                void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){_error_("not implemented yet");};
                 void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
                 void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);