Unifying the standard pressure, PorousPressure, and DamagedPressure p…

…olicies

src/CRKSPH/SolidCRKSPHHydroBase.cc

@@ -26,7 +26,6 @@
 #include "DataBase/IncrementBoundedState.hh"
 #include "DataBase/ReplaceState.hh"
 #include "DataBase/ReplaceBoundedState.hh"
-#include "Damage/DamagedPressurePolicy.hh"
 #include "ArtificialViscosity/ArtificialViscosity.hh"
 #include "DataBase/DataBase.hh"
 #include "Field/FieldList.hh"
@@ -203,9 +202,7 @@ registerState(DataBase<Dimension>& dataBase,
 
   // Grab the normal Hydro's registered version of the sound speed and pressure.
   auto cs = state.fields(HydroFieldNames::soundSpeed, 0.0);
-  auto P = state.fields(HydroFieldNames::pressure, 0.0);
   CHECK(cs.numFields() == dataBase.numFluidNodeLists());
-  CHECK(P.numFields() == dataBase.numFluidNodeLists());
 
   // Register the deviatoric stress and plastic strain to be evolved.
   auto ps = dataBase.solidPlasticStrain();
@@ -220,7 +217,6 @@ registerState(DataBase<Dimension>& dataBase,
 
   // Override the policies for the sound speed and pressure.
   state.enroll(cs, std::make_shared<StrengthSoundSpeedPolicy<Dimension>>());
-  state.enroll(P, std::make_shared<DamagedPressurePolicy<Dimension>>());
 
   // Register the damage with a default no-op update.
   // If there are any damage models running they can override this choice.

src/Damage/CMakeLists.txt

@@ -17,7 +17,6 @@ set(Damage_inst
     JohnsonCookDamage
     JohnsonCookFailureStrainPolicy
     JohnsonCookDamagePolicy
-    DamagedPressurePolicy
     PairMaxDamageNodeCoupling
     ThreePointDamagedNodeCoupling
     DamageGradientNodeCoupling
@@ -31,8 +30,6 @@ set(Damage_headers
     DamageGradientPolicy.hh
     DamageModel.hh
     DamageModelInline.hh
-    DamagedPressurePolicy.hh
-    DamagedSoundSpeedPolicy.hh
     EffectiveTensorDamagePolicy.hh
     GradyKippScalarDamage.hh
     JohnsonCookDamage.hh

src/Hydro/PressurePolicy.cc

@@ -6,13 +6,15 @@
 //----------------------------------------------------------------------------//
 
 #include "PressurePolicy.hh"
-#include "HydroFieldNames.hh"
+#include "Hydro/HydroFieldNames.hh"
+#include "Strength/SolidFieldNames.hh"
 #include "FSISPH/FSIFieldNames.hh"
 #include "DataBase/State.hh"
 #include "DataBase/StateDerivatives.hh"
 #include "Field/Field.hh"
 #include "NodeList/FluidNodeList.hh"
 #include "Material/EquationOfState.hh"
+#include "SolidMaterial/SolidEquationOfState.hh"
 #include "Utilities/DBC.hh"
 
 namespace Spheral {
@@ -54,17 +56,66 @@ update(const KeyType& key,
            fieldKey == FSIFieldNames::rawPressure));
   auto& P = state.field(key, Scalar());
 
-  // Get the mass density and specific thermal energy fields from the state.
-  const auto& massDensity = state.field(State<Dimension>::buildFieldKey(HydroFieldNames::massDensity, nodeListKey), Scalar());
-  const auto& eps = state.field(State<Dimension>::buildFieldKey(HydroFieldNames::specificThermalEnergy, nodeListKey), Scalar());
-
   // Get the eos.  This cast is ugly, but is a work-around for now.
   const auto* fluidNodeListPtr = dynamic_cast<const FluidNodeList<Dimension>*>(P.nodeListPtr());
   CHECK(fluidNodeListPtr != nullptr);
   const auto& eos = fluidNodeListPtr->equationOfState();
 
-  // Now set the pressure for this field.
-  eos.setPressure(P, massDensity, eps);
+  // Grab some of our required state
+  const auto  buildKey = [&](const std::string& fkey) { return StateBase<Dimension>::buildFieldKey(fkey, nodeListKey); };
+  const auto& eps = state.field(buildKey(HydroFieldNames::specificThermalEnergy), 0.0);
+
+  // Check if this material has porosity.
+  if (state.registered(buildKey(SolidFieldNames::porosityAlpha))) {
+
+    // Yep, there's porosity
+    const auto& rhoS = state.field(buildKey(SolidFieldNames::porositySolidDensity), 0.0);
+    const auto& alpha = state.field(buildKey(SolidFieldNames::porosityAlpha), 0.0);
+
+    // Check if we need to update the pressure derivatives by seeing if they're registered for this NodeList
+    const auto dPduKey = State<Dimension>::buildFieldKey(HydroFieldNames::partialPpartialEps, nodeListKey);
+    const auto dPdrKey = State<Dimension>::buildFieldKey(HydroFieldNames::partialPpartialRho, nodeListKey);
+    if (state.registered(dPduKey)) {
+      CHECK(state.registered(dPdrKey));
+      auto& dPdu = state.field(dPduKey, 0.0);
+      auto& dPdr = state.field(dPdrKey, 0.0);
+      eos.setPressureAndDerivs(P, dPdu, dPdr, rhoS, eps);
+    } else {
+      eos.setPressure(P, rhoS, eps);
+    }
+
+    // Scale the pressure to the bulk value.
+    P /= alpha;
+
+  } else {
+
+    // No porosity, so just set the pressure using the straight density
+    const auto& rho = state.field(buildKey(HydroFieldNames::massDensity), 0.0);
+    eos.setPressure(P, rho, eps);
+
+  }
+
+  // If there's damage for this material, apply it to the pressure
+  if (state.registered(buildKey(SolidFieldNames::tensorDamage))) {
+    const auto& D = state.field(buildKey(SolidFieldNames::tensorDamage), SymTensor::zero);
+
+    // Check for the appropriate minimum pressure
+    const auto* solidEOSptr = dynamic_cast<const SolidEquationOfState<Dimension>*>(&eos);
+    const auto Pmin = eos.minimumPressure();
+    const auto PminDamage = (solidEOSptr != nullptr ?
+                             solidEOSptr->minimumPressureDamage() :
+                             Pmin);
+
+    // Scale by the damage.
+    const auto ni = P.numInternalElements();
+#pragma omp parallel for
+    for (auto i = 0u; i < ni; ++i) {
+      const auto Di = std::max(0.0, std::min(1.0, D(i).eigenValues().maxElement()));
+      CHECK(Di >= 0.0 and Di <= 1.0);
+      P(i) = std::max(Pmin, (1.0 - Di)*P(i)) + std::max(PminDamage, Di*P(i));
+    }
+  }
+
 }
 
 //------------------------------------------------------------------------------

src/Hydro/PressurePolicy.hh

@@ -25,6 +25,7 @@ public:
   //--------------------------- Public Interface ---------------------------//
   // Useful typedefs
   using Scalar = typename Dimension::Scalar;
+  using SymTensor = typename Dimension::SymTensor;
   using KeyType = typename FieldUpdatePolicy<Dimension>::KeyType;
 
   // Constructors, destructor.

src/Porosity/CMakeLists.txt

@@ -4,7 +4,6 @@ set(Porosity_inst
     StrainPorosity
     PalphaPorosity
     PorositySolidMassDensityPolicy
-    PorousPressurePolicy
     PorousSoundSpeedPolicy
     PorousGammaPolicy
     PorousBulkModulusPolicy
@@ -24,7 +23,6 @@ set(Porosity_headers
     StrainPorosity.hh
     StrainPorosityInline.hh
     PalphaPorosity.hh
-    PorousPressurePolicy.hh
     PorousSoundSpeedPolicy.hh
     PorousGammaPolicy.hh
     PorousBulkModulusPolicy.hh

src/Porosity/PorosityModel.cc

@@ -7,7 +7,6 @@
 #include "Porosity/PorosityModel.hh"
 #include "Porosity/PorositySolidMassDensityPolicy.hh"
 #include "Porosity/PorousBulkModulusPolicy.hh"
-#include "Porosity/PorousPressurePolicy.hh"
 #include "Porosity/PorousSoundSpeedPolicy.hh"
 #include "Porosity/PorousGammaPolicy.hh"
 #include "Porosity/PorousEntropyPolicy.hh"
@@ -170,10 +169,6 @@ registerState(DataBase<Dimension>& dataBase,
   // Initial sound speed
   state.enroll(mc0);
 
-  // Override the pressure policy
-  auto& P = state.field(buildKey(HydroFieldNames::pressure), 0.0);
-  state.enroll(P, std::make_shared<PorousPressurePolicy<Dimension>>());
-
   // Check what other state is registered which needs to be overridden for porosity
   auto optionalOverridePolicy = [&](const std::string& fkey, std::shared_ptr<UpdatePolicyBase<Dimension>> policy) -> void {
                                   const auto fullkey = buildKey(fkey);