Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18919) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18920) @@ -306,7 +306,4 @@ } /*}}}*/ -void Element::DeleteMaterials(void){/*{{{*/ - delete this->material; -}/*}}}*/ void Element::DeepEcho(void){/*{{{*/ @@ -343,38 +340,13 @@ } /*}}}*/ -void Element::Echo(void){/*{{{*/ - _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n"); - _printf_(" id : "<id <<"\n"); - _printf_(" sid: "<sid<<"\n"); - if(vertices){ - int numvertices = this->GetNumberOfVertices(); - for(int i=0;iEcho(); - } - else _printf_("vertices = NULL\n"); - - if(nodes){ - int numnodes = this->GetNumberOfNodes(); - for(int i=0;iEcho(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->Echo(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->Echo(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->Echo(); - else _printf_("inputs=NULL\n"); -} -/*}}}*/ +void Element::DeleteInput(int input_enum){/*{{{*/ + + inputs->DeleteInput(input_enum); + +} +/*}}}*/ +void Element::DeleteMaterials(void){/*{{{*/ + delete this->material; +}/*}}}*/ IssmDouble Element::Divergence(void){/*{{{*/ /*Compute element divergence*/ @@ -419,4 +391,58 @@ return divergence; }/*}}}*/ +void Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble eps_eff; + IssmDouble eps0=1.e-27; + + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); + } + else{ + /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]); + } + /*Get viscosity*/ + material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ +void Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ + IssmDouble epsilon2d[2];/* epsilon=[exx,eyy,exy]; */ + IssmDouble eps_eff; + IssmDouble eps0=1.e-27; + + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); + } + else{ + /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ + this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1] + eps0*eps0); + } + /*Get viscosity*/ + material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ void Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ @@ -446,58 +472,4 @@ } /*}}}*/ -void Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dmudB; - IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ - IssmDouble epsilon2d[2];/* epsilon=[exx,eyy,exy]; */ - IssmDouble eps_eff; - IssmDouble eps0=1.e-27; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); - } - else{ - /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ - this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1] + eps0*eps0); - } - /*Get viscosity*/ - material->GetViscosity_B(&dmudB,eps_eff); - - /*Assign output pointer*/ - *pdmudB=dmudB; - -} -/*}}}*/ -void Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dmudB; - IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ - IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ - IssmDouble eps_eff; - IssmDouble eps0=1.e-27; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); - eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); - } - else{ - /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ - this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]); - } - /*Get viscosity*/ - material->GetViscosity_B(&dmudB,eps_eff); - - /*Assign output pointer*/ - *pdmudB=dmudB; - -} -/*}}}*/ void Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ @@ -527,15 +499,46 @@ } /*}}}*/ -void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ - matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); +void Element::Echo(void){/*{{{*/ + _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n"); + _printf_(" id : "<id <<"\n"); + _printf_(" sid: "<sid<<"\n"); + if(vertices){ + int numvertices = this->GetNumberOfVertices(); + for(int i=0;iEcho(); + } + else _printf_("vertices = NULL\n"); + + if(nodes){ + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iEcho(); + else _printf_("material = NULL\n"); + + if (matpar) matpar->Echo(); + else _printf_("matpar = NULL\n"); + + _printf_(" parameters\n"); + if (parameters) parameters->Echo(); + else _printf_("parameters = NULL\n"); + + _printf_(" inputs\n"); + if (inputs) inputs->Echo(); + else _printf_("inputs=NULL\n"); +} +/*}}}*/ +IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ + return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); +}/*}}}*/ +IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/ + return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); }/*}}}*/ void Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure); -}/*}}}*/ -IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ - return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); -}/*}}}*/ -IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/ - return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); }/*}}}*/ void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/ @@ -574,4 +577,28 @@ } /*}}}*/ +void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/ + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); + + /*First, figure out size of doflist and create it: */ + int numberofdofs=0; + for(int i=vnumnodes;iGetNumberOfDofs(FSApproximationEnum,setenum); + + /*Allocate output*/ + int* doflist=xNew(numberofdofs); + + /*Populate: */ + int count=0; + for(int i=vnumnodes;iGetDofList(doflist+count,FSApproximationEnum,setenum); + count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum); + } + + /*Assign output pointers:*/ + *pdoflist=doflist; +} +/*}}}*/ void Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/ @@ -597,28 +624,165 @@ } /*}}}*/ -void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/ - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - /*First, figure out size of doflist and create it: */ - int numberofdofs=0; - for(int i=vnumnodes;iGetNumberOfDofs(FSApproximationEnum,setenum); - - /*Allocate output*/ - int* doflist=xNew(numberofdofs); - - /*Populate: */ - int count=0; - for(int i=vnumnodes;iGetDofList(doflist+count,FSApproximationEnum,setenum); - count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ +Input* Element::GetInput(int inputenum){/*{{{*/ + return inputs->GetInput(inputenum); +}/*}}}*/ +void Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/ + + _assert_(pvalue); + + Input *input = this->GetInput(enumtype); + int numnodes = this->GetNumberOfNodes(); + + /* Start looping on the number of vertices: */ + if(input){ + Gauss* gauss=this->NewGauss(); + for(int iv=0;ivGaussNode(this->FiniteElement(),iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; + } + else{ + for(int iv=0;ivGetNumberOfNodes(); + Input *input = this->GetInput(enumtype); + if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + + /* Start looping on the number of vertices: */ + Gauss* gauss=this->NewGauss(); + for(int iv=0;ivGaussNode(this->FiniteElement(),iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; +} +/*}}}*/ +void Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/ + + _assert_(pvalue); + + int numnodes = this->NumberofNodesVelocity(); + Input *input = this->GetInput(enumtype); + if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + + /* Start looping on the number of vertices: */ + Gauss* gauss=this->NewGauss(); + for(int iv=0;ivGaussNode(this->VelocityInterpolation(),iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; +} +/*}}}*/ +void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/ + + /*Recover input*/ + Input* input=this->GetInput(enumtype); + if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + + /*Fetch number vertices for this element*/ + int numvertices = this->GetNumberOfVertices(); + + /*Checks in debugging mode*/ + _assert_(pvalue); + + /* Start looping on the number of vertices: */ + Gauss*gauss=this->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&pvalue[iv],gauss); + } + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/ + + /*Recover input*/ + Input* input=this->GetInput(enumtype); + + /*Checks in debugging mode*/ + _assert_(pvalue); + + /*Fetch number vertices for this element*/ + int numvertices = this->GetNumberOfVertices(); + + /* Start looping on the number of vertices: */ + if (input){ + Gauss* gauss=this->NewGauss(); + for (int iv=0;ivGaussVertex(iv); + input->GetInputValue(&pvalue[iv],gauss); + } + delete gauss; + } + else{ + for(int iv=0;ivGetNumberOfNodes(); + + /*Some checks to avoid segmentation faults*/ + _assert_(ug); + _assert_(numnodes>0); + _assert_(nodes); + + /*Get element minimum/maximum*/ + IssmDouble input_min = ug[nodes[0]->GetDof(0,GsetEnum)]; + IssmDouble input_max = input_min; + for(int i=1;iGetDof(0,GsetEnum)] < input_min) input_min = ug[nodes[i]->GetDof(0,GsetEnum)]; + if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)]; + } + + + /*Second loop to reassign min and max with local extrema*/ + for(int i=0;iGetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min; + if(max[nodes[i]->GetDof(0,GsetEnum)]GetDof(0,GsetEnum)] = input_max; + } +} +/*}}}*/ +void Element::GetInputValue(bool* pvalue,int inputenum){/*{{{*/ + + Input* input=inputs->GetInput(inputenum); + if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); + input->GetInputValue(pvalue); + +}/*}}}*/ +void Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/ + + Input* input=inputs->GetInput(inputenum); + if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); + input->GetInputValue(pvalue); + +}/*}}}*/ +void Element::GetInputValue(IssmDouble* pvalue,int inputenum){/*{{{*/ + + Input* input=inputs->GetInput(inputenum); + if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); + input->GetInputValue(pvalue); + +}/*}}}*/ +void Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/ + + Input* input=inputs->GetInput(inputenum); + if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); + input->GetInputValue(pvalue,gauss); + +}/*}}}*/ IssmDouble Element::GetMaterialParameter(int enum_in){/*{{{*/ @@ -635,4 +799,24 @@ } }/*}}}*/ +void Element::GetNodesLidList(int* lidlist){/*{{{*/ + + _assert_(lidlist); + _assert_(nodes); + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iLid(); + } +} +/*}}}*/ +void Element::GetNodesSidList(int* sidlist){/*{{{*/ + + _assert_(sidlist); + _assert_(nodes); + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iSid(); + } +} +/*}}}*/ void Element::GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity){/*{{{*/ /*Compute deformational heating from epsilon and viscosity */ @@ -672,185 +856,4 @@ } /*}}}*/ -Input* Element::GetInput(int inputenum){/*{{{*/ - return inputs->GetInput(inputenum); -}/*}}}*/ -void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/ - - /*Recover input*/ - Input* input=this->GetInput(enumtype); - if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); - - /*Fetch number vertices for this element*/ - int numvertices = this->GetNumberOfVertices(); - - /*Checks in debugging mode*/ - _assert_(pvalue); - - /* Start looping on the number of vertices: */ - Gauss*gauss=this->NewGauss(); - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&pvalue[iv],gauss); - } - - /*clean-up*/ - delete gauss; -} -/*}}}*/ -void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/ - - /*Recover input*/ - Input* input=this->GetInput(enumtype); - - /*Checks in debugging mode*/ - _assert_(pvalue); - - /*Fetch number vertices for this element*/ - int numvertices = this->GetNumberOfVertices(); - - /* Start looping on the number of vertices: */ - if (input){ - Gauss* gauss=this->NewGauss(); - for (int iv=0;ivGaussVertex(iv); - input->GetInputValue(&pvalue[iv],gauss); - } - delete gauss; - } - else{ - for(int iv=0;ivGetInput(enumtype); - int numnodes = this->GetNumberOfNodes(); - - /* Start looping on the number of vertices: */ - if(input){ - Gauss* gauss=this->NewGauss(); - for(int iv=0;ivGaussNode(this->FiniteElement(),iv); - input->GetInputValue(&pvalue[iv],gauss); - } - delete gauss; - } - else{ - for(int iv=0;ivGetNumberOfNodes(); - Input *input = this->GetInput(enumtype); - if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); - - /* Start looping on the number of vertices: */ - Gauss* gauss=this->NewGauss(); - for(int iv=0;ivGaussNode(this->FiniteElement(),iv); - input->GetInputValue(&pvalue[iv],gauss); - } - delete gauss; -} -/*}}}*/ -void Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/ - - _assert_(pvalue); - - int numnodes = this->NumberofNodesVelocity(); - Input *input = this->GetInput(enumtype); - if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); - - /* Start looping on the number of vertices: */ - Gauss* gauss=this->NewGauss(); - for(int iv=0;ivGaussNode(this->VelocityInterpolation(),iv); - input->GetInputValue(&pvalue[iv],gauss); - } - delete gauss; -} -/*}}}*/ -void Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/ - - - /*Get number of nodes for this element*/ - int numnodes = this->GetNumberOfNodes(); - - /*Some checks to avoid segmentation faults*/ - _assert_(ug); - _assert_(numnodes>0); - _assert_(nodes); - - /*Get element minimum/maximum*/ - IssmDouble input_min = ug[nodes[0]->GetDof(0,GsetEnum)]; - IssmDouble input_max = input_min; - for(int i=1;iGetDof(0,GsetEnum)] < input_min) input_min = ug[nodes[i]->GetDof(0,GsetEnum)]; - if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)]; - } - - - /*Second loop to reassign min and max with local extrema*/ - for(int i=0;iGetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min; - if(max[nodes[i]->GetDof(0,GsetEnum)]GetDof(0,GsetEnum)] = input_max; - } -} -/*}}}*/ -void Element::GetInputValue(bool* pvalue,int inputenum){/*{{{*/ - - Input* input=inputs->GetInput(inputenum); - if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); - input->GetInputValue(pvalue); - -}/*}}}*/ -void Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/ - - Input* input=inputs->GetInput(inputenum); - if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); - input->GetInputValue(pvalue); - -}/*}}}*/ -void Element::GetInputValue(IssmDouble* pvalue,int inputenum){/*{{{*/ - - Input* input=inputs->GetInput(inputenum); - if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); - input->GetInputValue(pvalue); - -}/*}}}*/ -void Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/ - - Input* input=inputs->GetInput(inputenum); - if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); - input->GetInputValue(pvalue,gauss); - -}/*}}}*/ -void Element::GetNodesSidList(int* sidlist){/*{{{*/ - - _assert_(sidlist); - _assert_(nodes); - int numnodes = this->GetNumberOfNodes(); - for(int i=0;iSid(); - } -} -/*}}}*/ -void Element::GetNodesLidList(int* lidlist){/*{{{*/ - - _assert_(lidlist); - _assert_(nodes); - int numnodes = this->GetNumberOfNodes(); - for(int i=0;iLid(); - } -} -/*}}}*/ void Element::GetVectorFromInputs(Vector* vector,int input_enum){/*{{{*/ @@ -878,4 +881,10 @@ } /*}}}*/ +void Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/ + + int numvertices = this->GetNumberOfVertices(); + for(int i=0;ivertices[i]->Connectivity(); +} +/*}}}*/ void Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/ @@ -891,10 +900,4 @@ int numvertices = this->GetNumberOfVertices(); for(int i=0;ivertices[i]->Sid(); -} -/*}}}*/ -void Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/ - - int numvertices = this->GetNumberOfVertices(); - for(int i=0;ivertices[i]->Connectivity(); } /*}}}*/ @@ -1055,10 +1058,4 @@ /*Call inputs method*/ if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); - -} -/*}}}*/ -void Element::DeleteInput(int input_enum){/*{{{*/ - - inputs->DeleteInput(input_enum); } @@ -1298,14 +1295,14 @@ } /*}}}*/ +ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/ + return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); +} +/*}}}*/ +ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/ + return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum); +} +/*}}}*/ ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/ return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/ - return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); -} -/*}}}*/ -ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/ - return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum); } /*}}}*/ @@ -1393,4 +1390,12 @@ *pnodesperelement = input->GetResultNumberOfNodes(); }/*}}}*/ +void Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/ + + Input* input=this->inputs->GetInput(output_enum); + if(!input) _error_("input "<ResultToPatch(values,nodesperelement,this->Sid()); + +} /*}}}*/ void Element::ResultToVector(Vector* vector,int output_enum){/*{{{*/ @@ -1440,12 +1445,4 @@ _error_("interpolation "<GetResultInterpolation())<<" not supported yet"); } -} /*}}}*/ -void Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/ - - Input* input=this->inputs->GetInput(output_enum); - if(!input) _error_("input "<ResultToPatch(values,nodesperelement,this->Sid()); - } /*}}}*/ void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/ @@ -1525,7 +1522,119 @@ } /*}}}*/ -IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/ - _assert_(matpar); - return this->matpar->TMeltingPoint(pressure); +void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + /*Compute the 3d Strain Rate (6 components): + * + * epsilon=[exx eyy ezz exy exz eyz] + */ + + /*Intermediaries*/ + IssmDouble dvx[3]; + IssmDouble dvy[3]; + IssmDouble dvz[3]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !vz_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + epsilon[0] = dvx[0]; + epsilon[1] = dvy[1]; + epsilon[2] = dvz[2]; + epsilon[3] = 0.5*(dvx[1] + dvy[0]); + epsilon[4] = 0.5*(dvx[2] + dvz[0]); + epsilon[5] = 0.5*(dvy[2] + dvz[1]); + +}/*}}}*/ +void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Compute the 3d Blatter/HOStrain Rate (5 components): + * + * epsilon=[exx eyy exy exz eyz] + * + * with exz=1/2 du/dz + * eyz=1/2 dv/dz + * + * the contribution of vz is neglected + */ + + /*Intermediaries*/ + IssmDouble dvx[3]; + IssmDouble dvy[3]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + epsilon[0] = dvx[0]; + epsilon[1] = dvy[1]; + epsilon[2] = 0.5*(dvx[1] + dvy[0]); + epsilon[3] = 0.5*dvx[2]; + epsilon[4] = 0.5*dvy[2]; + +}/*}}}*/ +void Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Compute the 2d Blatter/HOStrain Rate (2 components): + * + * epsilon=[exx exz] + * + * with exz=1/2 du/dz + * + * the contribution of vz is neglected + */ + + /*Intermediaries*/ + IssmDouble dvx[3]; + + /*Check that both inputs have been found*/ + if (!vx_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + epsilon[0] = dvx[0]; + epsilon[1] = 0.5*dvx[1]; + +}/*}}}*/ +void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dvx[3]; + IssmDouble dvy[3]; + + /*Check that both inputs have been found*/ + if(!vx_input || !vy_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + epsilon[0] = dvx[0]; + epsilon[1] = dvy[1]; + epsilon[2] = 0.5*(dvx[1] + dvy[0]); + +}/*}}}*/ +void Element::StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dvx[3]; + + /*Check that both inputs have been found*/ + if (!vx_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + *epsilon = dvx[0]; + }/*}}}*/ void Element::StressMaxPrincipalCreateInput(void){/*{{{*/ @@ -1615,46 +1724,234 @@ } /*}}}*/ -void Element::ViscousHeatingCreateInput(void){/*{{{*/ - - /*Intermediaries*/ - IssmDouble phi; - IssmDouble viscosity; - IssmDouble epsilon[6]; - IssmDouble thickness; - IssmDouble *xyz_list = NULL; - - /*Fetch number vertices and allocate memory*/ - int numvertices = this->GetNumberOfVertices(); - IssmDouble* viscousheating = xNew(numvertices); - - /*Retrieve all inputs and parameters*/ - this->GetVerticesCoordinatesBase(&xyz_list); - Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = this->GetInput(VzEnum); _assert_(vz_input); - Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input); - - /*loop over vertices: */ - Gauss* gauss=this->NewGauss(); - for (int iv=0;ivGaussVertex(iv); - - thickness_input->GetInputValue(&thickness,gauss); - - this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); - this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); - this->GetPhi(&phi,&epsilon[0],viscosity); - - viscousheating[iv]=phi*thickness; - } - - /*Create PentaVertex input, which will hold the basal friction:*/ - this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum); - - /*Clean up and return*/ - xDelete(viscousheating); - delete gauss; -} -/*}}}*/ +void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ + matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); +}/*}}}*/ +IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/ + _assert_(matpar); + return this->matpar->TMeltingPoint(pressure); +}/*}}}*/ +void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ + + /*All nodes have the same Coordinate System*/ + int numnodes = this->GetNumberOfNodes(); + int* cs_array = xNew(numnodes); + for(int i=0;inodes,numnodes,cs_array); + + /*Clean-up*/ + xDelete(cs_array); +}/*}}}*/ +void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ + + int i,j; + int numdofs = 0; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; + + /*Get total number of dofs*/ + for(i=0;i(Ke->nrows*Ke->ncols); + for(i=0;inrows;i++) for(j=0;jncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j]; + + /*Get Coordinate Systems transform matrix*/ + CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array); + + /*Transform matrix: R*Ke*R^T */ + TripleMultiply(transform,numdofs,numdofs,0, + values,Ke->nrows,Ke->ncols,0, + transform,numdofs,numdofs,1, + &Ke->values[0],0); + + /*Free Matrix*/ + xDelete(transform); + xDelete(values); +}/*}}}*/ +void Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/ + + /*All nodes have the same Coordinate System*/ + int numnodes = this->GetNumberOfNodes(); + int* cs_array = xNew(numnodes); + for(int i=0;iTransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array); + + /*Clean-up*/ + xDelete(cs_array); +}/*}}}*/ +void Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/ + + this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array); + +}/*}}}*/ +void Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ + + int i; + int numdofs = 0; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; + + /*Get total number of dofs*/ + for(i=0;i(pe->nrows); + for(i=0;inrows;i++) values[i]=pe->values[i]; + + /*Get Coordinate Systems transform matrix*/ + CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array); + + /*Transform matrix: R^T*pe */ + MatrixMultiply(transform,numdofs,numdofs,1, + values,pe->nrows,1,0, + &pe->values[0],0); + + /*Free Matrices*/ + xDelete(transform); + xDelete(values); +}/*}}}*/ +void Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/ + + /*All nodes have the same Coordinate System*/ + int numnodes = this->GetNumberOfNodes(); + int* cs_array = xNew(numnodes); + for(int i=0;iTransformSolutionCoord(values,this->nodes,numnodes,cs_array); + + /*Clean-up*/ + xDelete(cs_array); +}/*}}}*/ +void Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/ + this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list); +}/*}}}*/ +void Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/ + + /*All nodes have the same Coordinate System*/ + int* cs_array = xNew(numnodes); + for(int i=0;iTransformSolutionCoord(values,this->nodes,numnodes,cs_array); + + /*Clean-up*/ + xDelete(cs_array); +}/*}}}*/ +void Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/ + this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array); +}/*}}}*/ +void Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/ + /*NOT NEEDED*/ + /*All nodes have the same Coordinate System*/ + int* cs_array = xNew(numnodes); + for(int i=0;iTransformSolutionCoord(values,nodes_list,numnodes,cs_array); + + /*Clean-up*/ + xDelete(cs_array); +}/*}}}*/ +void Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ + + int i; + int numdofs = 0; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; + + /*Get total number of dofs*/ + for(i=0;i(numdofs); + for(i=0;i(transform); + xDelete(values); +}/*}}}*/ +void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ + + /*All nodes have the same Coordinate System*/ + int numnodes = this->GetNumberOfNodes(); + int* cs_array = xNew(numnodes); + for(int i=0;iTransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array); + + /*Clean-up*/ + xDelete(cs_array); +}/*}}}*/ +void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/ + this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list); +}/*}}}*/ +void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ + + int numdofs = 0; + IssmDouble *transform = NULL; + IssmDouble *values = NULL; + + /*Get total number of dofs*/ + for(int i=0;i(Ke->nrows*Ke->ncols); + for(int i=0;inrows*Ke->ncols;i++) values[i]=Ke->values[i]; + + /*Get Coordinate Systems transform matrix*/ + CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array); + + /*Transform matrix: R^T*Ke*R */ + TripleMultiply(transform,numdofs,numdofs,1, + values,Ke->nrows,Ke->ncols,0, + transform,numdofs,numdofs,0, + &Ke->values[0],0); + + /*Free Matrix*/ + xDelete(transform); + xDelete(values); +}/*}}}*/ void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ /*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9, @@ -1696,4 +1993,35 @@ } /*}}}*/ +void Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ +void Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ + IssmDouble eps_eff; + + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]); + } + else{ + /* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/ + this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]); + } + + /*Get viscosity*/ + material->GetViscosity(&viscosity,eps_eff); + + /*Assign output pointer*/ + *pviscosity=viscosity; +}/*}}}*/ +void Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); +}/*}}}*/ void Element::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ /*Compute the L1L2 viscosity @@ -1755,29 +2083,4 @@ *pviscosity = viscosity; }/*}}}*/ -void Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble viscosity; - IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ - IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ - IssmDouble eps_eff; - - if(dim==3){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]); - } - else{ - /* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/ - this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]); - } - - /*Get viscosity*/ - material->GetViscosity(&viscosity,eps_eff); - - /*Assign output pointer*/ - *pviscosity=viscosity; -}/*}}}*/ void Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ @@ -1808,347 +2111,44 @@ this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon); }/*}}}*/ -void Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -}/*}}}*/ -void Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ - this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon); -}/*}}}*/ -void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ - /*Compute the 3d Strain Rate (6 components): - * - * epsilon=[exx eyy ezz exy exz eyz] - */ +void Element::ViscousHeatingCreateInput(void){/*{{{*/ /*Intermediaries*/ - IssmDouble dvx[3]; - IssmDouble dvy[3]; - IssmDouble dvz[3]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input || !vz_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); - epsilon[0] = dvx[0]; - epsilon[1] = dvy[1]; - epsilon[2] = dvz[2]; - epsilon[3] = 0.5*(dvx[1] + dvy[0]); - epsilon[4] = 0.5*(dvx[2] + dvz[0]); - epsilon[5] = 0.5*(dvy[2] + dvz[1]); - -}/*}}}*/ -void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - /*Compute the 3d Blatter/HOStrain Rate (5 components): - * - * epsilon=[exx eyy exy exz eyz] - * - * with exz=1/2 du/dz - * eyz=1/2 dv/dz - * - * the contribution of vz is neglected - */ - - /*Intermediaries*/ - IssmDouble dvx[3]; - IssmDouble dvy[3]; - - /*Check that both inputs have been found*/ - if (!vx_input || !vy_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - epsilon[0] = dvx[0]; - epsilon[1] = dvy[1]; - epsilon[2] = 0.5*(dvx[1] + dvy[0]); - epsilon[3] = 0.5*dvx[2]; - epsilon[4] = 0.5*dvy[2]; - -}/*}}}*/ -void Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - /*Compute the 2d Blatter/HOStrain Rate (2 components): - * - * epsilon=[exx exz] - * - * with exz=1/2 du/dz - * - * the contribution of vz is neglected - */ - - /*Intermediaries*/ - IssmDouble dvx[3]; - - /*Check that both inputs have been found*/ - if (!vx_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - epsilon[0] = dvx[0]; - epsilon[1] = 0.5*dvx[1]; - -}/*}}}*/ -void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dvx[3]; - IssmDouble dvy[3]; - - /*Check that both inputs have been found*/ - if(!vx_input || !vy_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); - epsilon[0] = dvx[0]; - epsilon[1] = dvy[1]; - epsilon[2] = 0.5*(dvx[1] + dvy[0]); - -}/*}}}*/ -void Element::StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input){/*{{{*/ - - /*Intermediaries*/ - IssmDouble dvx[3]; - - /*Check that both inputs have been found*/ - if (!vx_input){ - _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << "\n"); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); - *epsilon = dvx[0]; - -}/*}}}*/ -void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ - - /*All nodes have the same Coordinate System*/ - int numnodes = this->GetNumberOfNodes(); - int* cs_array = xNew(numnodes); - for(int i=0;inodes,numnodes,cs_array); - - /*Clean-up*/ - xDelete(cs_array); -}/*}}}*/ -void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ - - int i,j; - int numdofs = 0; - IssmDouble *transform = NULL; - IssmDouble *values = NULL; - - /*Get total number of dofs*/ - for(i=0;i(Ke->nrows*Ke->ncols); - for(i=0;inrows;i++) for(j=0;jncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j]; - - /*Get Coordinate Systems transform matrix*/ - CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array); - - /*Transform matrix: R*Ke*R^T */ - TripleMultiply(transform,numdofs,numdofs,0, - values,Ke->nrows,Ke->ncols,0, - transform,numdofs,numdofs,1, - &Ke->values[0],0); - - /*Free Matrix*/ - xDelete(transform); - xDelete(values); -}/*}}}*/ -void Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/ - - /*All nodes have the same Coordinate System*/ - int numnodes = this->GetNumberOfNodes(); - int* cs_array = xNew(numnodes); - for(int i=0;iTransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array); - - /*Clean-up*/ - xDelete(cs_array); -}/*}}}*/ -void Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/ - - this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array); - -}/*}}}*/ -void Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ - - int i; - int numdofs = 0; - IssmDouble *transform = NULL; - IssmDouble *values = NULL; - - /*Get total number of dofs*/ - for(i=0;i(pe->nrows); - for(i=0;inrows;i++) values[i]=pe->values[i]; - - /*Get Coordinate Systems transform matrix*/ - CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array); - - /*Transform matrix: R^T*pe */ - MatrixMultiply(transform,numdofs,numdofs,1, - values,pe->nrows,1,0, - &pe->values[0],0); - - /*Free Matrices*/ - xDelete(transform); - xDelete(values); -}/*}}}*/ -void Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/ - - /*All nodes have the same Coordinate System*/ - int numnodes = this->GetNumberOfNodes(); - int* cs_array = xNew(numnodes); - for(int i=0;iTransformSolutionCoord(values,this->nodes,numnodes,cs_array); - - /*Clean-up*/ - xDelete(cs_array); -}/*}}}*/ -void Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/ - this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list); -}/*}}}*/ -void Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/ - - /*All nodes have the same Coordinate System*/ - int* cs_array = xNew(numnodes); - for(int i=0;iTransformSolutionCoord(values,this->nodes,numnodes,cs_array); - - /*Clean-up*/ - xDelete(cs_array); -}/*}}}*/ -void Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/ - this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array); -}/*}}}*/ -void Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/ - /*NOT NEEDED*/ - /*All nodes have the same Coordinate System*/ - int* cs_array = xNew(numnodes); - for(int i=0;iTransformSolutionCoord(values,nodes_list,numnodes,cs_array); - - /*Clean-up*/ - xDelete(cs_array); -}/*}}}*/ -void Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ - - int i; - int numdofs = 0; - IssmDouble *transform = NULL; - IssmDouble *values = NULL; - - /*Get total number of dofs*/ - for(i=0;i(numdofs); - for(i=0;i(transform); - xDelete(values); -}/*}}}*/ -void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ - - /*All nodes have the same Coordinate System*/ - int numnodes = this->GetNumberOfNodes(); - int* cs_array = xNew(numnodes); - for(int i=0;iTransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array); - - /*Clean-up*/ - xDelete(cs_array); -}/*}}}*/ -void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/ - this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list); -}/*}}}*/ -void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/ - - int numdofs = 0; - IssmDouble *transform = NULL; - IssmDouble *values = NULL; - - /*Get total number of dofs*/ - for(int i=0;i(Ke->nrows*Ke->ncols); - for(int i=0;inrows*Ke->ncols;i++) values[i]=Ke->values[i]; - - /*Get Coordinate Systems transform matrix*/ - CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array); - - /*Transform matrix: R^T*Ke*R */ - TripleMultiply(transform,numdofs,numdofs,1, - values,Ke->nrows,Ke->ncols,0, - transform,numdofs,numdofs,0, - &Ke->values[0],0); - - /*Free Matrix*/ - xDelete(transform); - xDelete(values); -}/*}}}*/ + IssmDouble phi; + IssmDouble viscosity; + IssmDouble epsilon[6]; + IssmDouble thickness; + IssmDouble *xyz_list = NULL; + + /*Fetch number vertices and allocate memory*/ + int numvertices = this->GetNumberOfVertices(); + IssmDouble* viscousheating = xNew(numvertices); + + /*Retrieve all inputs and parameters*/ + this->GetVerticesCoordinatesBase(&xyz_list); + Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = this->GetInput(VzEnum); _assert_(vz_input); + Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input); + + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + for (int iv=0;ivGaussVertex(iv); + + thickness_input->GetInputValue(&thickness,gauss); + + this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input); + this->GetPhi(&phi,&epsilon[0],viscosity); + + viscousheating[iv]=phi*thickness; + } + + /*Create PentaVertex input, which will hold the basal friction:*/ + this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum); + + /*Clean up and return*/ + xDelete(viscousheating); + delete gauss; +} +/*}}}*/