Changeset 27660

Timestamp:

03/24/23 14:15:01 (2 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixed 3d problem, now working

File:

• issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (modified) (5 diffs)

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

@@ -539 +539 @@
         if(element->IsAllFloating() || !element->IsIceInElement()) return;
         if(!element->IsOnBase()) return;
+        Element* basalelement = element->SpawnBasalElement();
 
         /*Intermediaries */
@@ -552 +553 @@
 
         /*Retrieve all inputs and parameters*/
-        element->GetVerticesCoordinates(&xyz_list);
-        element->FindParam(&dt,TimesteppingTimeStepEnum);
-        IssmDouble  latentheat      = element->FindParam(MaterialsLatentheatEnum);
-        IssmDouble  g               = element->FindParam(ConstantsGEnum);
-        IssmDouble  rho_ice         = element->FindParam(MaterialsRhoIceEnum);
-        IssmDouble  rho_water       = element->FindParam(MaterialsRhoFreshwaterEnum);
-        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
-        Input* head_input           = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
-        Input* gap_input            = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
-        Input* thickness_input      = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
-        Input* base_input           = element->GetInput(BaseEnum);                       _assert_(base_input);
-        Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
-        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
-        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
-        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+        IssmDouble  latentheat      = basalelement->FindParam(MaterialsLatentheatEnum);
+        IssmDouble  g               = basalelement->FindParam(ConstantsGEnum);
+        IssmDouble  rho_ice         = basalelement->FindParam(MaterialsRhoIceEnum);
+        IssmDouble  rho_water       = basalelement->FindParam(MaterialsRhoFreshwaterEnum);
+        Input* geothermalflux_input = basalelement->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
+        Input* head_input           = basalelement->GetInput(HydrologyHeadEnum);              _assert_(head_input);
+        Input* gap_input            = basalelement->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
+        Input* thickness_input      = basalelement->GetInput(ThicknessEnum);                  _assert_(thickness_input);
+        Input* base_input           = basalelement->GetInput(BaseEnum);                       _assert_(base_input);
+        Input* B_input              = basalelement->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+        Input* n_input              = basalelement->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+        Input* lr_input             = basalelement->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
+        Input* br_input             = basalelement->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
 
         /*Get conductivity from inputs*/
-        IssmDouble conductivity = GetConductivity(element);
-
-        /*Build friction element, needed later: */
-        Friction* friction=new Friction(element,2);
+        IssmDouble conductivity = GetConductivity(basalelement);
+
+        /*Build friction basalelement, needed later: */
+        Friction* friction=new Friction(basalelement,2);
 
         /*Keep track of weights*/
@@ -578 +579 @@
 
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(2);
+        Gauss* gauss=basalelement->NewGauss(2);
         while(gauss->next()){
 
-                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 
                 geothermalflux_input->GetInputValue(&G,gauss);
@@ -621 +622 @@
                 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
 
-                element->AddBasalInput(DummyEnum,&meltrate,P0Enum);
-                element->AddBasalInput(EsaEmotionEnum,&frictionheat,P0Enum);
-                element->AddBasalInput(EsaNmotionEnum,&dissipation,P0Enum);
-                element->AddBasalInput(EsaUmotionEnum,&PMPheat,P0Enum);
+                //element->AddBasalInput(DummyEnum,&meltrate,P0Enum);
+                //element->AddBasalInput(EsaEmotionEnum,&frictionheat,P0Enum);
+                //element->AddBasalInput(EsaNmotionEnum,&dissipation,P0Enum);
+                //element->AddBasalInput(EsaUmotionEnum,&PMPheat,P0Enum);
 
                 newgap += gauss->weight*Jdet*(gap+dt*(
@@ -666 +667 @@
         delete friction;
         delete gauss;
+        if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void HydrologyShaktiAnalysis::UpdateEffectivePressure(FemModel* femmodel){/*{{{*/