Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16925) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16926) @@ -851,4 +851,6 @@ case SSAHOApproximationEnum: return CreatePVectorSSAHO(element); + case HOFSApproximationEnum: + return CreatePVectorHOFS(element); case NoneApproximationEnum: return NULL; @@ -3101,11 +3103,11 @@ int i,j; IssmDouble Jdet2d,alpha2; - IssmDouble *xyz_list_tria=NULL; - IssmDouble* L = xNewZeroInit(2*numdof); - IssmDouble DL[2][2] ={{ 0,0 },{0,0}}; //for basal drag + IssmDouble *xyz_list_tria = NULL; + IssmDouble* L = xNewZeroInit(2*numdof); + IssmDouble DL[2][2] = {{ 0,0 },{0,0}}; //for basal drag IssmDouble DL_scalar; - IssmDouble* Ke_gg = xNewZeroInit(numdof*numdof); - Node *node_list[2*numnodes]; - int* cs_list= xNew(2*numnodes); + IssmDouble* Ke_gg = xNewZeroInit(numdof*numdof); + Node **node_list = xNew(2*numnodes); + int* cs_list = xNew(2*numnodes); /*Initialize Element matrix and return if necessary*/ @@ -3163,4 +3165,5 @@ delete friction; xDelete(cs_list); + xDelete(node_list); xDelete(xyz_list_tria); xDelete(Ke_gg); @@ -3184,5 +3187,5 @@ IssmDouble D_scalar; IssmDouble* Ke_gg = xNewZeroInit(numdofp*numdofm); - Node *node_list[2*numnodes]; + Node **node_list = xNew(2*numnodes); int* cs_list= xNew(2*numnodes); @@ -3249,4 +3252,5 @@ xDelete(Ke_gg); xDelete(xyz_list); + xDelete(node_list); xDelete(cs_list); return Ke; @@ -3265,4 +3269,192 @@ return pe; }/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/ + + /*compute all load vectors for this element*/ + int init = element->FiniteElement(); + element->SetTemporaryElementType(P1Enum); + ElementVector* pe1=CreatePVectorHO(element); + element->SetTemporaryElementType(init); + ElementVector* pe2=CreatePVectorFS(element); + int indices[3]={18,19,20}; + element->SetTemporaryElementType(MINIcondensedEnum); + ElementMatrix* Ke = CreateKMatrixFS(element); + element->SetTemporaryElementType(init); + pe2->StaticCondensation(Ke,3,&indices[0]); + delete Ke; + ElementVector* pe3=CreatePVectorCouplingHOFS(element); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorCouplingHOFSViscous(element); + ElementVector* pe2=CreatePVectorCouplingHOFSFriction(element); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSFriction(Element* element){/*{{{*/ + + /*Intermediaries*/ + int i,j,approximation; + IssmDouble Jdet,Jdet2d,FSreconditioning; + IssmDouble bed_normal[3]; + IssmDouble viscosity, w, alpha2_gauss; + IssmDouble dw[3]; + IssmDouble *xyz_list_tria = NULL; + IssmDouble *xyz_list = NULL; + IssmDouble basis[6]; //for the six nodes of the penta + + /*Initialize Element vector and return if necessary*/ + if(!element->IsOnBed() || element->IsFloating()) return NULL; + element->GetInputValue(&approximation,ApproximationEnum); + if(approximation!=HOFSApproximationEnum) return NULL; + + int vnumnodes = element->GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + int numnodes = vnumnodes+pnumnodes; + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + Node **node_list = xNew(vnumnodes+pnumnodes); + for(i=0;iGetNode(i); + } + for(i=0;iGetNode(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* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_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();igend();ig++){ + + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss); + element->NodalFunctionsP1(basis, gauss); + + vzHO_input->GetInputValue(&w, gauss); + vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); + + element->NormalBase(&bed_normal[0],xyz_list_tria); + element->ViscosityFS(&viscosity,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); + + pe->Echo(); + element->CreatePVector(); + /*Clean up and return*/ + xDelete(cs_list); + xDelete(node_list); + xDelete(xyz_list); + xDelete(xyz_list_tria); + delete gauss; + delete friction; + return pe; +}/*}}}*/ +ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFSViscous(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!=HOFSApproximationEnum) return NULL; + int vnumnodes = element->GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + Node **node_list = xNew(vnumnodes+pnumnodes); + for(i=0;iGetNode(i); + } + for(i=0;iGetNode(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* vzHO_input=element->GetInput(VzHOEnum); _assert_(vzHO_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(5); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet, xyz_list,gauss); + element->NodalFunctionsP1(&basis[0],gauss); + element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); + + element->ViscosityFS(&viscosity,xyz_list,gauss,vx_input,vy_input,vz_input); + vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss); + + for(i=0;ivalues[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(cs_list); + xDelete(node_list); + xDelete(xyz_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. Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16925) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16926) @@ -76,4 +76,8 @@ ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element); ElementVector* CreatePVectorSSAHO(Element* element); + ElementVector* CreatePVectorHOFS(Element* element); + ElementVector* CreatePVectorCouplingHOFS(Element* element); + ElementVector* CreatePVectorCouplingHOFSFriction(Element* element); + ElementVector* CreatePVectorCouplingHOFSViscous(Element* element); void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 16925) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 16926) @@ -314,4 +314,7 @@ ::TransformSolutionCoord(values,this->nodes,numnodes,transformenum_list); }/*}}}*/ +void Element::TransformLoadVectorCoord(ElementVector* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/ + ::TransformLoadVectorCoord(Ke,nodes_list,numnodes,transformenum_list); +}/*}}}*/ void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16925) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16926) @@ -77,4 +77,5 @@ void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");}; /*Tiling only*/ void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");}; /*Tiling only*/ + void TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* transformenum_list); void TransformSolutionCoord(IssmDouble* values,int transformenum); void TransformSolutionCoord(IssmDouble* values,int* transformenum_list); @@ -106,7 +107,9 @@ virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0; virtual void NodalFunctions(IssmDouble* basis,Gauss* gauss)=0; + virtual void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss)=0; virtual void NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss)=0; virtual void NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss)=0; virtual void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0; + virtual void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0; virtual void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0; virtual void NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0; @@ -193,4 +196,5 @@ virtual void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0; virtual void ResetCoordinateSystem()=0; + virtual void SetTemporaryElementType(int element_type_in)=0; virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0; virtual void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16925) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16926) @@ -2506,4 +2506,12 @@ } /*}}}*/ +/*FUNCTION Penta::NodalFunctionsP1{{{*/ +void Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetNodalFunctionsP1(basis,(GaussPenta*)gauss); + +} +/*}}}*/ /*FUNCTION Penta::NodalFunctionsDerivatives{{{*/ void Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){ @@ -2511,4 +2519,12 @@ _assert_(gauss->Enum()==GaussPentaEnum); this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss); + +} +/*}}}*/ +/*FUNCTION Penta::NodalFunctionsP1Derivatives{{{*/ +void Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetNodalFunctionsP1Derivatives(dbasis,xyz_list,(GaussPenta*)gauss); } @@ -2847,4 +2863,9 @@ if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); else this->nodes=NULL; +} +/*}}}*/ +/*FUNCTION Penta::SetTemporaryElementType{{{*/ +void Penta::SetTemporaryElementType(int element_type_in){ + this->element_type=element_type_in; } /*}}}*/ @@ -8234,4 +8255,6 @@ ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list); + pe->Echo(); + _error_("stop"); /*Clean up and return*/ xDelete(cs_list); @@ -8396,4 +8419,5 @@ /*compute all load vectors for this element*/ int init = this->element_type; + printf("init init %i\n",init); this->element_type=P1Enum; ElementVector* pe1=CreatePVectorStressbalanceHO(); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16925) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16926) @@ -115,4 +115,5 @@ void ResetCoordinateSystem(void); void SmbGradients(); + void SetTemporaryElementType(int element_type_in); Element* SpawnBasalElement(void); IssmDouble SurfaceArea(void); @@ -257,5 +258,7 @@ ElementMatrix* NewElementMatrix(int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss); void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); + void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16925) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16926) @@ -108,7 +108,9 @@ IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); + void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; bool NoIceInElement(){_error_("not implemented yet");}; @@ -181,4 +183,5 @@ void ResetCoordinateSystem(void){_error_("not implemented yet");}; void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");}; + void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");}; void SmbGradients(){_error_("not implemented yet");}; IssmDouble SurfaceArea(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16925) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16926) @@ -288,9 +288,12 @@ ElementMatrix* NewElementMatrix(int approximation_enum); void NodalFunctions(IssmDouble* basis,Gauss* gauss); + void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); + void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss); void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss); void SetClone(int* minranks); + void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");}; Seg* SpawnSeg(int index1,int index2); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};