Merge pull request #5960 from akva2/add_mech_container

Add MechContainer

CMakeLists_files.cmake

@@ -107,6 +107,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/flow/LogOutputHelper.cpp
   opm/simulators/flow/Main.cpp
   opm/simulators/flow/MainDispatchDynamic.cpp
+  opm/simulators/flow/MechContainer.cpp
   opm/simulators/flow/MixingRateControls.cpp
   opm/simulators/flow/NonlinearSolver.cpp
   opm/simulators/flow/partitionCells.cpp
@@ -856,6 +857,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/KeywordValidation.hpp
   opm/simulators/flow/LogOutputHelper.hpp
   opm/simulators/flow/Main.hpp
+  opm/simulators/flow/MechContainer.hpp
   opm/simulators/flow/MixingRateControls.hpp
   opm/simulators/flow/NewTranFluxModule.hpp
   opm/simulators/flow/NonlinearSolver.hpp

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -513,22 +513,11 @@ assignToSolution(data::Solution& sol)
     const auto extendedSolutionArrays = std::array {
         DataEntry{"BIOFILM",  UnitSystem::measure::identity,           cBiofilm_},
         DataEntry{"CALCITE",  UnitSystem::measure::identity,           cCalcite_},
-        DataEntry{"DELSTRXX", UnitSystem::measure::pressure,           delstressXX_},
-        DataEntry{"DELSTRYY", UnitSystem::measure::pressure,           delstressYY_},
-        DataEntry{"DELSTRZZ", UnitSystem::measure::pressure,           delstressZZ_},
-        DataEntry{"DELSTRXY", UnitSystem::measure::pressure,           delstressXY_},
-        DataEntry{"DELSTRXZ", UnitSystem::measure::pressure,           delstressXZ_},
-        DataEntry{"DELSTRYZ", UnitSystem::measure::pressure,           delstressYZ_},
-        DataEntry{"DISPX",    UnitSystem::measure::length,             dispX_},
-        DataEntry{"DISPY",    UnitSystem::measure::length,             dispY_},
-        DataEntry{"DISPZ",    UnitSystem::measure::length,             dispZ_},
         DataEntry{"DRSDTCON", UnitSystem::measure::gas_oil_ratio_rate, drsdtcon_},
-        DataEntry{"MECHPOTF", UnitSystem::measure::pressure,           mechPotentialForce_},
         DataEntry{"MICROBES", UnitSystem::measure::density,            cMicrobes_},
         DataEntry{"OXYGEN",   UnitSystem::measure::density,            cOxygen_},
         DataEntry{"PERMFACT", UnitSystem::measure::identity,           permFact_},
         DataEntry{"PORV_RC",  UnitSystem::measure::identity,           rockCompPorvMultiplier_},
-        DataEntry{"PRESPOTF", UnitSystem::measure::pressure,           mechPotentialPressForce_},
         DataEntry{"PRES_OVB", UnitSystem::measure::pressure,           overburdenPressure_},
         DataEntry{"RSW",      UnitSystem::measure::gas_oil_ratio,      rsw_},
         DataEntry{"RSWSAT",   UnitSystem::measure::gas_oil_ratio,      gasDissolutionFactorInWater_},
@@ -542,31 +531,6 @@ assignToSolution(data::Solution& sol)
         DataEntry{"STD_CO2",  UnitSystem::measure::identity,           mFracCo2_},
         DataEntry{"STD_GAS",  UnitSystem::measure::identity,           mFracGas_},
         DataEntry{"STD_OIL",  UnitSystem::measure::identity,           mFracOil_},
-        DataEntry{"STRAINXX", UnitSystem::measure::identity,           strainXX_},
-        DataEntry{"STRAINYY", UnitSystem::measure::identity,           strainYY_},
-        DataEntry{"STRAINZZ", UnitSystem::measure::identity,           strainZZ_},
-        DataEntry{"STRAINXY", UnitSystem::measure::identity,           strainXY_},
-        DataEntry{"STRAINXZ", UnitSystem::measure::identity,           strainXZ_},
-        DataEntry{"STRAINYZ", UnitSystem::measure::identity,           strainYZ_},
-        DataEntry{"STRESSXX", UnitSystem::measure::pressure,           stressXX_},
-        DataEntry{"STRESSYY", UnitSystem::measure::pressure,           stressYY_},
-        DataEntry{"STRESSZZ", UnitSystem::measure::pressure,           stressZZ_},
-        DataEntry{"STRESSXY", UnitSystem::measure::pressure,           stressXY_},
-        DataEntry{"STRESSXZ", UnitSystem::measure::pressure,           stressXZ_},
-        DataEntry{"STRESSYZ", UnitSystem::measure::pressure,           stressYZ_},
-        DataEntry{"LINSTRXX", UnitSystem::measure::pressure,           linstressXX_},
-        DataEntry{"LINSTRYY", UnitSystem::measure::pressure,           linstressYY_},
-        DataEntry{"LINSTRZZ", UnitSystem::measure::pressure,           linstressZZ_},
-        DataEntry{"LINSTRXY", UnitSystem::measure::pressure,           linstressXY_},
-        DataEntry{"LINSTRXZ", UnitSystem::measure::pressure,           linstressXZ_},
-        DataEntry{"LINSTRYZ", UnitSystem::measure::pressure,           linstressYZ_},
-        DataEntry{"FRCSTRXX", UnitSystem::measure::pressure,           fracstressXX_},
-        DataEntry{"FRCSTRYY", UnitSystem::measure::pressure,           fracstressYY_},
-        DataEntry{"FRCSTRZZ", UnitSystem::measure::pressure,           fracstressZZ_},
-        DataEntry{"FRCSTRXY", UnitSystem::measure::pressure,           fracstressXY_},
-        DataEntry{"FRCSTRXZ", UnitSystem::measure::pressure,           fracstressXZ_},
-        DataEntry{"FRCSTRYZ", UnitSystem::measure::pressure,           fracstressYZ_},
-        DataEntry{"TEMPPOTF", UnitSystem::measure::pressure,           mechPotentialTempForce_},
         DataEntry{"TMULT_RC", UnitSystem::measure::identity,           rockCompTransMultiplier_},
         DataEntry{"UREA",     UnitSystem::measure::density,            cUrea_},
     };
