Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16849) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16850) @@ -816,4 +816,6 @@ element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ + case SIAApproximationEnum: + return NULL; case SSAApproximationEnum: return CreateKMatrixSSA(element); @@ -833,4 +835,6 @@ element->GetInputValue(&approximation,ApproximationEnum); switch(approximation){ + case SIAApproximationEnum: + return NULL; case SSAApproximationEnum: return CreatePVectorSSA(element); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16849) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (revision 16850) @@ -204,5 +204,91 @@ ElementVector* StressbalanceSIAAnalysis::CreatePVector3D(Element* element){/*{{{*/ - _error_("not implemented"); + /*Intermediaries */ + int nodeup,nodedown,numsegments; + IssmDouble ub,vb,slope2,drag,surface,thickness,constant_part,z,Jdet; + IssmDouble slope[2],connectivity[2],xyz_list_line[2][3]; + IssmDouble *xyz_list = NULL; + int *pairindices = NULL; + + /*Fetch number vertices for this element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize Element vector*/ + ElementVector* pe=element->NewElementVector(); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble n = element->GetMaterialParameter(MaterialsRheologyNEnum); + IssmDouble B = element->GetMaterialParameter(MaterialsRheologyBEnum); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + 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); + + /*Get Vertical segment indices*/ + element->VerticalSegmentIndices(&pairindices,&numsegments); + for(int is=0;isVertexConnectivity(nodedown); + connectivity[1]=(IssmDouble)element->VertexConnectivity(nodeup); + for(int i=0;i<3;i++){ + xyz_list_line[0][i]=xyz_list[nodedown*3+i]; + xyz_list_line[1][i]=xyz_list[nodeup*3+i]; + } + + Gauss* gauss=element->NewGaussLine(nodedown,nodeup,3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + slopex_input->GetInputValue(&slope[0],gauss); + slopey_input->GetInputValue(&slope[1],gauss); + surface_input->GetInputValue(&surface,gauss); + thickness_input->GetInputValue(&thickness,gauss); + + slope2=slope[0]*slope[0]+slope[1]*slope[1]; + constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.)); + + z = element->GetZcoord(gauss); + element->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss); + + if(element->IsOnSurface()){ + pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1]; + pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1]; + } + else{/*connectivity is too large, should take only half on it*/ + pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1]; + pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1]; + } + } + + /*Deal with basal velocities*/ + if(element->IsOnBed()){ + drag_input->GetInputValue(&drag,gauss); + + /*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*nodedown+0]+=ub/connectivity[0]; + pe->values[2*nodedown+1]+=vb/connectivity[0]; + } + delete gauss; + } + + /*Clean up and return*/ + xDelete(pairindices); + xDelete(xyz_list); + return pe; + }/*}}}*/ void StressbalanceSIAAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16849) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16850) @@ -85,4 +85,5 @@ 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; virtual void JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0; virtual void JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0; @@ -140,4 +141,5 @@ virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert)=0; virtual Gauss* NewGaussBase(int order)=0; + virtual Gauss* NewGaussLine(int vertex1,int vertex2,int order)=0; virtual Gauss* NewGaussTop(int order)=0; virtual ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum)=0; @@ -184,4 +186,5 @@ virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0; virtual int VertexConnectivity(int vertexindex)=0; + virtual void VerticalSegmentIndices(int** pindices,int* pnumseg)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16849) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16850) @@ -2467,4 +2467,12 @@ } /*}}}*/ +/*FUNCTION Penta::JacobianDeterminantLine{{{*/ +void Penta::JacobianDeterminantLine(IssmDouble* pJdet,IssmDouble* xyz_list_line,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetSegmentJacobianDeterminant(pJdet,xyz_list_line,(GaussPenta*)gauss); + +} +/*}}}*/ /*FUNCTION Penta::JacobianDeterminantTop{{{*/ void Penta::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){ @@ -2574,4 +2582,9 @@ Gauss* Penta::NewGaussBase(int order){ return new GaussPenta(0,1,2,order); +} +/*}}}*/ +/*FUNCTION Penta::NewGaussLine(int vertex1,int vertex2,int order){{{*/ +Gauss* Penta::NewGaussLine(int vertex1,int vertex2,int order){ + return new GaussPenta(vertex1,vertex2,order); } /*}}}*/ @@ -3466,4 +3479,18 @@ _assert_(this->vertices); return this->vertices[vertexindex]->Connectivity(); +} +/*}}}*/ +/*FUNCTION Penta::VerticalSegmentIndices{{{*/ +void Penta::VerticalSegmentIndices(int** pindices,int* pnumseg){ + + /*output*/ + int *indices = xNew(3*2); + indices[0*2 + 0] = 0; indices[0*2 + 1] = 3; + indices[1*2 + 0] = 1; indices[1*2 + 1] = 4; + indices[2*2 + 0] = 2; indices[2*2 + 1] = 5; + + /*Assign output pointers*/ + *pindices = indices; + *pnumseg = 3; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16849) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16850) @@ -131,4 +131,5 @@ IssmDouble TimeAdapt(); int VertexConnectivity(int vertexindex); + void VerticalSegmentIndices(int** pindices,int* pnumseg); void ViscousHeatingCreateInput(void); void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); @@ -250,4 +251,5 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags); void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); + void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); @@ -259,4 +261,5 @@ Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); Gauss* NewGaussBase(int order); + Gauss* NewGaussLine(int vertex1,int vertex2,int order); Gauss* NewGaussTop(int order); ElementVector* NewElementVector(int approximation_enum); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16849) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16850) @@ -116,5 +116,6 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; - void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; @@ -151,8 +152,10 @@ Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");}; Gauss* NewGaussBase(int order){_error_("not implemented yet");}; + Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order){_error_("not implemented yet");}; 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 VerticalSegmentIndices(int** pindices,int* pnumseg){_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.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16849) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16850) @@ -138,4 +138,5 @@ IssmDouble TimeAdapt(); int VertexConnectivity(int vertexindex); + void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum); bool IsZeroLevelset(int levelset_enum); @@ -283,4 +284,5 @@ bool IsInput(int name); void JacobianDeterminant(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); + void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; @@ -292,4 +294,5 @@ Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");}; Gauss* NewGaussBase(int order){_error_("not implemented yet");}; + Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order){_error_("not implemented yet");}; ElementVector* NewElementVector(int approximation_enum);