Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17393) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17394) @@ -6,7 +6,7 @@ #include "../solutionsequences/solutionsequences.h" #include "../cores/cores.h" -//#define FSANALYTICAL 101 +//#define FSANALYTICAL 4 /*Model processing*/ int StressbalanceAnalysis::DofsPerNode(int** pdoftype,int meshtype,int approximation){/*{{{*/ Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17393) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17394) @@ -532,6 +532,7 @@ StressTensoryyEnum, StressTensoryzEnum, StressTensorzzEnum, + DivergenceEnum, GiaCrossSectionShapeEnum, GiadWdtEnum, GiaWEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17393) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17394) @@ -521,6 +521,7 @@ case StressTensoryyEnum : return "StressTensoryy"; case StressTensoryzEnum : return "StressTensoryz"; case StressTensorzzEnum : return "StressTensorzz"; + case DivergenceEnum : return "Divergence"; case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape"; case GiadWdtEnum : return "GiadWdt"; case GiaWEnum : return "GiaW"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17393) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17394) @@ -533,6 +533,7 @@ else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; + else if (strcmp(name,"Divergence")==0) return DivergenceEnum; else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; else if (strcmp(name,"GiaW")==0) return GiaWEnum; @@ -627,11 +628,11 @@ else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum; else if (strcmp(name,"Regular")==0) return RegularEnum; - else if (strcmp(name,"Scaled")==0) return ScaledEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Separate")==0) return SeparateEnum; + if (strcmp(name,"Scaled")==0) return ScaledEnum; + else if (strcmp(name,"Separate")==0) return SeparateEnum; else if (strcmp(name,"Sset")==0) return SsetEnum; else if (strcmp(name,"Verbose")==0) return VerboseEnum; else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 17393) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 17394) @@ -431,6 +431,7 @@ switch (response_descriptor_enum){ + case DivergenceEnum: this->Divergencex(responses); break; case IceVolumeEnum: this->IceVolumex(responses); break; case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break; case MinVelEnum: this->MinVelx(responses); break; @@ -503,6 +504,7 @@ switch(output_enum){ /*Scalar output*/ + case DivergenceEnum: this->Divergencex(&double_result); break; case IceVolumeEnum: this->IceVolumex(&double_result); break; case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break; case MinVelEnum: this->MinVelx(&double_result); break; @@ -1041,6 +1043,22 @@ *pSmb=total_smb; }/*}}}*/ +void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/ + + IssmDouble local_divergence=0; + IssmDouble total_divergence; + + for(int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + local_divergence+=element->Divergence(); + } + ISSM_MPI_Reduce(&local_divergence,&total_divergence,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_divergence,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pdiv=total_divergence; + +}/*}}}*/ void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/ IssmDouble local_ice_volume = 0; Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17393) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17394) @@ -55,6 +55,7 @@ /* bool AnyActive(void); */ void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array); void DeleteMaterials(void); + IssmDouble Divergence(void); void FindParam(bool* pvalue,int paramenum); void FindParam(int* pvalue,int paramenum); void FindParam(IssmDouble* pvalue,int paramenum); Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17393) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17394) @@ -129,6 +129,38 @@ void Element::DeleteMaterials(void){/*{{{*/ delete this->material; }/*}}}*/ +IssmDouble Element::Divergence(void){/*{{{*/ + /*Compute element divergence*/ + + /*Intermediaries*/ + IssmDouble Jdet; + IssmDouble divergence=0.; + IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble *xyz_list = NULL; + + /*Get inputs and parameters*/ + Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = this->GetInput(VzEnum); _assert_(vz_input); + this->GetVerticesCoordinates(&xyz_list); + + Gauss* gauss=this->NewGaussBase(5); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + this->JacobianDeterminant(&Jdet,xyz_list,gauss); + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + + divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet; + } + + /*Clean up and return*/ + delete gauss; + return divergence; +}/*}}}*/ void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/ this->parameters->FindParam(pvalue,paramenum); }/*}}}*/ Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 17393) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 17394) @@ -71,6 +71,7 @@ void MinVyx(IssmDouble* presponse); void MinVzx(IssmDouble* presponse); void TotalSmbx(IssmDouble* pSmb); + void Divergencex(IssmDouble* pdiv); void IceVolumex(IssmDouble* pV); void IceVolumeAboveFloatationx(IssmDouble* pV); void ElementResponsex(IssmDouble* presponse,int response_enum);