Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 23961) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 23962) @@ -328,4 +328,32 @@ /*Clean up*/ xDelete(triabasis); +} +/*}}}*/ +void TriaRef::GetSegmentNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list_tria,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/ + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + _assert_(index1>=0 && index1<3); + _assert_(index2>=0 && index2<3); + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(finiteelement); + + /*Get nodal functions*/ + IssmDouble* dtriabasis=xNew(2*numnodes); + GetNodalFunctionsDerivatives(dtriabasis,xyz_list_tria,gauss,finiteelement); + + switch(finiteelement){ + case P1Enum: case P1DGEnum: + dbasis[2*0+0] = dtriabasis[numnodes*0+index1]; + dbasis[2*0+1] = dtriabasis[numnodes*1+index1]; + dbasis[2*1+0] = dtriabasis[numnodes*0+index2]; + dbasis[2*1+1] = dtriabasis[numnodes*1+index2]; + break; + default: + _error_("Element type "<(dtriabasis); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 23961) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 23962) @@ -29,4 +29,5 @@ void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void GetSegmentNodalFunctions(IssmDouble* basis,Gauss* gauss, int index1,int index2,int finiteelement); + void GetSegmentNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list_tria,Gauss* gauss, int index1,int index2,int finiteelement); void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*do nothing */}; void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement); Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23961) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 23962) @@ -14,5 +14,6 @@ #include "../classes.h" /*}}}*/ -#define NUMNODES 2 +#define NUMNODES 2 +#define NUMVERTICES 2 /*Channel constructors and destructor*/ @@ -33,4 +34,7 @@ this->nodes = NULL; + /*Set channel cross section to 0*/ + this->S = 0.; + /*Get edge info*/ int i1 = iomodel->faces[4*index+0]; @@ -71,4 +75,5 @@ /*copy fields: */ channel->id=this->id; + channel->S=this->S; /*point parameters: */ @@ -92,4 +97,5 @@ _printf_("Channel:\n"); _printf_(" id: " << id << "\n"); + _printf_(" S: " << S << "\n"); hnodes->DeepEcho(); hvertices->DeepEcho(); @@ -105,4 +111,5 @@ _printf_("Channel:\n"); _printf_(" id: " << id << "\n"); + _printf_(" S: " << S << "\n"); hnodes->Echo(); hvertices->Echo(); @@ -122,4 +129,5 @@ MARSHALLING_ENUM(ChannelEnum); MARSHALLING(id); + MARSHALLING(S); if(marshall_direction==MARSHALLING_BACKWARD){ @@ -141,9 +149,6 @@ /*}}}*/ int Channel::ObjectEnum(void){/*{{{*/ - return ChannelEnum; - -} -/*}}}*/ +}/*}}}*/ /*Load virtual functions definitions:*/ @@ -323,28 +328,155 @@ /*Channel specific functions*/ -ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/ - - _error_("not implemented :( "); +ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/ /*Initialize Element matrix and return if necessary*/ Tria* tria=(Tria*)element; if(!tria->IsIceInElement()) return NULL; - ElementVector* pe=new ElementVector(nodes,NUMNODES,this->parameters); + _assert_(tria->FiniteElement()==P1Enum); + int index1=tria->GetNodeIndex(nodes[0]); + int index2=tria->GetNodeIndex(nodes[1]); + + /*Intermediaries */ + IssmDouble Jdet,v1; + IssmDouble A,B,n,phi_old,phi,phi_0; + IssmDouble H,h,b; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[3][3]; + const int numnodes = NUMNODES; + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters); + IssmDouble basis[NUMNODES]; + IssmDouble dbasisdx[2*NUMNODES]; + IssmDouble dbasisds[NUMNODES]; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0] ,this->vertices,NUMVERTICES); + GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); + IssmDouble L = element->FindParam(MaterialsLatentheatEnum); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum); + IssmDouble g = element->FindParam(ConstantsGEnum); + Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); + Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + + /*Get tangent vector*/ + IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0]; + IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1]; + tx = tx/sqrt(tx*tx+ty*ty); + ty = ty/sqrt(tx*tx+ty*ty); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=new GaussTria(index1,index2,2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement()); + + /*Get input values at gauss points*/ + h_input->GetInputValue(&h,gauss); + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + phi_input->GetInputValue(&phi,gauss); + b_input->GetInputValue(&b,gauss); + H_input->GetInputValue(&H,gauss); + + /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/ + phi_0 = rho_water*g*b + rho_ice*g*H; + A=pow(B,-n); + v1 = 2./pow(n,n)*A*S*(pow(fabs(phi_0 - phi),n-1.)*( - n)); + for(int i=0;ivalues[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j]; + } + } + + /*Transient term if dt>0*/ + if(dt>0.){ + } + } /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +ElementVector* Channel::CreatePVectorHydrologyGlaDS(void){/*{{{*/ + + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + _assert_(tria->FiniteElement()==P1Enum); + int index1=tria->GetNodeIndex(nodes[0]); + int index2=tria->GetNodeIndex(nodes[1]); + + /*Intermediaries */ + IssmDouble Jdet,v2; + IssmDouble A,B,n,phi_old,phi,phi_0; + IssmDouble H,h,b; + IssmDouble xyz_list[NUMVERTICES][3]; + const int numnodes = NUMNODES; + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = new ElementVector(this->nodes,NUMNODES,this->parameters); + IssmDouble basis[NUMNODES]; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + IssmDouble L = element->FindParam(MaterialsLatentheatEnum); + IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum); + IssmDouble g = element->FindParam(ConstantsGEnum); + Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* b_input = element->GetInput(BedEnum); _assert_(b_input); + Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum); _assert_(phiold_input); + Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=new GaussTria(index1,index2,2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); + + /*Get input values at gauss points*/ + h_input->GetInputValue(&h,gauss); + B_input->GetInputValue(&B,gauss); + n_input->GetInputValue(&n,gauss); + phi_input->GetInputValue(&phi,gauss); + b_input->GetInputValue(&b,gauss); + H_input->GetInputValue(&H,gauss); + + /*Compute closing rate*/ + /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/ + phi_0 = rho_water*g*b + rho_ice*g*H; + A=pow(B,-n); + v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi)); + + for(int i=0;ivalues[i]+= - Jdet*gauss->weight*(-v2)*basis[i]; + + /*Transient term if dt>0*/ + if(dt>0.){ + //phiold_input->GetInputValue(&phi_old,gauss); + //for(int i=0;ivalues[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i]; + } + } + + /*Clean up and return*/ + delete gauss; return pe; } /*}}}*/ -ElementMatrix* Channel::CreateKMatrixHydrologyGlaDS(void){/*{{{*/ - - _error_("not implemented :( "); - - /*Initialize Element matrix and return if necessary*/ - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES,this->parameters); - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.h (revision 23961) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.h (revision 23962) @@ -17,4 +17,7 @@ class Channel: public Load { + + private: + IssmDouble S; public: