source: issm/oecreview/Archive/16554-17801/ISSM-17393-17394.diff@ 17802

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -6 +6 @@
 #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){/*{{{*/

../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -532 +532 @@
         StressTensoryyEnum,
         StressTensoryzEnum,
         StressTensorzzEnum,
+        DivergenceEnum,
         GiaCrossSectionShapeEnum,
         GiadWdtEnum,
         GiaWEnum,

../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -521 +521 @@
                 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";

../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -533 +533 @@
               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 +628 @@
               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;

../trunk-jpl/src/c/classes/FemModel.cpp

@@ -431 +431 @@
 
         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 +504 @@
                         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 +1043 @@
         *pSmb=total_smb;
 
 }/*}}}*/
+void FemModel::Divergencex(IssmDouble* pdiv){/*{{{*/
+
+        IssmDouble local_divergence=0;
+        IssmDouble total_divergence;
+
+        for(int i=0;i<this->elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(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;

../trunk-jpl/src/c/classes/Elements/Element.h

@@ -55 +55 @@
                 /* 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);

../trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -129 +129 @@
 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();ig<gauss->end();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);
 }/*}}}*/

../trunk-jpl/src/c/classes/FemModel.h

@@ -71 +71 @@
                 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);