Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16848) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16849) @@ -139,5 +139,70 @@ }/*}}}*/ ElementVector* StressbalanceSIAAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + + int meshtype; + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + return CreatePVector2D(element); + case Mesh3DEnum: + return CreatePVector3D(element); + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize Element vector*/ + ElementVector* pe=element->NewElementVector(); + + /*Retrieve all inputs and parameters*/ + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble n = element->GetMaterialParameter(MaterialsRheologyNEnum); + IssmDouble B = element->GetMaterialParameter(MaterialsRheologyBbarEnum); + Input* slopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); + Input* slopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* drag_input = element->GetInput(FrictionCoefficientEnum); _assert_(drag_input); + + Gauss* gauss=element->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + + connectivity=(IssmDouble)element->VertexConnectivity(iv); + + thickness_input->GetInputValue(&thickness,gauss); + drag_input->GetInputValue(&drag,gauss); + slopex_input->GetInputValue(&slope[0],gauss); + slopey_input->GetInputValue(&slope[1],gauss); + slope2=slope[0]*slope[0]+slope[1]*slope[1]; + + /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10 m/s/Pa*/ + ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0]; + vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1]; + ///*Ritz et al. 1996*/ + //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2); + //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2); + ///*Rutt et al. 2009*/ + //ub=-drag*rho_ice*gravity*thickness*slope[0]; + //vb=-drag*rho_ice*gravity*thickness*slope[1]; + + pe->values[2*iv+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity; + pe->values[2*iv+1]=(vb-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/connectivity; + } + + /*Clean up and return*/ + delete gauss; + return pe; +}/*}}}*/ +ElementVector* StressbalanceSIAAnalysis::CreatePVector3D(Element* element){/*{{{*/ + + _error_("not implemented"); }/*}}}*/ void StressbalanceSIAAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h (revision 16848) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.h (revision 16849) @@ -23,4 +23,6 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); + ElementVector* CreatePVector2D(Element* element); + ElementVector* CreatePVector3D(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16848) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16849) @@ -183,4 +183,5 @@ virtual IssmDouble TotalSmb(void)=0; virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0; + virtual int VertexConnectivity(int vertexindex)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16848) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16849) @@ -1022,12 +1022,11 @@ _assert_(this->matpar); - switch(enum_in){ // FIXME: change this to material - case MaterialsRheologyNEnum: - return this->material->GetN(); - case MaterialsRheologyBEnum: - return this->material->GetB(); - default: - return this->matpar->GetMaterialParameter(enum_in); + case MaterialsRheologyNEnum: + return this->material->GetN(); + case MaterialsRheologyBEnum: + return this->material->GetB(); + default: + return this->matpar->GetMaterialParameter(enum_in); } } @@ -3461,4 +3460,10 @@ break; } +} +/*}}}*/ +/*FUNCTION Penta::VertexConnectivity{{{*/ +int Penta::VertexConnectivity(int vertexindex){ + _assert_(this->vertices); + return this->vertices[vertexindex]->Connectivity(); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16848) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16849) @@ -130,4 +130,5 @@ int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum); IssmDouble TimeAdapt(); + int VertexConnectivity(int vertexindex); void ViscousHeatingCreateInput(void); void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16848) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16849) @@ -154,4 +154,5 @@ ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");}; ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");}; + int VertexConnectivity(int vertexindex){_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){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16848) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16849) @@ -1160,5 +1160,14 @@ _assert_(this->matpar); - return this->matpar->GetMaterialParameter(enum_in); + switch(enum_in){ // FIXME: change this to material + case MaterialsRheologyNEnum: + return this->material->GetN(); + case MaterialsRheologyBEnum: + return this->material->GetB(); + case MaterialsRheologyBbarEnum: + return this->material->GetBbar(); + default: + return this->matpar->GetMaterialParameter(enum_in); + } } /*}}}*/ @@ -2937,4 +2946,10 @@ delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::VertexConnectivity{{{*/ +int Tria::VertexConnectivity(int vertexindex){ + _assert_(this->vertices); + return this->vertices[vertexindex]->Connectivity(); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16848) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16849) @@ -137,4 +137,5 @@ void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); IssmDouble TimeAdapt(); + int VertexConnectivity(int vertexindex); void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); bool IsZeroLevelset(int levelset_enum);