Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17296) @@ -142,4 +142,74 @@ this->parameters->FindParam(pvalues,psize,paramenum); }/*}}}*/ +void Element::GetDofList(int** pdoflist,int approximation_enum,int setenum){/*{{{*/ + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->GetNumberOfNodes(); + + /*First, figure out size of doflist and create it: */ + int numberofdofs=0; + for(int i=0;iGetNumberOfDofs(approximation_enum,setenum); + + /*Allocate output*/ + int* doflist=xNew(numberofdofs); + + /*Populate: */ + int count=0; + for(int i=0;iGetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + /*Assign output pointers:*/ + *pdoflist=doflist; +} +/*}}}*/ +void Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/ + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->GetNumberOfNodesVelocity(); + + /*First, figure out size of doflist and create it: */ + int numberofdofs=0; + for(int i=0;iGetNumberOfDofs(FSvelocityEnum,setenum); + + /*Allocate output*/ + int* doflist=xNew(numberofdofs); + + /*Populate: */ + int count=0; + for(int i=0;iGetDofList(doflist+count,FSvelocityEnum,setenum); + count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum); + } + + /*Assign output pointers:*/ + *pdoflist=doflist; +} +/*}}}*/ +void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/ + + /*Fetch number of nodes and dof for this finite element*/ + int vnumnodes = this->GetNumberOfNodesVelocity(); + int pnumnodes = this->GetNumberOfNodesPressure(); + + /*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; +} +/*}}}*/ IssmDouble Element::GetMaterialParameter(int enum_in){/*{{{*/ @@ -556,15 +626,22 @@ *pviscosity=viscosity; }/*}}}*/ -void Element::ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ +void Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ /*Intermediaries*/ IssmDouble viscosity; - IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ IssmDouble eps_eff; - /*Get effective strain rate - * eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/ - this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + epsilon[0]*epsilon[1]); + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = exx^2*/ + this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d); + } /*Get viscosity*/ @@ -684,4 +761,19 @@ }/*}}}*/ +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){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17296) @@ -60,4 +60,7 @@ void FindParam(IssmDouble* pvalue,int paramenum); void FindParam(int** pvalues,int* psize,int paramenum); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); + void GetDofListVelocity(int** pdoflist,int setenum); + void GetDofListPressure(int** pdoflist,int setenum); Input* GetInput(int inputenum); IssmDouble GetMaterialParameter(int enum_in); @@ -71,4 +74,5 @@ void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum); void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input); void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); @@ -96,5 +100,5 @@ void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); - void ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); void ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); @@ -135,7 +139,5 @@ virtual Element* GetSurfaceElement(void)=0; virtual Element* GetBasalElement(void)=0; - virtual void GetDofList(int** pdoflist,int approximation_enum,int setenum)=0; - virtual void GetDofListVelocity(int** pdoflist,int setenum)=0; - virtual void GetDofListPressure(int** pdoflist,int setenum)=0; + virtual void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; virtual void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17296) @@ -631,77 +631,4 @@ Element* element=this->GetBasalPenta(); return element; -} -/*}}}*/ -/*FUNCTION Penta::GetDofList {{{*/ -void Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){ - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*First, figure out size of doflist and create it: */ - int numberofdofs=0; - for(int i=0;iGetNumberOfDofs(approximation_enum,setenum); - - /*Allocate output*/ - int* doflist=xNew(numberofdofs); - - /*Populate: */ - int count=0; - for(int i=0;iGetDofList(doflist+count,approximation_enum,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ -/*FUNCTION Penta::GetDofListVelocity{{{*/ -void Penta::GetDofListVelocity(int** pdoflist,int setenum){ - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodesVelocity(); - - /*First, figure out size of doflist and create it: */ - int numberofdofs=0; - for(int i=0;iGetNumberOfDofs(FSvelocityEnum,setenum); - - /*Allocate output*/ - int* doflist=xNew(numberofdofs); - - /*Populate: */ - int count=0; - for(int i=0;iGetDofList(doflist+count,FSvelocityEnum,setenum); - count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ -/*FUNCTION Penta::GetDofListPressure{{{*/ -void Penta::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; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17296) @@ -82,7 +82,4 @@ Element* GetSurfaceElement(void); Element* GetBasalElement(void); - void GetDofList(int** pdoflist,int approximation_enum,int setenum); - void GetDofListVelocity(int** pdoflist,int setenum); - void GetDofListPressure(int** pdoflist,int setenum); void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17296) @@ -79,7 +79,4 @@ Element* GetSurfaceElement(void){_error_("not implemented yet");}; Element* GetBasalElement(void){_error_("not implemented yet");}; - void GetDofList(int** pdoflist,int approximation_enum,int setenum){_error_("not implemented yet");}; - void GetDofListVelocity(int** pdoflist,int setenum){_error_("not implemented yet");}; - void GetDofListPressure(int** pdoflist,int setenum){_error_("not implemented yet");}; void GetInputListOnNodes(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");}; void GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17296) @@ -223,5 +223,5 @@ /*Compute strain rate viscosity and pressure: */ this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - this->ViscositySSA(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input); + this->ViscositySSA(&viscosity,2,&xyz_list[0][0],gauss,vx_input,vy_input); pressure_input->GetInputValue(&pressure,gauss); @@ -488,77 +488,4 @@ } } -} -/*}}}*/ -/*FUNCTION Tria::GetDofList {{{*/ -void Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){ - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodes(); - - /*First, figure out size of doflist and create it: */ - int numberofdofs=0; - for(int i=0;iGetNumberOfDofs(approximation_enum,setenum); - - /*Allocate output*/ - int* doflist=xNew(numberofdofs); - - /*Populate: */ - int count=0; - for(int i=0;iGetDofList(doflist+count,approximation_enum,setenum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ -/*FUNCTION Tria::GetDofListVelocity{{{*/ -void Tria::GetDofListVelocity(int** pdoflist,int setenum){ - - /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->NumberofNodesVelocity(); - - /*First, figure out size of doflist and create it: */ - int numberofdofs=0; - for(int i=0;iGetNumberOfDofs(FSvelocityEnum,setenum); - - /*Allocate output*/ - int* doflist=xNew(numberofdofs); - - /*Populate: */ - int count=0; - for(int i=0;iGetDofList(doflist+count,FSvelocityEnum,setenum); - count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum); - } - - /*Assign output pointers:*/ - *pdoflist=doflist; -} -/*}}}*/ -/*FUNCTION Tria::GetDofListPressure{{{*/ -void Tria::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; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17295) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17296) @@ -76,7 +76,4 @@ Element* GetSurfaceElement(void){_error_("not implemented yet");}; Element* GetBasalElement(void){_error_("not implemented yet");}; - void GetDofList(int** pdoflist,int approximation_enum,int setenum); - void GetDofListVelocity(int** pdoflist,int setenum); - void GetDofListPressure(int** pdoflist,int setenum); void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating); IssmDouble GetGroundedPortion(IssmDouble* xyz_list);