Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16838) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16839) @@ -1136,5 +1136,5 @@ element->GetVerticesCoordinates(&xyz_list); element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); - element->NormalBase(&normal[0],xyz_list_front); + element->NormalSection(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ @@ -1570,5 +1570,66 @@ ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/ - return NULL; + /*If no front, return NULL*/ + if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL; + + /*Intermediaries*/ + IssmDouble rho_ice,rho_water,gravity; + IssmDouble Jdet,surface,z_g,water_pressure,ice_pressure; + IssmDouble surface_under_water,base_under_water,pressure; + IssmDouble* xyz_list = NULL; + IssmDouble* xyz_list_front = NULL; + IssmDouble normal[3]; + Gauss* gauss = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = element->NewElementVector(HOApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input); + rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + gravity = element->GetMaterialParameter(ConstantsGEnum); + element->GetVerticesCoordinates(&xyz_list); + element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum); + element->NormalSection(&normal[0],xyz_list_front); + + /*Initialize gauss points*/ + IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<6;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2]; + IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<6;i++) if(xyz_list[i*3+2]0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity + else gauss=element->NewGauss(xyz_list,xyz_list_front,3,3); + + /* Start looping on the number of gaussian points: */ + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + surface_input->GetInputValue(&surface,gauss); + z_g=element->GetZcoord(gauss); + element->NodalFunctions(basis,gauss); + element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); + + water_pressure = rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level + ice_pressure = rho_ice*gravity*(surface-z_g); + pressure = ice_pressure + water_pressure; + + for (int i=0;ivalues[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; + pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,XYEnum); + + /*Clean up and return*/ + xDelete(basis); + xDelete(xyz_list); + xDelete(xyz_list_front); + delete gauss; + return pe; }/*}}}*/ void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ @@ -1937,5 +1998,5 @@ element->NodalFunctionsVelocity(vbasis,gauss); - element->NormalBase(&normal[0],xyz_list_base); + element->NormalSection(&normal[0],xyz_list_base); _assert_(normal[dim-1]<0.); bed_input->GetInputValue(&bed, gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16838) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16839) @@ -62,5 +62,5 @@ virtual void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0; virtual void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0; - virtual void NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0; + virtual void NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0; virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0; virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0; @@ -111,4 +111,5 @@ virtual void GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0; virtual void GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0; + virtual IssmDouble GetZcoord(Gauss* gauss)=0; virtual void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; @@ -132,4 +133,5 @@ virtual Gauss* NewGauss(int order)=0; virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0; + 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* NewGaussTop(int order)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16838) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16839) @@ -1511,6 +1511,6 @@ } /*}}}*/ -/*FUNCTION Penta::GetQuadNormal {{{*/ -void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list){ +/*FUNCTION Penta::NormalSection{{{*/ +void Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){ /*Build unit outward pointing vector*/ @@ -1634,5 +1634,5 @@ /*}}}*/ /*FUNCTION Penta::GetZcoord {{{*/ -IssmDouble Penta::GetZcoord(GaussPenta* gauss){ +IssmDouble Penta::GetZcoord(Gauss* gauss){ int i; @@ -2466,4 +2466,12 @@ } /*}}}*/ +/*FUNCTION Penta::JacobianDeterminantSurface{{{*/ +void Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussPentaEnum); + this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss); + +} +/*}}}*/ /*FUNCTION Penta::NoIceInElement {{{*/ bool Penta::NoIceInElement(){ @@ -2542,4 +2550,14 @@ Gauss* Penta::NewGauss(int order){ return new GaussPenta(order,order); +} +/*}}}*/ +/*FUNCTION Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){{{*/ +Gauss* Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){ + + IssmDouble area_coordinates[4][3]; + + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4); + + return new GaussPenta(area_coordinates,order_horiz,order_vert); } /*}}}*/ @@ -3123,28 +3141,4 @@ return S; } -} -/*}}}*/ -/*FUNCTION Penta::SurfaceNormal {{{*/ -void Penta::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){ - - int i; - IssmDouble v13[3],v23[3]; - IssmDouble normal[3]; - IssmDouble normal_norm; - - for (i=0;i<3;i++){ - v13[i]=xyz_list[0][i]-xyz_list[2][i]; - v23[i]=xyz_list[1][i]-xyz_list[2][i]; - } - - normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; - normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; - normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - - normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) ); - - *(surface_normal)=normal[0]/normal_norm; - *(surface_normal+1)=normal[1]/normal_norm; - *(surface_normal+2)=normal[2]/normal_norm; } /*}}}*/ @@ -6072,5 +6066,5 @@ /*Get normal vector to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list_tria); + NormalTop(&surface_normal[0],&xyz_list_tria[0][0]); bed_normal[0]=-surface_normal[0]; //Function is for upper surface, so the normal to the bed is the opposite of the result @@ -8123,5 +8117,5 @@ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); for(i=0;i(basis); + xDelete(xyz_list_front); delete gauss; return pe; @@ -8883,5 +8878,5 @@ ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4); - GetQuadNormal(&normal[0],xyz_list_front); + NormalSection(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16838) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16839) @@ -97,5 +97,5 @@ IssmDouble GetMaterialParameter(int enum_in); void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); - IssmDouble GetZcoord(GaussPenta* gauss); + IssmDouble GetZcoord(Gauss* gauss); void GetVectorFromInputs(Vector* vector,int name_enum); void GetVerticesCoordinates(IssmDouble** pxyz_list); @@ -199,4 +199,5 @@ void AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum); void NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list); + void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); void NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list); ElementMatrix* CreateBasalMassMatrix(void); @@ -232,5 +233,4 @@ void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); - void GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list); void GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input); void GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input); @@ -247,5 +247,5 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags); void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); - void JacobianDeterminantSurface(IssmDouble* pJdet, 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); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); @@ -254,4 +254,5 @@ Gauss* NewGauss(int order); Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); Gauss* NewGaussBase(int order); Gauss* NewGaussTop(int order); @@ -266,5 +267,4 @@ void SetClone(int* minranks); Tria* SpawnTria(int location); - void SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); void TransformLoadVectorCoord(ElementVector* pe,int transformenum); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16838) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16839) @@ -124,5 +124,5 @@ void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; bool NoIceInElement(){_error_("not implemented yet");}; - void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; + void NormalSection(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; @@ -141,7 +141,9 @@ void GetInputValue(IssmDouble* pvalue,int enum_type){_error_("not implemented yet");}; void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");}; + IssmDouble GetZcoord(Gauss* gauss){_error_("not implemented yet");}; Gauss* NewGauss(void){_error_("not implemented yet");}; Gauss* NewGauss(int order){_error_("not implemented yet");}; Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; + 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* NewGaussTop(int order){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16838) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16839) @@ -1168,6 +1168,6 @@ }/*}}}*/ -/*FUNCTION Tria::NormalBase {{{*/ -void Tria::NormalBase(IssmDouble* normal,IssmDouble* xyz_list){ +/*FUNCTION Tria::NormalSection{{{*/ +void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){ /*Build unit outward pointing vector*/ @@ -2697,27 +2697,4 @@ } /*}}}*/ -/*FUNCTION Tria::SurfaceNormal{{{*/ -void Tria::SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]){ - - IssmDouble v13[3],v23[3]; - IssmDouble normal[3]; - IssmDouble normal_norm; - - for(int i=0;i<3;i++){ - v13[i]=xyz_list[0][i]-xyz_list[2][i]; - v23[i]=xyz_list[1][i]-xyz_list[2][i]; - } - - normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; - normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; - normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - - normal_norm=sqrt( normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); - - *(surface_normal+0) = normal[0]/normal_norm; - *(surface_normal+1) = normal[1]/normal_norm; - *(surface_normal+2) = normal[2]/normal_norm; -} -/*}}}*/ /*FUNCTION Tria::TimeAdapt{{{*/ IssmDouble Tria::TimeAdapt(void){ @@ -3976,5 +3953,5 @@ ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2); - NormalBase(&normal[0],xyz_list_front); + NormalSection(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ @@ -4122,5 +4099,5 @@ GetNodalFunctionsVelocity(vbasis, gauss); - NormalBase(&normal[0],&xyz_list_seg[0][0]); + NormalSection(&normal[0],&xyz_list_seg[0][0]); _assert_(normal[1]<0.); bed_input->GetInputValue(&bed, gauss); @@ -4242,5 +4219,5 @@ ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum); GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2); - NormalBase(&normal[0],xyz_list_front); + NormalSection(&normal[0],xyz_list_front); /*Start looping on Gaussian points*/ @@ -7220,5 +7197,5 @@ this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); for(int i=0;iEdgeOnBedIndices(&indices[0],&indices[1]); for(int i=0;i