@@ -617,6 +581,8 @@ assignToSolution(data::Solution& sol)
         doInsert(array, data::TargetType::RESTART_OPM_EXTENDED);
     }
 
+    this->mech_.outputRestart(sol);
+
     if (! this->temperature_.empty())
     {
         sol.insert("TEMP", UnitSystem::measure::temperature,
@@ -992,93 +958,7 @@ doAllocBuffers(const unsigned bufferSize,
     rstKeywords["PRESSURE"] = 0;
 
     if (enableMech_ && eclState_.runspec().mech()) {
-        this->mechPotentialForce_.resize(bufferSize,0.0);
+        this->mech_.allocate(bufferSize, rstKeywords);
-        rstKeywords["MECHPOTF"] = 0;
-        this->mechPotentialTempForce_.resize(bufferSize,0.0);
-        rstKeywords["TEMPPOTF"] = 0;
-        this->mechPotentialPressForce_.resize(bufferSize,0.0);
-        rstKeywords["PRESPOTF"] = 0;
-
-        this->dispX_.resize(bufferSize,0.0);
-        rstKeywords["DISPX"] = 0;
-        this->dispY_.resize(bufferSize,0.0);
-        rstKeywords["DISPY"] = 0;
-        this->dispZ_.resize(bufferSize,0.0);
-        rstKeywords["DISPZ"] = 0;
-        this->stressXX_.resize(bufferSize,0.0);
-        rstKeywords["STRESSXX"] = 0;
-        this->stressYY_.resize(bufferSize,0.0);
-        rstKeywords["STRESSYY"] = 0;
-        this->stressZZ_.resize(bufferSize,0.0);
-        rstKeywords["STRESSZZ"] = 0;
-        this->stressXY_.resize(bufferSize,0.0);
-        rstKeywords["STRESSXY"] = 0;
-        this->stressXZ_.resize(bufferSize,0.0);
-        rstKeywords["STRESSXZ"] = 0;
-        this->stressXY_.resize(bufferSize,0.0);
-        rstKeywords["STRESSXY"] = 0;
-        this->stressYZ_.resize(bufferSize,0.0);
-        rstKeywords["STRESSYZ"] = 0;
-
-        this->strainXX_.resize(bufferSize,0.0);
-        rstKeywords["STRAINXX"] = 0;
-        this->strainYY_.resize(bufferSize,0.0);
-        rstKeywords["STRAINYY"] = 0;
-        this->strainZZ_.resize(bufferSize,0.0);
-        rstKeywords["STRAINZZ"] = 0;
-        this->strainXY_.resize(bufferSize,0.0);
-        rstKeywords["STRAINXY"] = 0;
-        this->strainXZ_.resize(bufferSize,0.0);
-        rstKeywords["STRAINXZ"] = 0;
-        this->strainXY_.resize(bufferSize,0.0);
-        rstKeywords["STRAINXY"] = 0;
-        this->strainYZ_.resize(bufferSize,0.0);
-        rstKeywords["STRAINYZ"] = 0;
-
-        this->delstressXX_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRXX"] = 0;
-        this->delstressYY_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRYY"] = 0;
-        this->delstressZZ_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRZZ"] = 0;
-        this->delstressXY_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRXY"] = 0;
-        this->delstressXZ_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRXZ"] = 0;
-        this->delstressXY_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRXY"] = 0;
-        this->delstressYZ_.resize(bufferSize,0.0);
-        rstKeywords["DELSTRYZ"] = 0;
-
-        this->fracstressXX_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRXX"] = 0;
-        this->fracstressYY_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRYY"] = 0;
-        this->fracstressZZ_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRZZ"] = 0;
-        this->fracstressXY_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRXY"] = 0;
-        this->fracstressXZ_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRXZ"] = 0;
-        this->fracstressXY_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRXY"] = 0;
-        this->fracstressYZ_.resize(bufferSize,0.0);
-        rstKeywords["FRCSTRYZ"] = 0;
-
-        this->linstressXX_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRXX"] = 0;
-        this->linstressYY_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRYY"] = 0;
-        this->linstressZZ_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRZZ"] = 0;
-        this->linstressXY_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRXY"] = 0;
-        this->linstressXZ_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRXZ"] = 0;
-        this->linstressXY_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRXY"] = 0;
-        this->linstressYZ_.resize(bufferSize,0.0);
-        rstKeywords["LINSTRYZ"] = 0;
     }
 
     // If TEMP is set in RPTRST we output temperature even if THERMAL

opm/simulators/flow/GenericOutputBlackoilModule.hpp

@@ -36,6 +36,7 @@
 #include <opm/simulators/flow/FlowsData.hpp>
 #include <opm/simulators/flow/InterRegFlows.hpp>
 #include <opm/simulators/flow/LogOutputHelper.hpp>
+#include <opm/simulators/flow/MechContainer.hpp>
 #include <opm/simulators/flow/RegionPhasePVAverage.hpp>
 
 #include <opm/simulators/utils/ParallelCommunication.hpp>
@@ -503,43 +504,7 @@ protected:
     ScalarBuffer pcog_;
 
     // buffers for mechanical output
-    ScalarBuffer mechPotentialForce_;
+    MechContainer<Scalar> mech_;
-    ScalarBuffer mechPotentialPressForce_;
-    ScalarBuffer mechPotentialTempForce_;
-
-    ScalarBuffer dispX_;
-    ScalarBuffer dispY_;
-    ScalarBuffer dispZ_;
-    ScalarBuffer stressXX_;
-    ScalarBuffer stressYY_;
-    ScalarBuffer stressZZ_;
-    ScalarBuffer stressXY_;
-    ScalarBuffer stressXZ_;
-    ScalarBuffer stressYZ_;
-    ScalarBuffer delstressXX_;
-    ScalarBuffer delstressYY_;
-    ScalarBuffer delstressZZ_;
-    ScalarBuffer delstressXY_;
-    ScalarBuffer delstressXZ_;
-    ScalarBuffer delstressYZ_;
-    ScalarBuffer linstressXX_;
-    ScalarBuffer linstressYY_;
-    ScalarBuffer linstressZZ_;
-    ScalarBuffer linstressXY_;
-    ScalarBuffer linstressXZ_;
-    ScalarBuffer linstressYZ_;
-    ScalarBuffer fracstressXX_;
-    ScalarBuffer fracstressYY_;
-    ScalarBuffer fracstressZZ_;
-    ScalarBuffer fracstressXY_;
-    ScalarBuffer fracstressXZ_;
-    ScalarBuffer fracstressYZ_;
-    ScalarBuffer strainXX_;
-    ScalarBuffer strainYY_;
-    ScalarBuffer strainZZ_;
-    ScalarBuffer strainXY_;
-    ScalarBuffer strainXZ_;
-    ScalarBuffer strainYZ_;
 
     std::array<ScalarBuffer, numPhases> saturation_;
     std::array<ScalarBuffer, numPhases> invB_;

opm/simulators/flow/MechContainer.cpp

@@ -0,0 +1,214 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+
+#include <config.h>
+#include <opm/simulators/flow/MechContainer.hpp>
+
+#include <opm/common/utility/Visitor.hpp>
+
+#include <opm/output/data/Solution.hpp>
+
+#include <algorithm>
+#include <array>
+#include <tuple>
+#include <variant>
+
+namespace Opm {
+
+template<class Scalar>
+void MechContainer<Scalar>::
+allocate(const std::size_t bufferSize,
+         std::map<std::string, int>& rstKeywords)
+{
+    this->potentialForce_.resize(bufferSize, 0.0);
+    rstKeywords["MECHPOTF"] = 0;
+    this->potentialTempForce_.resize(bufferSize, 0.0);
+    rstKeywords["TEMPPOTF"] = 0;
+    this->potentialPressForce_.resize(bufferSize, 0.0);
+    rstKeywords["PRESPOTF"] = 0;
+
+    std::for_each(disp_.begin(), disp_.end(),
+                  [suffix = "XYZ", is = 0, bufferSize, &rstKeywords](auto& v) mutable
+                  {
+                      v.resize(bufferSize, 0.0);
+                      rstKeywords[std::string{"DISP"} + suffix[is++]] = 0;
+                  });
+
+    auto resizeAndRegister =
+        [&rstKeywords,bufferSize](auto& tensor,
+                                  const std::string& name)
+    {
+        static constexpr auto fields = std::array{
+            "XX", "YY", "ZZ", "YZ", "XZ", "XY",
+        };
+        tensor.resize(bufferSize);
+        for (const auto& f : fields) {
+            rstKeywords[name + f] = 0;
+        }
+    };
+
+    resizeAndRegister(delstress_, "DELSTR");
+    resizeAndRegister(fracstress_, "FRCSTR");
+    resizeAndRegister(linstress_, "LINSTR");
+    resizeAndRegister(strain_, "STRAIN");
+    resizeAndRegister(stress_, "STRESS");
+
+    allocated_ = true;
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignDelStress(const unsigned globalDofIdx,
+                const Dune::FieldVector<Scalar,6>& delStress)
+{
+    this->delstress_.assign(globalDofIdx, VoigtContainer<Scalar>(delStress));
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignDisplacement(const unsigned globalDofIdx,
+                   const Dune::FieldVector<Scalar,3>& disp)
+{
+    this->disp_[0][globalDofIdx] = disp[0];
+    this->disp_[1][globalDofIdx] = disp[1];
+    this->disp_[2][globalDofIdx] = disp[2];
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignFracStress(const unsigned globalDofIdx,
+                 const Dune::FieldVector<Scalar,6>& fracStress)
+{
+    this->fracstress_.assign(globalDofIdx, VoigtContainer<Scalar>(fracStress));
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignLinStress(const unsigned globalDofIdx,
+                const Dune::FieldVector<Scalar,6>& linStress)
+{
+    this->linstress_.assign(globalDofIdx, VoigtContainer<Scalar>(linStress));
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignPotentialForces(const unsigned globalDofIdx,
+                      const Scalar force,
+                      const Scalar pressForce,
+                      const Scalar tempForce)
+{
+    potentialForce_[globalDofIdx] = force;
+    potentialPressForce_[globalDofIdx] = pressForce;
+    potentialTempForce_[globalDofIdx] = tempForce;
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignStrain(const unsigned globalDofIdx,
+             const Dune::FieldVector<Scalar,6>& strain)
+{
+    this->strain_.assign(globalDofIdx, VoigtContainer<Scalar>(strain));
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+assignStress(const unsigned globalDofIdx,
+             const Dune::FieldVector<Scalar,6>& stress)
+{
+    this->stress_.assign(globalDofIdx, VoigtContainer<Scalar>(stress));
+}
+
+template<class Scalar>
+void MechContainer<Scalar>::
+outputRestart(data::Solution& sol) const
+{
+    if (!allocated_) {
+        return;
+    }
+    using DataEntry = std::tuple<std::string,
+                                 UnitSystem::measure,
+                                 std::variant<const std::vector<Scalar>*,
+                                              const std::array<std::vector<Scalar>,3>*,
+                                              const VoigtArray<Scalar>*>>;
+
+    auto doInsert = [&sol](const std::string& name,
+                           const UnitSystem::measure& measure,
+                           const std::vector<Scalar>& entry)
+    {
+        if (!entry.empty()) {
+            sol.insert(name, measure, entry, data::TargetType::RESTART_OPM_EXTENDED);
+        }
+    };
+
+    const auto solutionVectors = std::array{
+        DataEntry{"DELSTR",   UnitSystem::measure::pressure, &delstress_},
+        DataEntry{"DISP",     UnitSystem::measure::length,   &disp_},
+        DataEntry{"FRCSTR",   UnitSystem::measure::pressure, &fracstress_},
+        DataEntry{"LINSTR",   UnitSystem::measure::pressure, &linstress_},
+        DataEntry{"MECHPOTF", UnitSystem::measure::pressure, &potentialForce_},
+        DataEntry{"PRESPOTF", UnitSystem::measure::pressure, &potentialPressForce_},
+        DataEntry{"STRAIN",   UnitSystem::measure::identity, &strain_},
+        DataEntry{"STRESS",   UnitSystem::measure::length,   &stress_},
+        DataEntry{"TEMPPOTF", UnitSystem::measure::pressure, &potentialTempForce_},
+    };
+
+    std::for_each(solutionVectors.begin(), solutionVectors.end(),
+                  [&doInsert](const auto& array)
+                  {
+                      std::visit(VisitorOverloadSet{
+                                    [&array, &doInsert](const std::vector<Scalar>* v)
+                                    {
+                                        doInsert(std::get<0>(array), std::get<1>(array), *v);
+                                    },
+                                    [&array, &doInsert](const std::array<std::vector<Scalar>,3>* V)
+                                    {
+                                        const auto& v = *V;
+                                        const auto& name = std::get<0>(array);
+                                        const auto& measure = std::get<1>(array);
+                                        doInsert(name + "X", measure, v[0]);
+                                        doInsert(name + "Y", measure, v[1]);
+                                        doInsert(name + "Z", measure, v[2]);
+                                    },
+                                    [&array, &doInsert](const VoigtArray<Scalar>* V)
+                                    {
+                                        const auto& v = *V;
+                                        const auto& name = std::get<0>(array);
+                                        const auto& measure = std::get<1>(array);
+                                        doInsert(name + "XX", measure, v[VoigtIndex::XX]);
+                                        doInsert(name + "YY", measure, v[VoigtIndex::YY]);
+                                        doInsert(name + "ZZ", measure, v[VoigtIndex::ZZ]);
+                                        doInsert(name + "YZ", measure, v[VoigtIndex::YZ]);
+                                        doInsert(name + "XZ", measure, v[VoigtIndex::XZ]);
+                                        doInsert(name + "XY", measure, v[VoigtIndex::XY]);
+                                    }
+                                 }, std::get<2>(array));
+                  }
+    );
+}
+
+template class MechContainer<double>;
+
+#if FLOW_INSTANTIATE_FLOAT
+template class MechContainer<float>;
+#endif
+
+} // namespace Opm

opm/simulators/flow/MechContainer.hpp

@@ -0,0 +1,96 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ * \copydoc Opm::OutputBlackOilModule
+ */
+#ifndef OPM_MECH_CONTAINER_HPP
+#define OPM_MECH_CONTAINER_HPP
+
+#include <dune/common/fvector.hh>
+
+#include <opm/simulators/utils/VoigtArray.hpp>
+
+#include <array>
+#include <cstddef>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace Opm {
+
+namespace data { class Solution; }
+
+template<class Scalar>
+class MechContainer
+{
+    using ScalarBuffer = std::vector<Scalar>;
+
+public:
+    void allocate(const std::size_t bufferSize,
+                  std::map<std::string, int>& rstKeywords);
+
+    void assignDisplacement(const unsigned globalDofIdx,
+                            const Dune::FieldVector<Scalar,3>& disp);
+
+    void assignDelStress(const unsigned globalDofIdx,
+                         const Dune::FieldVector<Scalar,6>& delStress);
+
+    void assignPotentialForces(const unsigned globalDofIdx,
+                               const Scalar force,
+                               const Scalar pressForce,
+                               const Scalar tempForce);
+
+    void assignFracStress(const unsigned globalDofIdx,
+                          const Dune::FieldVector<Scalar,6>& fracStress);
+
+    void assignLinStress(const unsigned globalDofIdx,
+                         const Dune::FieldVector<Scalar,6>& linStress);
+
+    void assignStrain(const unsigned globalDofIdx,
+                      const Dune::FieldVector<Scalar,6>& strain);
+
+    void assignStress(const unsigned globalDofIdx,
+                      const Dune::FieldVector<Scalar,6>& stress);
+
+    void outputRestart(data::Solution& sol) const;
+
+    bool allocated() const
+    { return allocated_; }
+
+private:
+    bool allocated_ = false;
+    ScalarBuffer potentialForce_;
+    ScalarBuffer potentialPressForce_;
+    ScalarBuffer potentialTempForce_;
+
+    std::array<ScalarBuffer,3> disp_;
+    VoigtArray<Scalar> delstress_;
+    VoigtArray<Scalar> fracstress_;
+    VoigtArray<Scalar> linstress_;
+    VoigtArray<Scalar> strain_;
+    VoigtArray<Scalar> stress_;
+};
+
+} // namespace Opm
+
+#endif // OPM_MECH_CONTAINER_HPP

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -196,12 +196,10 @@ public:
     void processElementMech(const ElementContext& elemCtx)
     {
         if constexpr (getPropValue<TypeTag, Properties::EnableMech>()) {
-            if (this->mechPotentialForce_.empty()) {
+            if (!this->mech_.allocated()) {
                 return;
             }
 
-            enum Voigt { XX = 0, YY = 1, ZZ = 2, YZ = 3, XZ = 4, XY = 5, };
-
             const auto& problem = elemCtx.simulator().problem();
             const auto& model = problem.geoMechModel();
 
@@ -209,57 +207,30 @@ public:
                 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
 
                 // Assume all mechanical things should be written
-                this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx);
+                this->mech_.assignDelStress(globalDofIdx,
-                this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx);
+                                            model.delstress(globalDofIdx));
-                this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx);
 
-                const auto disp = model.disp(globalDofIdx, /*include_fracture*/true);
+                this->mech_.assignDisplacement(globalDofIdx,
-                this->dispX_[globalDofIdx] = disp[Voigt::XX];
+                                               model.disp(globalDofIdx, /*include_fracture*/true));
-                this->dispY_[globalDofIdx] = disp[Voigt::YY];
-                this->dispZ_[globalDofIdx] = disp[Voigt::ZZ];
-
-                // Total stress is not stored but calculated result is Voigt notation
-                const auto stress = model.stress(globalDofIdx, /*include_fracture*/true);
-                this->stressXX_[globalDofIdx] = stress[Voigt::XX];
-                this->stressYY_[globalDofIdx] = stress[Voigt::YY];
-                this->stressZZ_[globalDofIdx] = stress[Voigt::ZZ];
-                this->stressXY_[globalDofIdx] = stress[Voigt::XY];
-                this->stressXZ_[globalDofIdx] = stress[Voigt::XZ];
-                this->stressYZ_[globalDofIdx] = stress[Voigt::YZ];
-
-                const auto strain = model.strain(globalDofIdx, /*include_fracture*/true);
-                this->strainXX_[globalDofIdx] = strain[Voigt::XX];
-                this->strainYY_[globalDofIdx] = strain[Voigt::YY];
-                this->strainZZ_[globalDofIdx] = strain[Voigt::ZZ];
-                this->strainXY_[globalDofIdx] = strain[Voigt::XY];
-                this->strainXZ_[globalDofIdx] = strain[Voigt::XZ];
-                this->strainYZ_[globalDofIdx] = strain[Voigt::YZ];
-
-                // Not including fracture
-                const auto delstress = model.delstress(globalDofIdx);
-                this->delstressXX_[globalDofIdx] = delstress[Voigt::XX];
-                this->delstressYY_[globalDofIdx] = delstress[Voigt::YY];
-                this->delstressZZ_[globalDofIdx] = delstress[Voigt::ZZ];
-                this->delstressXY_[globalDofIdx] = delstress[Voigt::XY];
-                this->delstressXZ_[globalDofIdx] = delstress[Voigt::XZ];
-                this->delstressYZ_[globalDofIdx] = delstress[Voigt::YZ];
-
-                const auto linstress = model.linstress(globalDofIdx);
-                this->linstressXX_[globalDofIdx] = linstress[Voigt::XX];
-                this->linstressYY_[globalDofIdx] = linstress[Voigt::YY];
-                this->linstressZZ_[globalDofIdx] = linstress[Voigt::ZZ];
-                this->linstressXY_[globalDofIdx] = linstress[Voigt::XY];
-                this->linstressXZ_[globalDofIdx] = linstress[Voigt::XZ];
-                this->linstressYZ_[globalDofIdx] = linstress[Voigt::YZ];
 
                 // is the tresagii stress which make rock fracture
-                const auto fracstress = model.fractureStress(globalDofIdx);
+                this->mech_.assignFracStress(globalDofIdx,
-                this->fracstressXX_[globalDofIdx] = fracstress[Voigt::XX];
+                                             model.fractureStress(globalDofIdx));
-                this->fracstressYY_[globalDofIdx] = fracstress[Voigt::YY];
+
-                this->fracstressZZ_[globalDofIdx] = fracstress[Voigt::ZZ];
+                this->mech_.assignLinStress(globalDofIdx,
-                this->fracstressXY_[globalDofIdx] = fracstress[Voigt::XY];
+                                            model.linstress(globalDofIdx));
-                this->fracstressXZ_[globalDofIdx] = fracstress[Voigt::XZ];
+
-                this->fracstressYZ_[globalDofIdx] = fracstress[Voigt::YZ];
+                this->mech_.assignPotentialForces(globalDofIdx,
+                                                  model.mechPotentialForce(globalDofIdx),
+                                                  model.mechPotentialPressForce(globalDofIdx),
+                                                  model.mechPotentialTempForce(globalDofIdx));
+
+                this->mech_.assignStrain(globalDofIdx,
+                                         model.strain(globalDofIdx, /*include_fracture*/true));
+
+                // Total stress is not stored but calculated result is Voigt notation
+                this->mech_.assignStress(globalDofIdx,
+                                         model.stress(globalDofIdx, /*include_fracture*/true));;
             }
         }
     }