Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16906) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16907) @@ -826,4 +826,6 @@ case FSApproximationEnum: return CreateKMatrixFS(element); + case SSAHOApproximationEnum: + return CreateKMatrixSSAHO(element); case NoneApproximationEnum: return NULL; @@ -847,4 +849,6 @@ case FSApproximationEnum: return CreatePVectorFS(element); + case SSAHOApproximationEnum: + return CreatePVectorSSAHO(element); case NoneApproximationEnum: return NULL; @@ -3059,4 +3063,163 @@ /*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(3*numdofp); + IssmDouble* Bprime = xNew(3*numdofm); + IssmDouble D[3][3]={0.0}; // material matrix, simple scalar matrix. + IssmDouble D_scalar; + IssmDouble* Ke_gg = xNewZeroInit(numdofp*numdofm); + IssmDouble* Ke_gg_gaussian = xNew(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;iGetNode(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();igend();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; ivalues[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofp+j]; + for(i=0;ivalues[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(3*numnodes); + + /*Get dbasis in actual coordinate system: */ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); +}/*}}}*/ void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16906) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16907) @@ -71,4 +71,10 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16906) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16907) @@ -84,4 +84,5 @@ 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,4 +121,5 @@ 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; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16906) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16907) @@ -1326,4 +1326,9 @@ } /*}}}*/ +/*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,4 +3234,11 @@ ::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); } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16906) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16907) @@ -244,4 +244,5 @@ 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,4 +296,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16906) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16907) @@ -155,4 +155,5 @@ 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,4 +227,5 @@ 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");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16906) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16907) @@ -286,4 +286,5 @@ 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,4 +327,5 @@ 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);