Changeset 17394

Timestamp:

03/08/14 13:37:59 (11 years ago)

Author:

seroussi

Message:

NEW: added divergence calculation in elements

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (1 diff)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (3 diffs)
• classes/FemModel.h (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (2 diffs)

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

@@ -7 +7 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 101
+//#define FSANALYTICAL 4
 
 /*Model processing*/

issm/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){/*{{{*/

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

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

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

@@ -432 +432 @@
         switch (response_descriptor_enum){
 
+                case DivergenceEnum:               this->Divergencex(responses); break;
                 case IceVolumeEnum:                this->IceVolumex(responses); break;
                 case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break;
@@ -504 +505 @@
 
                                 /*Scalar output*/
+                                case DivergenceEnum:               this->Divergencex(&double_result);               break;
                                 case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
                                 case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
@@ -1042 +1044 @@
 
 }/*}}}*/
+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){/*{{{*/
 

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

@@ -72 +72 @@
                 void MinVzx(IssmDouble* presponse);
                 void TotalSmbx(IssmDouble* pSmb);
+                void Divergencex(IssmDouble* pdiv);
                 void IceVolumex(IssmDouble* pV);
                 void IceVolumeAboveFloatationx(IssmDouble* pV);

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

@@ -533 +533 @@
         StressTensoryzEnum,
         StressTensorzzEnum,
+        DivergenceEnum,
         GiaCrossSectionShapeEnum,
         GiadWdtEnum,

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

@@ -522 +522 @@
                 case StressTensoryzEnum : return "StressTensoryz";
                 case StressTensorzzEnum : return "StressTensorzz";
+                case DivergenceEnum : return "Divergence";
                 case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape";
                 case GiadWdtEnum : return "GiadWdt";

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

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