Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16806) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16807) @@ -431,4 +431,14 @@ iomodel->FetchData(&z,NULL,NULL,MeshZEnum); switch(finiteelement){ + case P1Enum: + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + if(reCast(vertices_type[i])==NoneApproximationEnum){ + constraints->AddObject(new SpcStatic(count+1,iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,1,g*rho_ice*(surface[i]-z[i])/FSreconditioning,StressbalanceAnalysisEnum)); + count++; + } + } + } + break; case P1bubbleEnum: for(i=0;inumberofvertices;i++){ @@ -909,5 +919,71 @@ ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/ - return NULL; + int i,meshtype,dim; + IssmDouble Jdet,water_pressure,bed; + IssmDouble normal[3]; + IssmDouble *xyz_list_base = NULL; + + /*Get basal element*/ + if(!element->IsOnBed() || !element->IsFloating()) return NULL; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + + /*Prepare coordinate system list*/ + int* cs_list = xNew(vnumnodes+pnumnodes); + if(dim==2) for(i=0;iNewElementVector(FSvelocityEnum); + IssmDouble* vbasis = xNew(vnumnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + Input* bed_input=element->GetInput(BedEnum); _assert_(bed_input); + IssmDouble rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsVelocity(vbasis,gauss); + + element->NormalBase(&normal[0],xyz_list_base); + _assert_(normal[dim-1]<0.); + bed_input->GetInputValue(&bed, gauss); + water_pressure=gravity*rho_water*bed; + + for(i=0;ivalues[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0]; + pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1]; + if(dim==3){ + pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2]; + } + } + } + + /*Transform coordinate system*/ + element->TransformLoadVectorCoord(pe,cs_list); + + /*Clean up and return*/ + delete gauss; + xDelete(cs_list); + xDelete(vbasis); + xDelete(xyz_list_base); + return pe; }/*}}}*/ ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16806) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16807) @@ -56,4 +56,5 @@ virtual void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss)=0; virtual void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0; + virtual void NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0; virtual void TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0; virtual void TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16806) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16807) @@ -162,6 +162,6 @@ /*}}}*/ -/*FUNCTION Penta::BedNormal {{{*/ -void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){ +/*FUNCTION Penta::NormalBase {{{*/ +void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ int i; @@ -171,6 +171,6 @@ 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]; + v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; + v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; } @@ -181,7 +181,7 @@ /*Bed normal is opposite to surface normal*/ - *(bed_normal)=-normal[0]/normal_norm; - *(bed_normal+1)=-normal[1]/normal_norm; - *(bed_normal+2)=-normal[2]/normal_norm; + bed_normal[0]=-normal[0]/normal_norm; + bed_normal[1]=-normal[1]/normal_norm; + bed_normal[2]=-normal[2]/normal_norm; } /*}}}*/ @@ -304,5 +304,5 @@ /*Get normal vector to the bed */ - BedNormal(&bed_normal[0],xyz_list_tria); + NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); /*basalforce*/ @@ -4788,5 +4788,5 @@ /*geothermal heatflux*/ - BedNormal(&normal_base[0],xyz_list_tria); + NormalBase(&normal_base[0],&xyz_list_tria[0][0]); heatflux=0.; for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; @@ -6888,5 +6888,5 @@ material->GetViscosity3dFS(&viscosity,&epsilon[0]); - BedNormal(&bed_normal[0],xyz_list_tria); + NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); @@ -8070,5 +8070,5 @@ vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - BedNormal(&bed_normal[0],xyz_list_tria); + NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -8235,5 +8235,5 @@ vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - BedNormal(&bed_normal[0],xyz_list_tria); + NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); @@ -8934,5 +8934,5 @@ GetNodalFunctionsVelocity(vbasis, gauss); - BedNormal(&bed_normal[0],xyz_list_tria); + NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); bed_input->GetInputValue(&bed, gauss); if(shelf_dampening){ //add dampening to avoid too high vertical velocities when not in hydrostatic equilibrium Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16806) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16807) @@ -192,5 +192,5 @@ void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); void AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum); - void BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]); + void NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list); ElementMatrix* CreateBasalMassMatrix(void); ElementMatrix* CreateKMatrix(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16806) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16807) @@ -114,4 +114,5 @@ void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; bool NoIceInElement(){_error_("not implemented yet");}; + void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; int NumberofNodesPressure(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16806) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16807) @@ -1162,6 +1162,6 @@ }/*}}}*/ -/*FUNCTION Tria::GetSegmentNormal {{{*/ -void Tria:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]){ +/*FUNCTION Tria::NormalBase {{{*/ +void Tria::NormalBase(IssmDouble* normal,IssmDouble* xyz_list){ /*Build unit outward pointing vector*/ @@ -1169,6 +1169,6 @@ IssmDouble norm; - vector[0]=xyz_list[1][0] - xyz_list[0][0]; - vector[1]=xyz_list[1][1] - xyz_list[0][1]; + vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0]; + vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1]; norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]); @@ -1976,5 +1976,15 @@ bool onbed; - inputs->GetInputValue(&onbed,MeshElementonbedEnum); + int meshtype; + this->parameters->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: + return HasEdgeOnBed(); + break; + case Mesh2DhorizontalEnum: + inputs->GetInputValue(&onbed,MeshElementonbedEnum); + break; + default: _error_("mesh "<GetInputValue(&bed, gauss); @@ -4157,5 +4167,5 @@ GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum); GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2); - GetSegmentNormal(&normal[0],xyz_list_front); + NormalBase(&normal[0],&xyz_list_front[0][0]); /*Start looping on Gaussian points*/ @@ -7125,5 +7135,5 @@ this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); for(int i=0;iEdgeOnBedIndices(&indices[0],&indices[1]); for(int i=0;i