diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh index 2b69a04218f76320fcd30b894887bf6c67d469cc..bc59b27fda3b24ad71f4d44bd50ada830319bde5 100644 --- a/dumux/geomechanics/elastic/model.hh +++ b/dumux/geomechanics/elastic/model.hh @@ -26,6 +26,8 @@ #ifndef DUMUX_GEOMECHANICS_ELASTIC_MODEL_HH #define DUMUX_GEOMECHANICS_ELASTIC_MODEL_HH +#include <dune/common/fvector.hh> + #include <dumux/common/properties.hh> #include <dumux/common/properties/model.hh> @@ -66,14 +68,16 @@ struct ElasticModelTraits * \brief Traits class for the volume variables of the elastic model. * * \tparam PV The type used for primary variables + * \tparam DV The type used for displacement vectors * \tparam MT The model traits * \tparam SST The solid state * \tparam SSY The solid system */ -template<class PV, class MT, class SST, class SSY> +template<class PV, class DV, class MT, class SST, class SSY> struct ElasticVolumeVariablesTraits { using PrimaryVariables = PV; + using DisplacementVector = DV; using ModelTraits = MT; using SolidState = SST; using SolidSystem = SSY; @@ -95,11 +99,13 @@ SET_TYPE_PROP(Elastic, ModelTraits, ElasticModelTraits< GET_PROP_TYPE(TypeTag, G SET_PROP(Elastic, VolumeVariables) { private: + static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using DV = Dune::FieldVector<typename PV::value_type, dim>; using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits); using SST = typename GET_PROP_TYPE(TypeTag, SolidState); using SSY = typename GET_PROP_TYPE(TypeTag, SolidSystem); - using Traits = ElasticVolumeVariablesTraits<PV, MT, SST, SSY>; + using Traits = ElasticVolumeVariablesTraits<PV, DV, MT, SST, SSY>; public: using type = ElasticVolumeVariables<Traits>; }; diff --git a/dumux/geomechanics/elastic/volumevariables.hh b/dumux/geomechanics/elastic/volumevariables.hh index 9e805a0b3470884742a639ebb1b52a05896b3ac2..038d87c5e97fb5b311b5f6006e566bb879f43371 100644 --- a/dumux/geomechanics/elastic/volumevariables.hh +++ b/dumux/geomechanics/elastic/volumevariables.hh @@ -52,6 +52,8 @@ class ElasticVolumeVariables public: //! export the type used for the primary variables using PrimaryVariables = typename Traits::PrimaryVariables; + //! export the type used for displacement vectors + using DisplacementVector = typename Traits::DisplacementVector; //! export the type encapsulating primary variable indices using Indices = typename ModelTraits::Indices; //! export type of solid state @@ -89,10 +91,19 @@ public: Scalar solidDensity() const { return solidState_.density(); } - //! Returns the permeability within the control volume in \f$[m^2]\f$. + //! Returns the permeability within the control volume in \f$[m]\f$. Scalar displacement(unsigned int dir) const { return priVars_[ Indices::momentum(dir) ]; } + //! Returns the displacement vector within the scv in \f$[m]\f$. + DisplacementVector displacement() const + { + DisplacementVector d; + for (int dir = 0; dir < d.size(); ++dir) + d[dir] = displacement(dir); + return d; + } + //! Return a component of primary variable vector for a given index Scalar priVar(const int pvIdx) const { return priVars_[pvIdx]; } diff --git a/dumux/geomechanics/poroelastic/model.hh b/dumux/geomechanics/poroelastic/model.hh index 370b8629b59f2f15efedbef69eeb35c766b22e6b..4ea6e0494ffef92aad7cc1505d813fdf6735d089 100644 --- a/dumux/geomechanics/poroelastic/model.hh +++ b/dumux/geomechanics/poroelastic/model.hh @@ -26,6 +26,8 @@ #ifndef DUMUX_GEOMECHANICS_POROELASTIC_MODEL_HH #define DUMUX_GEOMECHANICS_POROELASTIC_MODEL_HH +#include <dune/common/fvector.hh> + #include <dumux/common/properties.hh> #include <dumux/geomechanics/elastic/indices.hh> #include <dumux/geomechanics/elastic/model.hh> @@ -90,13 +92,15 @@ public: SET_PROP(PoroElastic, VolumeVariables) { private: + static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension; using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using DV = Dune::FieldVector<typename PV::value_type, dim>; using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits); using SST = typename GET_PROP_TYPE(TypeTag, SolidState); using SSY = typename GET_PROP_TYPE(TypeTag, SolidSystem); // we reuse the elastic volume variable traits here - using Traits = ElasticVolumeVariablesTraits<PV, MT, SST, SSY>; + using Traits = ElasticVolumeVariablesTraits<PV, DV, MT, SST, SSY>; public: using type = PoroElasticVolumeVariables<Traits>; }; diff --git a/dumux/geomechanics/poroelastic/volumevariables.hh b/dumux/geomechanics/poroelastic/volumevariables.hh index c9c5d504f92b1f55bf8bcbf787d9ed98e50075cf..cd285e7501561c78b991fd47a2cf2baedaacdf47 100644 --- a/dumux/geomechanics/poroelastic/volumevariables.hh +++ b/dumux/geomechanics/poroelastic/volumevariables.hh @@ -46,6 +46,8 @@ class PoroElasticVolumeVariables public: //! export the type used for the primary variables using PrimaryVariables = typename Traits::PrimaryVariables; + //! export the type used for displacement vectors + using DisplacementVector = typename Traits::DisplacementVector; //! export the type encapsulating primary variable indices using Indices = typename ModelTraits::Indices; //! export type of solid state @@ -97,10 +99,19 @@ public: Scalar divU() const { return divU_; } - //! Returns the permeability within the scv in \f$[m^2]\f$. + //! Returns the permeability within the scv in \f$[m]\f$. Scalar displacement(unsigned int dir) const { return priVars_[ Indices::momentum(dir) ]; } + //! Returns the displacement vector within the scv in \f$[m]\f$. + DisplacementVector displacement() const + { + DisplacementVector d; + for (int dir = 0; dir < d.size(); ++dir) + d[dir] = displacement(dir); + return d; + } + //! Return a component of primary variable vector for a given index Scalar priVar(const int pvIdx) const { return priVars_[pvIdx]; } diff --git a/dumux/geomechanics/poroelastic/vtkoutputfields.hh b/dumux/geomechanics/poroelastic/vtkoutputfields.hh index 21d1f9a36fa805cadf6f95c8e4c51c910bf83366..4143ce3b2c163bded5c2e334454f450d7695f336 100644 --- a/dumux/geomechanics/poroelastic/vtkoutputfields.hh +++ b/dumux/geomechanics/poroelastic/vtkoutputfields.hh @@ -36,6 +36,7 @@ public: template <class VtkOutputModule> static void init(VtkOutputModule& vtk) { + vtk.addVolumeVariable([](const auto& volVars){ return volVars.displacement(); }, "u"); vtk.addVolumeVariable([](const auto& volVars){ return volVars.divU(); }, "divU"); vtk.addVolumeVariable([](const auto& volVars){ return volVars.porosity(); }, "porosity"); } diff --git a/test/geomechanics/poroelastic/test_poroelastic.cc b/test/geomechanics/poroelastic/test_poroelastic.cc index 61acaa7dd76f80067305b64e0b77a9cfbc7bbae8..dba038fe79d08e4335790a5630071d8c2405989c 100644 --- a/test/geomechanics/poroelastic/test_poroelastic.cc +++ b/test/geomechanics/poroelastic/test_poroelastic.cc @@ -149,8 +149,7 @@ int main(int argc, char** argv) try VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); VtkOutputFields::init(vtkWriter); - // also, add displacement and exact solution to the output - vtkWriter.addField(x, "u"); + // also, add exact solution to the output SolutionVector xExact(fvGridGeometry->numDofs()); for (const auto& v : vertices(leafGridView)) xExact[ fvGridGeometry->vertexMapper().index(v) ] = problem->exactSolution(v.geometry().center());