diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
index 696feee855797e07d660154970618c81379ec199..61918c6fcfcd81959d5039aba1e750ef7553645d 100644
--- a/dumux/CMakeLists.txt
+++ b/dumux/CMakeLists.txt
@@ -2,6 +2,7 @@ add_subdirectory("adaptive")
add_subdirectory("assembly")
add_subdirectory("common")
add_subdirectory("discretization")
+add_subdirectory("geomechanics")
add_subdirectory("freeflow")
add_subdirectory("io")
add_subdirectory("linear")
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 78a0ac7c6a2757a28e150e5ca664fadf05ac78a4..08a752fc93c31bac0b28ce476f279a14eafbb1a0 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -138,6 +138,11 @@ NEW_PROP_TAG(EnableWaterDiffusionInAir); //!< Property for turning Richards into
//////////////////////////////////////////////////////////////
NEW_PROP_TAG(OnlyGasPhaseCanDisappear); //!< reduces the phasestates to threePhases and wnPhaseOnly
+/////////////////////////////////////////////////////////////
+// Properties used by geomechanical models:
+/////////////////////////////////////////////////////////////
+NEW_PROP_TAG(StressType); //!< The type used for the evaluation of stress tensors and forces
+
/////////////////////////////////////////////////////////////
// Properties used by the staggered-grid discretization method
/////////////////////////////////////////////////////////////
diff --git a/dumux/discretization/box/fluxvariablescache.hh b/dumux/discretization/box/fluxvariablescache.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c69fb4c9d22834eb3e4ce279708d308afbfaa64c
--- /dev/null
+++ b/dumux/discretization/box/fluxvariablescache.hh
@@ -0,0 +1,101 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Discretization
+ * \brief Flux variables cache class for the box scheme
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_FLUXVARIABLES_CACHE_HH
+#define DUMUX_DISCRETIZATION_BOX_FLUXVARIABLES_CACHE_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Discretization
+ * \ingroup BoxDiscretization
+ * \brief Flux variables cache class for the box scheme.
+ * For the box scheme, this class does not contain any physics-/process-dependent
+ * data. It solely stores disretization-/grid-related data.
+ */
+template< class Scalar, class FVGridGeometry >
+class BoxFluxVariablesCache
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using CoordScalar = typename GridView::ctype;
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
+ using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
+ using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
+
+public:
+
+ template< class Problem, class ElementVolumeVariables >
+ void update(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace &scvf)
+ {
+ const auto geometry = element.geometry();
+ const auto& localBasis = fvGeometry.feLocalBasis();
+
+ // evaluate shape functions and gradients at the integration point
+ const auto ipLocal = geometry.local(scvf.ipGlobal());
+ jacInvT_ = geometry.jacobianInverseTransposed(ipLocal);
+ localBasis.evaluateJacobian(ipLocal, shapeJacobian_);
+ localBasis.evaluateFunction(ipLocal, shapeValues_); // shape values for rho
+
+ // compute the gradN at for every scv/dof
+ gradN_.resize(fvGeometry.numScv());
+ for (const auto& scv: scvs(fvGeometry))
+ jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.indexInElement()]);
+ }
+
+ //! returns the shape function gradients in local coordinates at the scvf integration point
+ const std::vector<ShapeJacobian>& shapeJacobian() const { return shapeJacobian_; }
+ //! returns the shape function values at the scvf integration point
+ const std::vector<ShapeValue>& shapeValues() const { return shapeValues_; }
+ //! returns inverse transposed jacobian at the scvf integration point
+ const JacobianInverseTransposed& jacInvT() const { return jacInvT_; }
+ //! returns the shape function gradients in global coordinates at the scvf integration point
+ const GlobalPosition& gradN(unsigned int scvIdxInElement) const { return gradN_[scvIdxInElement]; }
+
+private:
+ std::vector<GlobalPosition> gradN_;
+ std::vector<ShapeJacobian> shapeJacobian_;
+ std::vector<ShapeValue> shapeValues_;
+ JacobianInverseTransposed jacInvT_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/box/hookeslaw.hh b/dumux/discretization/box/hookeslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0749d4b05eb2cce2da6a73e2ae5e7faa81ee46da
--- /dev/null
+++ b/dumux/discretization/box/hookeslaw.hh
@@ -0,0 +1,126 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Specialization of Hooke's law for the box scheme. This computes
+ * the stress tensor and surface forces resulting from mechanical deformation.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_HOOKES_LAW_HH
+#define DUMUX_DISCRETIZATION_BOX_HOOKES_LAW_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoxDiscretization
+ * \brief Hooke's law for box scheme
+ * \tparam Scalar the scalar type for scalar physical quantities
+ * \tparam FVGridGeometry the grid geometry
+ */
+template<class Scalar, class FVGridGeometry>
+class HookesLaw<Scalar, FVGridGeometry, DiscretizationMethod::box>
+{
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using CoordScalar = typename GridView::ctype;
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ static_assert(dim == dimWorld, "Hookes Law not implemented for network/surface grids");
+
+public:
+ //! export the type used for the stress tensor
+ using StressTensor = Dune::FieldMatrix<Scalar, dim, dimWorld>;
+ //! export the type used for force vectors
+ using ForceVector = typename StressTensor::row_type;
+
+ //! computes the force acting on a sub-control volume face
+ template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
+ static ForceVector force(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVarsCache& elemFluxVarCache)
+ {
+ const auto sigma = stressTensor(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarCache);
+
+ ForceVector scvfForce(0.0);
+ sigma.mv(scvf.unitOuterNormal(), scvfForce);
+ scvfForce *= scvf.area();
+
+ return scvfForce;
+ }
+
+ //! assembles the stress tensor at the integration point of an scvf
+ template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
+ static StressTensor stressTensor(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVarsCache& elemFluxVarCache)
+ {
+ const auto& fluxVarCache = elemFluxVarCache[scvf];
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+
+ // TODO What to do with extrusion factor??
+ // TODO this only works for element-wise constant lame params
+ // How could this be done for heterogeneous params within the element?
+ const auto& lameParams = insideVolVars.lameParams();
+
+ // evaluate displacement gradient
+ StressTensor gradU(0.0);
+ for (int dir = 0; dir < dim; ++dir)
+ for (const auto& scv : scvs(fvGeometry))
+ gradU[dir].axpy(elemVolVars[scv.indexInElement()].displacement(dir), fluxVarCache.gradN(scv.indexInElement()));
+
+ // evaluate strain tensor
+ StressTensor epsilon;
+ for (int i = 0; i < dim; ++i)
+ for (int j = 0; j < dimWorld; ++j)
+ epsilon[i][j] = 0.5*(gradU[i][j] + gradU[j][i]);
+
+ // calculate sigma
+ StressTensor sigma(0.0);
+ const auto traceEpsilon = trace(epsilon);
+ for (int i = 0; i < dim; ++i)
+ {
+ sigma[i][i] = lameParams.lambda()*traceEpsilon;
+ for (int j = 0; j < dimWorld; ++j)
+ sigma[i][j] += 2.0*lameParams.mu()*epsilon[i][j];
+ }
+
+ return sigma;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/hookeslaw.hh b/dumux/discretization/hookeslaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f3abc81ec7593fad57cfa2791e10969d3f5d9e19
--- /dev/null
+++ b/dumux/discretization/hookeslaw.hh
@@ -0,0 +1,45 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Discretization
+ * \brief Hooke's law specialized for different discretization schemes.
+ * This computes the stress tensor and surface forces resulting from mechanical deformation.
+ */
+#ifndef DUMUX_DISCRETIZATION_HOOKES_LAW_HH
+#define DUMUX_DISCRETIZATION_HOOKES_LAW_HH
+
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Discretization
+ * \brief Evaluates the normal component of the Darcy velocity on a (sub)control volume face.
+ * \note Specializations are provided for the different discretization methods.
+ * These specializations are found in the headers included below.
+ */
+template <class Scalar, class FVGridGeometry, DiscretizationMethod dm = FVGridGeometry::discMethod>
+class HookesLaw;
+
+} // end namespace Dumux
+
+#include <dumux/discretization/box/hookeslaw.hh>
+
+#endif
diff --git a/dumux/discretization/staggered/freeflow/velocityoutput.hh b/dumux/discretization/staggered/freeflow/velocityoutput.hh
index 05ed28c28bde9c98f83afc9e84437d9524b62cd2..fadcbd7115836ea63febccab7ee5a1611df33cdf 100644
--- a/dumux/discretization/staggered/freeflow/velocityoutput.hh
+++ b/dumux/discretization/staggered/freeflow/velocityoutput.hh
@@ -78,6 +78,14 @@ public:
bool enableOutput()
{ return velocityOutput_; }
+ // returns the name of the phase for a given index
+ static std::string phaseName(int phaseIdx)
+ { return GET_PROP_TYPE(TypeTag, FluidSystem)::phaseName(phaseIdx); }
+
+ // returns the number of phase velocities computed by this class
+ static constexpr int numPhaseVelocities()
+ { return GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases(); }
+
//! Return the problem boundary types
auto problemBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
{ return problem_.boundaryTypes(element, scvf); }
diff --git a/dumux/geomechanics/CMakeLists.txt b/dumux/geomechanics/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6cf6436721977b020d0603c6f855ca96e38e7a03
--- /dev/null
+++ b/dumux/geomechanics/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_subdirectory("elastic")
+
+install(FILES
+fvproblem.hh
+lameparams.hh
+properties.hh
+stressvariablescache.hh
+velocityoutput.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics)
diff --git a/dumux/geomechanics/elastic/CMakeLists.txt b/dumux/geomechanics/elastic/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d37ff68eeca3061bb5dd71a5dc787b4fa583c77a
--- /dev/null
+++ b/dumux/geomechanics/elastic/CMakeLists.txt
@@ -0,0 +1,6 @@
+install(FILES
+indices.hh
+localresidual.hh
+model.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics/elastic)
diff --git a/dumux/geomechanics/elastic/indices.hh b/dumux/geomechanics/elastic/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..782d06c20e525474e1505d8b4c3619c3069f5723
--- /dev/null
+++ b/dumux/geomechanics/elastic/indices.hh
@@ -0,0 +1,58 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Defines the indices for the elastic model
+ */
+#ifndef DUMUX_ELASTIC_INDICES_HH
+#define DUMUX_ELASTIC_INDICES_HH
+
+namespace Dumux {
+// \{
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief The indices for the linear elasticity model.
+ */
+struct ElasticIndices
+{
+ // returns the equation index for a given space direction
+ static constexpr int momentum(int dirIdx) { return dirIdx; };
+
+ // returns the primary variable index for a given space direction
+ static constexpr int u(int dirIdx) { return dirIdx; };
+
+ // Equation indices
+ static constexpr int momentumXEqIdx = 0;
+ static constexpr int momentumYEqIdx = 1;
+ static constexpr int momentumZEqIdx = 2;
+
+ // primary variable indices
+ static constexpr int uxIdx = 0;
+ static constexpr int uyIdx = 1;
+ static constexpr int uzIdx = 2;
+};
+
+// \}
+} // end namespace
+
+#endif
diff --git a/dumux/geomechanics/elastic/localresidual.hh b/dumux/geomechanics/elastic/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..99f6cf0fc7916b412f7c97b475443cc77b6dae74
--- /dev/null
+++ b/dumux/geomechanics/elastic/localresidual.hh
@@ -0,0 +1,151 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Element-wise calculation of the local residual for problems
+ * using the elastic model considering linear elasticity.
+ */
+#ifndef DUMUX_GEOMECHANICS_ELASTIC_LOCAL_RESIDUAL_HH
+#define DUMUX_GEOMECHANICS_ELASTIC_LOCAL_RESIDUAL_HH
+
+#include <vector>
+#include <dumux/common/properties.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Element-wise calculation of the local residual for problems
+ * using the elastic model considering linear elasticity.
+ */
+template<class TypeTag>
+class ElasticLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
+{
+ using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
+ using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
+
+ // class assembling the stress tensor
+ using StressType = typename GET_PROP_TYPE(TypeTag, StressType);
+
+public:
+ using ParentType::ParentType;
+
+ /*!
+ * \brief For the elastic model the storage term is zero since
+ * we neglect inertia forces.
+ *
+ * \param problem The problem
+ * \param scv The sub control volume
+ * \param volVars The current or previous volVars
+ */
+ NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ return NumEqVector(0.0);
+ }
+
+ /*!
+ * \brief Evaluates the force in all grid directions acting
+ * on a face of a sub-control volume.
+ *
+ * \param problem The problem
+ * \param element The current element.
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars The volume variables of the current element
+ * \param scvf The sub control volume face to compute the flux on
+ * \param elemFluxVarsCache The cache related to flux compuation
+ */
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
+ {
+ // obtain force on the face from stress type
+ return StressType::force(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+ }
+
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param problem The problem to solve
+ * \param element The DUNE Codim<0> entity for which the residual
+ * ought to be calculated
+ * \param fvGeometry The finite-volume geometry of the element
+ * \param elemVolVars The volume variables associated with the element stencil
+ * \param scv The sub-control volume over which we integrate the source term
+ * \note This is the default implementation for geomechanical models adding to
+ * the user defined sources the source stemming from the gravitational acceleration.
+ *
+ */
+ NumEqVector computeSource(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ NumEqVector source(0.0);
+
+ // add contributions from volume flux sources
+ source += problem.source(element, fvGeometry, elemVolVars, scv);
+
+ // add contribution from possible point sources
+ source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
+
+ // maybe add gravitational acceleration
+ static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+ if (gravity)
+ {
+ const auto& g = problem.gravityAtPos(scv.center());
+ for (int dir = 0; dir < GridView::dimensionworld; ++dir)
+ source[Indices::momentum(dir)] += elemVolVars[scv].solidDensity()*g[dir];
+ }
+
+ return source;
+ }
+
+protected:
+ Implementation *asImp_()
+ { return static_cast<Implementation *> (this); }
+
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *> (this); }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5a6dad5b4142bdb1e3eda7702a6274454c97262e
--- /dev/null
+++ b/dumux/geomechanics/elastic/model.hh
@@ -0,0 +1,113 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Properties
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Defines a type tag and some properties for the elastic geomechanical model
+ */
+#ifndef DUMUX_GEOMECHANICS_ELASTIC_MODEL_HH
+#define DUMUX_GEOMECHANICS_ELASTIC_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/properties/model.hh>
+
+#include <dumux/geomechanics/properties.hh>
+#include <dumux/discretization/hookeslaw.hh>
+
+#include "indices.hh"
+#include "localresidual.hh"
+#include "volumevariables.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Specifies a number properties of the elastic model
+ */
+template< int dim, int numSolidComp >
+struct ElasticModelTraits
+{
+ //! export the type encapsulating indices
+ using Indices = ElasticIndices;
+ //! the number of equations is equal to grid dimension
+ static constexpr int numEq() { return dim; }
+ //! This model does not consider fluid phases
+ static constexpr int numFluidPhases() { return 0; }
+ //! This model does not consider fluid phases
+ static constexpr int numFluidComponents() { return 0; }
+ //! We have one solid phase here
+ static constexpr int numSolidComponents() { return numSolidComp; }
+
+ //! Energy balance not yet implemented
+ static constexpr bool enableEnergyBalance() { return false; }
+};
+
+/*!
+ * \ingroup Elastic
+ * \brief Traits class for the volume variables of the elastic model.
+ *
+ * \tparam PV The type used for primary variables
+ * \tparam MT The model traits
+ * \tparam LP The class used for storing the lame parameters
+ * \tparam SST The solid state
+ * \tparam SSY The solid system
+ */
+template<class PV, class MT, class LP, class SST, class SSY>
+struct ElasticVolumeVariablesTraits
+{
+ using PrimaryVariables = PV;
+ using ModelTraits = MT;
+ using LameParams = LP;
+ using SolidState = SST;
+ using SolidSystem = SSY;
+};
+
+namespace Properties {
+
+//! Type tag for the elastic geomechanical model
+NEW_TYPE_TAG(Elastic, INHERITS_FROM(Geomechanics));
+
+//! Use the local residual of the elastic model
+SET_TYPE_PROP(Elastic, LocalResidual, ElasticLocalResidual<TypeTag>);
+
+//! By default, we use hooke's law for stress evaluations
+SET_TYPE_PROP(Elastic, ModelTraits, ElasticModelTraits< GET_PROP_TYPE(TypeTag, GridView)::dimension,
+ GET_PROP_TYPE(TypeTag, SolidSystem)::numComponents >);
+
+//! Set the volume variables property
+SET_PROP(Elastic, VolumeVariables)
+{
+private:
+ using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ using LP = typename GET_PROP_TYPE(TypeTag, SpatialParams)::LameParams;
+ using SST = typename GET_PROP_TYPE(TypeTag, SolidState);
+ using SSY = typename GET_PROP_TYPE(TypeTag, SolidSystem);
+ using Traits = ElasticVolumeVariablesTraits<PV, MT, LP, SST, SSY>;
+public:
+ using type = ElasticVolumeVariables<Traits>;
+};
+
+} // namespace Properties
+} // namespace Dumux
+
+ #endif
diff --git a/dumux/geomechanics/elastic/volumevariables.hh b/dumux/geomechanics/elastic/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..80ef4b02410defa075e7c848fb82e13c40d31fdb
--- /dev/null
+++ b/dumux/geomechanics/elastic/volumevariables.hh
@@ -0,0 +1,121 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Quantities required by the elastic model defined on a sub-control volume.
+ */
+#ifndef DUMUX_ELASTIC_VOLUME_VARIABLES_HH
+#define DUMUX_ELASTIC_VOLUME_VARIABLES_HH
+
+#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
+
+namespace Dumux {
+
+/*!
+* \ingroup Geomechanics
+* \ingroup Elastic
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the elastic model.
+ *
+ * \tparam Traits Class encapsulating types to be used by the vol vars
+ */
+template<class Traits>
+class ElasticVolumeVariables
+{
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+ using LameParams = typename Traits::LameParams;
+ using ModelTraits = typename Traits::ModelTraits;
+
+ //! The elastic model only makes sense with inert solid systems
+ static_assert(Traits::SolidSystem::isInert(), "Elastic model can only be used with inert solid systems");
+public:
+ //! export the type used for the primary variables
+ using PrimaryVariables = typename Traits::PrimaryVariables;
+ //! export the type encapsulating primary variable indices
+ using Indices = typename ModelTraits::Indices;
+ //! export type of solid state
+ using SolidState = typename Traits::SolidState;
+ //! export the solid system used
+ using SolidSystem = typename Traits::SolidSystem;
+
+ /*!
+ * \brief Update all quantities for a given control volume
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to
+ * be simulated
+ * \param element An element which contains part of the control volume
+ * \param scv The sub-control volume
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ void update(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv)
+ {
+ priVars_ = elemSol[scv.localDofIndex()];
+ extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol);
+ lameParams_ = problem.spatialParams().lameParams(element, elemSol);
+
+ //! set the volume fractions of the solid components
+ updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, /*numFluidComps=*/0);
+ // set the temperature of the solid phase
+ setSolidTemperature_(problem, elemSol);
+ // update the density of the solid phase
+ solidState_.setDensity(SolidSystem::density(solidState_));
+ };
+
+ //! Return the average porosity \f$\mathrm{[-]}\f$ within the control volume.
+ Scalar solidDensity() const { return solidState_.density(); }
+ //! Returns the permeability within the control volume in \f$[m^2]\f$.
+ Scalar displacement(unsigned int dir) const { return priVars_[ Indices::momentum(dir) ]; }
+ //! Return a component of primary variable vector for a given index
+ Scalar priVar(const int pvIdx) const { return priVars_[pvIdx]; }
+ //! Return the vector of primary variables
+ const PrimaryVariables& priVars() const { return priVars_; }
+ //! Return the object containing the lame params of this scv
+ const LameParams& lameParams() const { return lameParams_; }
+ //! TODO We don't know yet how to interpret extrusion for mechanics
+ static constexpr Scalar extrusionFactor() { return 1.0; }
+
+private:
+ //! sets the temperature in the solid state for non-isothermal models
+ template< class Problem, class ElemSol,
+ bool EB = ModelTraits::enableEnergyBalance(), std::enable_if_t<EB, int> = 0 >
+ void setSolidTemperature_(const Problem& problem, const ElemSol& elemSol)
+ { DUNE_THROW(Dune::InvalidStateException, "Non-isothermal elastic model."); }
+
+ //! sets the temperature in the solid state for isothermal models
+ template< class Problem, class ElemSol,
+ bool EB = ModelTraits::enableEnergyBalance(), std::enable_if_t<!EB, int> = 0 >
+ void setSolidTemperature_(const Problem& problem, const ElemSol& elemSol)
+ { solidState_.setTemperature(problem.temperature()); }
+
+ // data members
+ Scalar extrusionFactor_;
+ PrimaryVariables priVars_;
+ LameParams lameParams_;
+ SolidState solidState_;
+};
+
+}
+
+#endif
diff --git a/dumux/geomechanics/fvproblem.hh b/dumux/geomechanics/fvproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dd180b1e8e90db75e6dff3332fa04f7ac02f14ea
--- /dev/null
+++ b/dumux/geomechanics/fvproblem.hh
@@ -0,0 +1,49 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for all geomechanical problems
+ */
+#ifndef DUMUX_GEOMECHANICS_FV_PROBLEM_HH
+#define DUMUX_GEOMECHANICS_FV_PROBLEM_HH
+
+#include <dumux/porousmediumflow/problem.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geomechanics
+ * \brief Base class for all geomechanical problems
+ * \note We actually require the same functionality as the
+ * porous medium flow problem, which is why we simply
+ * inherit from that here.
+ */
+template<class TypeTag>
+class GeomechanicsFVProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+public:
+ //! pull up the constructor of the parent class
+ using ParentType::ParentType;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/geomechanics/lameparams.hh b/dumux/geomechanics/lameparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..04b9fea971c11a3c4810891ce4a68d52d42a64d5
--- /dev/null
+++ b/dumux/geomechanics/lameparams.hh
@@ -0,0 +1,57 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief \copydoc Dumux::LameParams
+ */
+#ifndef DUMUX_GEOMECHANICS_LAME_PARAMS_HH
+#define DUMUX_GEOMECHANICS_LAME_PARAMS_HH
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Structure encapsulating the lame parameters
+ */
+template<class Scalar>
+struct LameParams
+{
+ //! Default constructor
+ LameParams() = default;
+ //! Constructor taking lambda and mu directly
+ LameParams(Scalar lambda, Scalar mu) : lambda_(lambda), mu_(mu) {}
+
+ //! Return the first lame parameter
+ Scalar lambda() const { return lambda_; }
+ //! Return the second lame parameter
+ Scalar mu() const { return mu_; }
+
+ //! set the first lame parameter
+ void setLambda(Scalar lambda) { lambda_ = lambda; }
+ //! set the second lame parameter
+ void setMu(Scalar mu) { mu_ = mu; }
+
+private:
+ Scalar lambda_;
+ Scalar mu_;
+};
+}
+#endif
diff --git a/dumux/geomechanics/properties.hh b/dumux/geomechanics/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..34b73e92e9be32512d42878ee4de4f3018629d0e
--- /dev/null
+++ b/dumux/geomechanics/properties.hh
@@ -0,0 +1,78 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Properties
+ * \ingroup Geomechanics
+ * \brief Defines a type tag and some properties for geomechanical DuMuX models.
+ */
+
+#ifndef DUMUX_GEOMECHANICS_PROPERTIES_HH
+#define DUMUX_GEOMECHANICS_PROPERTIES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/properties/model.hh>
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/solidstates/inertsolidstate.hh>
+#include <dumux/material/solidsystems/inertsolidphase.hh>
+#include <dumux/discretization/hookeslaw.hh>
+
+#include "stressvariablescache.hh"
+#include "velocityoutput.hh"
+
+namespace Dumux {
+namespace Properties {
+
+//! Type tag for geomechanical models
+NEW_TYPE_TAG(Geomechanics, INHERITS_FROM(ModelProperties));
+
+//! The flux variables cache class for models involving flow in porous media
+SET_TYPE_PROP(Geomechanics, FluxVariablesCache, StressVariablesCache< typename GET_PROP_TYPE(TypeTag, Scalar),
+ typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >);
+
+//! By default, we use hooke's law for stress evaluations
+SET_TYPE_PROP(Geomechanics, StressType, HookesLaw< typename GET_PROP_TYPE(TypeTag, Scalar),
+ typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >);
+
+//! The (currently empty) velocity output
+SET_TYPE_PROP(Geomechanics, VelocityOutput, GeomechanicsVelocityOutput);
+
+//! The solid state must be inert
+SET_PROP(Geomechanics, SolidState)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using SolidSystem = typename GET_PROP_TYPE(TypeTag, SolidSystem);
+public:
+ using type = InertSolidState<Scalar, SolidSystem>;
+};
+
+//! Per default we use one constant component in the inert solid system
+SET_PROP(Geomechanics, SolidSystem)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using InertComponent = Components::Constant<1, Scalar>;
+public:
+ using type = SolidSystems::InertSolidPhase<Scalar, InertComponent>;
+};
+} // namespace Properties
+} // namespace Dumux
+
+ #endif
diff --git a/dumux/geomechanics/stressvariablescache.hh b/dumux/geomechanics/stressvariablescache.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cd192aa4e7c07b3a08d5fca1433a65bf9e533662
--- /dev/null
+++ b/dumux/geomechanics/stressvariablescache.hh
@@ -0,0 +1,66 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Geomechanics
+ * \brief Base class for the stress variables cache
+ */
+#ifndef DUMUX_GEOMECHANICS_STRESSVARIABLESCACHE_HH
+#define DUMUX_GEOMECHANICS_STRESSVARIABLESCACHE_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/fluxvariablescaching.hh>
+#include <dumux/discretization/box/fluxvariablescache.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geomechanics
+ * \brief The stress variables cache classes for models involving geomechanics.
+ * Store data required for stress calculation.
+ */
+template< class Scalar, class FVGridGeometry, DiscretizationMethod dm = FVGridGeometry::discMethod >
+class StressVariablesCache;
+
+//! We only store discretization-related quantities for the box method.
+template< class Scalar, class FVGridGeometry >
+class StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::box>
+: public BoxFluxVariablesCache< Scalar, FVGridGeometry >
+{};
+
+// specialization for the cell centered tpfa method
+template< class Scalar, class FVGridGeometry >
+class StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::cctpfa> : public FluxVariablesCaching::_EmptyCache
+{
+public:
+ template<typename... Args>
+ void update(Args&&... args) { DUNE_THROW(Dune::NotImplemented, "Geomechanics with cell-centered schemes"); }
+};
+
+// specialization for the cell centered mpfa method
+template< class Scalar, class FVGridGeometry >
+class StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::ccmpfa>
+: public StressVariablesCache<Scalar, FVGridGeometry, DiscretizationMethod::cctpfa>
+{};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/geomechanics/velocityoutput.hh b/dumux/geomechanics/velocityoutput.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dd4571eb49f0ec7f4996c00a1e04785150e53db7
--- /dev/null
+++ b/dumux/geomechanics/velocityoutput.hh
@@ -0,0 +1,60 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Velocity output for geomechanical models
+ */
+#ifndef DUMUX_GEOMECHANICS_VELOCITYOUTPUT_HH
+#define DUMUX_GEOMECHANICS_VELOCITYOUTPUT_HH
+
+namespace Dumux {
+
+/*!
+ * \brief Velocity output for geomechanical models.
+ * This class could be used to compute the temporal derivative
+ * of the displacement. Currently this is not implemented and
+ * we simply define this here in order to be able to reuse the
+ * VtkOutputModule which expects a VelocityOutput class.
+ */
+class GeomechanicsVelocityOutput
+{
+public:
+ //! The constructor
+ template< typename... Args >
+ GeomechanicsVelocityOutput(Args&&... args) {}
+
+ //! Output is currently disabled (not implemented)
+ static constexpr bool enableOutput() { return false; }
+
+ //! There is always only one solid phase
+ static constexpr int numPhaseVelocities() { return 1; }
+
+ //! Returns the name of phase for which velocity is computed
+ static std::string phaseName(int phaseIdx)
+ { DUNE_THROW(Dune::NotImplemented, "Velocity output for geomechanical models."); }
+
+ //! Calculate the velocities for the scvs in the element
+ template< typename... Args >
+ void calculateVelocity(Args... args)
+ { DUNE_THROW(Dune::NotImplemented, "Velocity output for geomechanical models."); }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
index 708d5406bd1f5f72be137d8f49722b8ae1ec9514..da5cb25c25e4050995c44f9c03f81e6c6bf683a2 100644
--- a/dumux/io/vtkoutputmodule.hh
+++ b/dumux/io/vtkoutputmodule.hh
@@ -68,11 +68,9 @@ class VtkOutputModule
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
- using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
- static constexpr int numPhases = ModelTraits::numPhases();
+ static constexpr int numPhaseVelocities = VelocityOutput::numPhaseVelocities();
using VV = typename GridVariables::VolumeVariables;
- using FluidSystem = typename VV::FluidSystem;
using Scalar = typename GridVariables::Scalar;
using GridView = typename FVGridGeometry::GridView;
@@ -264,7 +262,7 @@ private:
// instatiate the velocity output
VelocityOutput velocityOutput(problem_, gridGeom_, gridVariables_, sol_);
- std::array<std::vector<VelocityVector>, numPhases> velocity;
+ std::array<std::vector<VelocityVector>, numPhaseVelocities> velocity;
// process rank
static bool addProcessRank = getParamFromGroup<bool>(paramGroup_, "Vtk.AddProcessRank");
@@ -292,7 +290,7 @@ private:
if (velocityOutput.enableOutput())
{
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
{
if(isBox && dim == 1)
velocity[phaseIdx].resize(numCells);
@@ -344,7 +342,7 @@ private:
// velocity output
if (velocityOutput.enableOutput())
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
//! the rank
@@ -379,17 +377,17 @@ private:
{
if (isBox && dim > 1)
{
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
sequenceWriter_.addVertexData( Field(gridGeom_.gridView(), gridGeom_.vertexMapper(), velocity[phaseIdx],
- "velocity_" + std::string(FluidSystem::phaseName(phaseIdx+phaseIdxOffset)) + " (m/s)",
+ "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
/*numComp*/dimWorld, /*codim*/dim).get() );
}
// cell-centered models
else
{
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
sequenceWriter_.addCellData( Field(gridGeom_.gridView(), gridGeom_.elementMapper(), velocity[phaseIdx],
- "velocity_" + std::string(FluidSystem::phaseName(phaseIdx+phaseIdxOffset)) + " (m/s)",
+ "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
/*numComp*/dimWorld, /*codim*/0).get() );
}
}
@@ -433,7 +431,7 @@ private:
// instatiate the velocity output
VelocityOutput velocityOutput(problem_, gridGeom_, gridVariables_, sol_);
- std::array<std::vector<VelocityVector>, numPhases> velocity;
+ std::array<std::vector<VelocityVector>, numPhaseVelocities> velocity;
// process rank
static bool addProcessRank = getParamFromGroup<bool>(paramGroup_, "Vtk.AddProcessRank");
@@ -463,7 +461,7 @@ private:
if (velocityOutput.enableOutput())
{
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
{
if(isBox && dim == 1)
velocity[phaseIdx].resize(numCells);
@@ -521,7 +519,7 @@ private:
// velocity output
if (velocityOutput.enableOutput())
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
//! the rank
@@ -547,16 +545,16 @@ private:
{
// node-wise velocities
if (dim > 1)
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
sequenceWriter_.addVertexData( Field(gridGeom_.gridView(), gridGeom_.vertexMapper(), velocity[phaseIdx],
- "velocity_" + std::string(FluidSystem::phaseName(phaseIdx+phaseIdxOffset)) + " (m/s)",
+ "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
/*numComp*/dimWorld, /*codim*/dim).get() );
// cell-wise velocities
else
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
sequenceWriter_.addCellData( Field(gridGeom_.gridView(), gridGeom_.elementMapper(), velocity[phaseIdx],
- "velocity_" + std::string(FluidSystem::phaseName(phaseIdx+phaseIdxOffset)) + " (m/s)",
+ "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
/*numComp*/dimWorld, /*codim*/0).get() );
}
diff --git a/dumux/material/spatialparams/fvelastic.hh b/dumux/material/spatialparams/fvelastic.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2a2430cdb3ef1a548617d25f3566c5ca08c1dbff
--- /dev/null
+++ b/dumux/material/spatialparams/fvelastic.hh
@@ -0,0 +1,143 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The base class for spatial parameters of linear elastic geomechanical problems
+ */
+#ifndef DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH
+#define DUMUX_GEOMECHANICS_ELASTIC_FV_SPATIAL_PARAMS_HH
+
+#include <dumux/common/typetraits/isvalid.hh>
+#include <dumux/geomechanics/lameparams.hh>
+
+namespace Dumux {
+
+#ifndef DOXYGEN
+namespace Detail {
+// helper struct detecting if the user-defined spatial params class has a lameParamsAtPos function
+// for g++ > 5.3, this can be replaced by a lambda
+template<class GlobalPosition>
+struct hasLameParamsAtPos
+{
+ template<class SpatialParams>
+ auto operator()(const SpatialParams& a)
+ -> decltype(a.lameParamsAtPos(std::declval<GlobalPosition>()))
+ {};
+};
+
+} // end namespace Detail
+#endif
+
+/*!
+ * \ingroup Geomechanics
+ * \brief The base class for spatial parameters of linear elastic geomechanical problems
+ */
+template<class Scalar, class FVGridGeometry, class Implementation>
+class FVSpatialParamsElastic
+{
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ //! The constructor
+ FVSpatialParamsElastic(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : fvGridGeometry_(fvGridGeometry)
+ {}
+
+ /*!
+ * \brief Function for defining the solid volume fraction.
+ * That is possibly solution dependent.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \param compIdx The solid component index
+ * \return the volume fraction of the inert solid component with index compIdx
+ */
+ template<class SolidSystem, class ElementSolution>
+ Scalar inertVolumeFraction(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol,
+ int compIdx) const
+ {
+ static_assert(SolidSystem::isInert(), "Elastic model can only be used with inert solid systems");
+
+ // when there is only one component, the volume fraction must be one
+ if (SolidSystem::numInertComponents == 1)
+ return 1.0;
+
+ // otherwise we require the user to define the solid composition
+ return asImp_().template inertVolumeFractionAtPos<SolidSystem>(scv.center(), compIdx);
+ }
+
+ /*!
+ * \brief Function for defining the solid volume fraction.
+ * That is possibly solution dependent.
+ *
+ * \param globalPos The global position
+ * \param compIdx The solid component index
+ * \return the volume fraction of the inert solid component with index compIdx
+ */
+ template<class SolidSystem>
+ Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const
+ { DUNE_THROW(Dune::InvalidStateException, "The spatial parameters do not provide inertVolumeFractionAtPos() method."); }
+
+ /*!
+ * \brief Define the Lame parameters
+ * \note It is possibly solution dependent.
+ *
+ * \param element The current element
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return lame parameters
+ * \todo TODO Could the lame parameters also be heterogeneous inside element
+ * (i.e. pass scv as additional argument to this function)?
+ * This would need appropriate adjustments in Hooke's law.
+ */
+ template<class ElementSolution>
+ decltype(auto) lameParams(const Element& element,
+ const ElementSolution& elemSol) const
+ {
+ static_assert(decltype(isValid(Detail::hasLameParamsAtPos<GlobalPosition>())(this->asImp_()))::value," \n\n"
+ " Your spatial params class has to either implement\n\n"
+ " const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const\n\n"
+ " or overload this function\n\n"
+ " template<class ElementSolution>\n"
+ " const LameParams& lameParams(const Element& element,\n"
+ " const ElementSolution& elemSol) const\n\n");
+
+ return asImp_().lameParamsAtPos(element.geometry().center());
+ }
+
+ //! The finite volume grid geometry
+ const FVGridGeometry& fvGridGeometry() const { return *fvGridGeometry_; }
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+private:
+ std::shared_ptr<const FVGridGeometry> fvGridGeometry_;
+};
+}
+#endif
diff --git a/dumux/porousmediumflow/fluxvariablescache.hh b/dumux/porousmediumflow/fluxvariablescache.hh
index 59654c98672385a444547ce25a66e3a17c756b66..2dc87b9d4e6aeb0c8de09720d47843f84261ecbf 100644
--- a/dumux/porousmediumflow/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/fluxvariablescache.hh
@@ -23,11 +23,10 @@
#ifndef DUMUX_POROUSMEDIUM_FLUXVARIABLESCACHE_HH
#define DUMUX_POROUSMEDIUM_FLUXVARIABLESCACHE_HH
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
-
#include <dumux/common/properties.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/discretization/fluxvariablescaching.hh>
+#include <dumux/discretization/box/fluxvariablescache.hh>
namespace Dumux {
@@ -51,75 +50,12 @@ class PorousMediumFluxVariablesCacheImplementation;
template<class TypeTag>
using PorousMediumFluxVariablesCache = PorousMediumFluxVariablesCacheImplementation<TypeTag, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
-//! We only store discretization-related quantities for the box method.
-//! Thus, we need no physics-dependent specialization.
+//! We only store discretization-related quantities for the box method. Thus, we need no
+//! physics-dependent specialization and simply inherit from the physics-independent implementation.
template<class TypeTag>
class PorousMediumFluxVariablesCacheImplementation<TypeTag, DiscretizationMethod::box>
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = std::vector<IndexType>;
- using TransmissibilityVector = std::vector<IndexType>;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-
- using CoordScalar = typename GridView::ctype;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
- using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
- using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
- using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
- using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
- using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
-
-public:
-
- void update(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace &scvf)
- {
- const auto geometry = element.geometry();
- const auto& localBasis = fvGeometry.feLocalBasis();
-
- // evaluate shape functions and gradients at the integration point
- const auto ipLocal = geometry.local(scvf.ipGlobal());
- jacInvT_ = geometry.jacobianInverseTransposed(ipLocal);
- localBasis.evaluateJacobian(ipLocal, shapeJacobian_);
- localBasis.evaluateFunction(ipLocal, shapeValues_); // shape values for rho
-
- // compute the gradN at for every scv/dof
- gradN_.resize(fvGeometry.numScv());
- for (const auto& scv: scvs(fvGeometry))
- jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.indexInElement()]);
- }
-
- const std::vector<ShapeJacobian>& shapeJacobian() const
- { return shapeJacobian_; }
-
- const std::vector<ShapeValue>& shapeValues() const
- { return shapeValues_; }
-
- const JacobianInverseTransposed& jacInvT() const
- { return jacInvT_; }
-
- const GlobalPosition& gradN(unsigned int scvIdxInElement) const
- { return gradN_[scvIdxInElement]; }
-
-private:
- std::vector<GlobalPosition> gradN_;
- std::vector<ShapeJacobian> shapeJacobian_;
- std::vector<ShapeValue> shapeValues_;
- JacobianInverseTransposed jacInvT_;
-};
+: public BoxFluxVariablesCache<typename GET_PROP_TYPE(TypeTag, Scalar), typename GET_PROP_TYPE(TypeTag, FVGridGeometry)>
+{};
// the following classes choose the cache type: empty if the law disabled and the law's cache if it's enabled
// if advections is disabled the advection type is still instatiated if we use std::conditional_t and has to be a full type
diff --git a/dumux/porousmediumflow/velocityoutput.hh b/dumux/porousmediumflow/velocityoutput.hh
index 65101339f5871303d31219f993f61ad94eadceab..a80fbcb5a5cabf2e9f063aeb512ff7d3a9f9b71e 100644
--- a/dumux/porousmediumflow/velocityoutput.hh
+++ b/dumux/porousmediumflow/velocityoutput.hh
@@ -48,6 +48,7 @@ class PorousMediumFlowVelocityOutput
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -106,6 +107,14 @@ public:
bool enableOutput()
{ return velocityOutput_; }
+ // returns the name of the phase for a given index
+ static std::string phaseName(int phaseIdx)
+ { return FluidSystem::phaseName(phaseIdx); }
+
+ // returns the number of phase velocities computed by this class
+ static constexpr int numPhaseVelocities()
+ { return GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases(); }
+
// The following SFINAE enable_if usage allows compilation, even if only a
//
// boundaryTypes(const Element&, const scv&)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 43b16a8a440257d19336aa9e0098d056584cc6ce..3a49a60b74a1378d9df53af83ee0fe6adca2ed04 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory("common")
+add_subdirectory("geomechanics")
add_subdirectory("freeflow")
add_subdirectory("io")
add_subdirectory("material")
diff --git a/test/geomechanics/CMakeLists.txt b/test/geomechanics/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f857d3401ac26c472d2eb6b650ee0177b4f58b5d
--- /dev/null
+++ b/test/geomechanics/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("elastic")
diff --git a/test/geomechanics/elastic/CMakeLists.txt b/test/geomechanics/elastic/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..77e6742c6055cab028604454f00ffe5fe8824cdb
--- /dev/null
+++ b/test/geomechanics/elastic/CMakeLists.txt
@@ -0,0 +1,19 @@
+dune_symlink_to_source_files(FILES "elastic.input")
+
+# using box and numeric differentiation
+dune_add_test(NAME test_elastic_box
+ SOURCES test_elastic.cc
+ COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleBox NUMDIFFMETHOD=DiffMethod::numeric
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/elasticbox-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/elastic-00001.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_elastic_box elastic.input")
+
+set(CMAKE_BUILD_TYPE Release)
+
+install(FILES
+problem.hh
+spatialparams.hh
+test_elastic.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/test/geomechanics/elastic)
diff --git a/test/geomechanics/elastic/elastic.input b/test/geomechanics/elastic/elastic.input
new file mode 100644
index 0000000000000000000000000000000000000000..5125b1e9690f4c673585cc523e36426208a0a592
--- /dev/null
+++ b/test/geomechanics/elastic/elastic.input
@@ -0,0 +1,13 @@
+[Grid]
+UpperRight = 1 1
+Cells = 10 10
+
+[Problem]
+Name = elastic
+EnableGravity = false
+
+[Assembly.NumericDifference]
+PriVarMagnitude = 1e5 1e5
+
+[Component]
+SolidDensity = 2700
diff --git a/test/geomechanics/elastic/problem.hh b/test/geomechanics/elastic/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fdba80616dca105057e20cbd449b0232afcae32c
--- /dev/null
+++ b/test/geomechanics/elastic/problem.hh
@@ -0,0 +1,199 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/**
+ * \file
+ * \brief Definition of a test problem for the linear elastic model
+ */
+#ifndef DUMUX_ELASTICPROBLEM_HH
+#define DUMUX_ELASTICPROBLEM_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/geomechanics/elastic/model.hh>
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/geomechanics/fvproblem.hh>
+
+#include "spatialparams.hh"
+
+namespace Dumux {
+
+template <class TypeTag>
+class ElasticProblem;
+
+namespace Properties {
+NEW_TYPE_TAG(ElasticTypeTag, INHERITS_FROM(BoxModel, Elastic));
+// Set the grid type
+SET_TYPE_PROP(ElasticTypeTag, Grid, Dune::YaspGrid<2>);
+// Set the problem property
+SET_TYPE_PROP(ElasticTypeTag, Problem, Dumux::ElasticProblem<TypeTag>);
+// The spatial parameters property
+SET_TYPE_PROP(ElasticTypeTag, SpatialParams, ElasticSpatialParams< typename GET_PROP_TYPE(TypeTag, Scalar),
+ typename GET_PROP_TYPE(TypeTag, FVGridGeometry) >);
+}
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ *
+ * \brief Problem definition for the deformation of an elastic body
+ */
+template<class TypeTag>
+class ElasticProblem : public GeomechanicsFVProblem<TypeTag>
+{
+ using ParentType = GeomechanicsFVProblem<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GradU = Dune::FieldMatrix<Scalar, dim, dimWorld>;
+
+public:
+ //! The constructor
+ ElasticProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry) {}
+
+ //! The problem name.
+ const std::string name() const { return "elastic";}
+ //! The temperature in the domain
+ static constexpr Scalar temperature() { return 273.15; }
+ //! Evaluate the initial value for a control volume.
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
+ //! Evaluate the boundary conditions for a Dirichlet boundary segment.
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
+ //! Evaluate the boundary conditions for a Neumannboundary segment.
+ PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const { return PrimaryVariables(0.0); }
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param globalPos The global position
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
+ {
+ BoundaryTypes values;
+ values.setAllDirichlet();
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ */
+ PrimaryVariables source(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume& scv) const
+ {
+ using std::sin;
+ using std::cos;
+
+ static const Scalar pi = 3.14159265358979323846;
+ const auto ipGlobal = scv.center();
+ const auto x = ipGlobal[0];
+ const auto y = ipGlobal[1];
+
+ // the lame parameters
+ const auto& lameParams = elemVolVars[scv].lameParams();
+ const auto lambda = lameParams.lambda();
+ const auto mu = lameParams.mu();
+
+ // precalculated products
+ const Scalar pi_2 = 2.0*pi;
+ const Scalar pi_2_square = pi_2*pi_2;
+ const Scalar cos_2pix = cos(pi_2*x);
+ const Scalar sin_2pix = sin(pi_2*x);
+ const Scalar cos_2piy = cos(pi_2*y);
+ const Scalar sin_2piy = sin(pi_2*y);
+
+ const Scalar dE11_dx = -2.0*sin_2piy;
+ const Scalar dE22_dx = pi_2_square*cos_2pix*cos_2piy;
+ const Scalar dE11_dy = pi_2*(1.0-2.0*x)*cos_2piy;
+ const Scalar dE22_dy = -1.0*pi_2_square*sin_2pix*sin_2piy;
+ const Scalar dE12_dy = 0.5*pi_2_square*(cos_2pix*cos_2piy - (x-x*x)*sin_2piy);
+ const Scalar dE21_dx = 0.5*((1.0-2*x)*pi_2*cos_2piy - pi_2_square*sin_2pix*sin_2piy);
+
+ // compute exact divergence of sigma
+ PrimaryVariables divSigma(0.0);
+ divSigma[Indices::momentum(/*x-dir*/0)] = lambda*(dE11_dx + dE22_dx) + 2*mu*(dE11_dx + dE12_dy);
+ divSigma[Indices::momentum(/*y-dir*/1)] = lambda*(dE11_dy + dE22_dy) + 2*mu*(dE21_dx + dE22_dy);
+ return divSigma;
+ }
+
+ /*!
+ * \brief Evaluate the exact displacement to this problem at a given position.
+ */
+ PrimaryVariables exactSolution(const GlobalPosition& globalPos) const
+ {
+ using std::sin;
+
+ static const Scalar pi = 3.14159265358979323846;
+ const auto x = globalPos[0];
+ const auto y = globalPos[1];
+
+ PrimaryVariables exact(0.0);
+ exact[Indices::momentum(/*x-dir*/0)] = (x-x*x)*sin(2*pi*y);
+ exact[Indices::momentum(/*y-dir*/1)] = sin(2*pi*x)*sin(2*pi*y);
+ return exact;
+ }
+
+ /*!
+ * \brief Evaluate the exact displacement gradient to this problem at a given position.
+ */
+ GradU exactGradient(const GlobalPosition& globalPos) const
+ {
+ using std::sin;
+ using std::cos;
+
+ static const Scalar pi = 3.14159265358979323846;
+ const auto x = globalPos[0];
+ const auto y = globalPos[1];
+
+ static constexpr int xIdx = Indices::momentum(/*x-dir*/0);
+ static constexpr int yIdx = Indices::momentum(/*y-dir*/1);
+
+ GradU exactGrad(0.0);
+ exactGrad[xIdx][xIdx] = (1-2*x)*sin(2*pi*y);
+ exactGrad[xIdx][yIdx] = (x - x*x)*2*pi*cos(2*pi*y);
+ exactGrad[yIdx][xIdx] = 2*pi*cos(2*pi*x)*sin(2*pi*y);
+ exactGrad[yIdx][yIdx] = 2*pi*sin(2*pi*x)*cos(2*pi*y);
+ return exactGrad;
+ }
+
+private:
+ static constexpr Scalar eps_ = 3e-6;
+};
+
+} //end namespace
+
+#endif
diff --git a/test/geomechanics/elastic/spatialparams.hh b/test/geomechanics/elastic/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4471085f92ba9d2033c11ef493c5b8f48dbbfe2e
--- /dev/null
+++ b/test/geomechanics/elastic/spatialparams.hh
@@ -0,0 +1,75 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Definition of the spatial parameters for the linear elasticity problem
+ */
+#ifndef DUMUX_ELASTIC_SPATIAL_PARAMS_HH
+#define DUMUX_ELASTIC_SPATIAL_PARAMS_HH
+
+#include <dumux/geomechanics/lameparams.hh>
+#include <dumux/material/spatialparams/fvelastic.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Geomechanics
+ * \ingroup Elastic
+ * \brief Definition of the spatial parameters for the linear elasticity problem
+ */
+template<class Scalar, class FVGridGeometry>
+class ElasticSpatialParams : public FVSpatialParamsElastic< Scalar,
+ FVGridGeometry,
+ ElasticSpatialParams<Scalar, FVGridGeometry> >
+{
+ using ThisType = ElasticSpatialParams<Scalar, FVGridGeometry>;
+ using ParentType = FVSpatialParamsElastic<Scalar, FVGridGeometry, ThisType>;
+
+ using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ //! export the type of the lame parameters
+ using LameParams = Dumux::LameParams<Scalar>;
+
+ //! The constructor
+ ElasticSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ lameParams_.setLambda(3e9);
+ lameParams_.setMu(3e9);
+ }
+
+ //! Define the Lame parameters
+ const LameParams& lameParamsAtPos(const GlobalPosition& globalPos) const
+ { return lameParams_; }
+
+ //! The solid phase consists of only one component here, thus, we return 1.0
+ template<class SolidSystem>
+ Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const
+ { return 1.0; }
+
+private:
+ LameParams lameParams_;
+};
+}
+#endif
diff --git a/test/geomechanics/elastic/test_elastic.cc b/test/geomechanics/elastic/test_elastic.cc
new file mode 100644
index 0000000000000000000000000000000000000000..005a0682eca9b145d91e60ce274f1a6787f03cac
--- /dev/null
+++ b/test/geomechanics/elastic/test_elastic.cc
@@ -0,0 +1,173 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the one-phase CC model
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include "problem.hh"
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+
+// main function
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(ElasticTypeTag);
+
+ // stop time for the entire computation
+ Dune::Timer timer;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x);
+
+ // intialize the vtk output module and add displacement
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ 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());
+ vtkWriter.addField(xExact, "u_exact");
+
+ // write initial solution
+ vtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // linearize & solve
+ nonLinearSolver.solve(x);
+
+ // the grid variables need to be up to date for subsequent output
+ gridVariables->update(x);
+
+ // write vtk output
+ vtkWriter.write(1.0);
+
+ // print time and say goodbye
+ const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+ if (mpiHelper.rank() == 0)
+ std::cout << "Simulation took " << timer.elapsed() << " seconds on "
+ << comm.size() << " processes.\n"
+ << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
+
+ // print parameters
+ if (mpiHelper.rank() == 0)
+ Parameters::print();
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException& e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException& e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception& e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/test/references/elasticbox-reference.vtu b/test/references/elasticbox-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..4a1cd64dddd8823444f82b8fd791d9c237103fb3
--- /dev/null
+++ b/test/references/elasticbox-reference.vtu
@@ -0,0 +1,183 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="100" NumberOfPoints="121">
+ <PointData Vectors="u">
+ <DataArray type="Float32" Name="u" NumberOfComponents="3" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0.0486269 0.363944 0
+ 0 0 0 0.0940417 0.589125 0 0 0 0 0.131474 0.588995 0
+ 0 0 0 0.156265 0.363961 0 0 0 0 0.164956 9.03744e-18 0
+ 0 0 0 0.156265 -0.363961 0 0 0 0 0.131474 -0.588995 0
+ 0 0 0 0.0940417 -0.589125 0 0 0 0 0.0486269 -0.363944 0
+ 0 0 0 0 0 0 0 0 0 0.084587 0.589637 0
+ 0.15322 0.951768 0 0.204761 0.950788 0 0.236895 0.587348 0 0.247828 3.63266e-17 0
+ 0.236895 -0.587348 0 0.204761 -0.950788 0 0.15322 -0.951768 0 0.084587 -0.589637 0
+ 0 0 0 0 0 0 0.0876187 0.590116 0 0.151736 0.951624 0
+ 0.194708 0.950165 0 0.219141 0.586836 0 0.227034 3.83134e-17 0 0.219141 -0.586836 0
+ 0.194708 -0.950165 0 0.151736 -0.951624 0 0.0876187 -0.590116 0 0 0 0
+ 0 0 0 0.0560947 0.365428 0 0.0900437 0.58885 0 0.106944 0.587692 0
+ 0.113462 0.362888 0 0.114957 7.95142e-18 0 0.113462 -0.362888 0 0.106944 -0.587692 0
+ 0.0900437 -0.58885 0 0.0560947 -0.365428 0 0 0 0 0 0 0
+ 0.00219115 0.00143062 0 -0.00779061 0.00174145 0 -0.0239508 0.00143597 0 -0.0381341 0.000793643 0
+ -0.043704 -6.0621e-17 0 -0.0381341 -0.000793643 0 -0.0239508 -0.00143597 0 -0.00779061 -0.00174145 0
+ 0.00219115 -0.00143062 0 0 0 0 0 0 0 -0.0534173 -0.362979 0
+ -0.104191 -0.585832 0 -0.147642 -0.585182 0 -0.177287 -0.361494 0 -0.187839 -1.27627e-16 0
+ -0.177287 0.361494 0 -0.147642 0.585182 0 -0.104191 0.585832 0 -0.0534173 0.362979 0
+ 0 0 0 0 0 0 -0.0895704 -0.588716 0 -0.162557 -0.949787 0
+ -0.21726 -0.948574 0 -0.251351 -0.585934 0 -0.262955 -1.50929e-16 0 -0.251351 0.585934 0
+ -0.21726 0.948574 0 -0.162557 0.949787 0 -0.0895704 0.588716 0 0 0 0
+ 0 0 0 -0.0925845 -0.589427 0 -0.161108 -0.951263 0 -0.207287 -0.950244 0
+ -0.233697 -0.587011 0 -0.242265 -1.26245e-16 0 -0.233697 0.587011 0 -0.207287 0.950244 0
+ -0.161108 0.951263 0 -0.0925845 0.589427 0 0 0 0 0 0 0
+ -0.0608212 -0.364415 0 -0.100277 -0.589403 0 -0.123287 -0.58911 0 -0.13469 -0.363994 0
+ -0.138052 -4.5718e-17 0 -0.13469 0.363994 0 -0.123287 0.58911 0 -0.100277 0.589403 0
+ -0.0608212 0.364415 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0
+ </DataArray>
+ <DataArray type="Float32" Name="u_exact" NumberOfComponents="3" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0.0529007 0.345491 0
+ 0 0 0 0.0940456 0.559017 0 0 0 0 0.123435 0.559017 0
+ 0 0 0 0.141068 0.345491 0 0 0 0 0.146946 7.19829e-17 0
+ 0 0 0 0.141068 -0.345491 0 0 -0 0 0.123435 -0.559017 0
+ 0 -0 0 0.0940456 -0.559017 0 0 -0 0 0.0529007 -0.345491 0
+ 0 -0 0 0 -1.43966e-16 0 0 0 0 0.0855951 0.559017 0
+ 0.152169 0.904508 0 0.199722 0.904508 0 0.228254 0.559017 0 0.237764 1.16471e-16 0
+ 0.228254 -0.559017 0 0.199722 -0.904508 0 0.152169 -0.904508 0 0.0855951 -0.559017 0
+ 0 -2.32942e-16 0 0 0 0 0.0855951 0.559017 0 0.152169 0.904508 0
+ 0.199722 0.904508 0 0.228254 0.559017 0 0.237764 1.16471e-16 0 0.228254 -0.559017 0
+ 0.199722 -0.904508 0 0.152169 -0.904508 0 0.0855951 -0.559017 0 0 -2.32942e-16 0
+ 0 0 0 0.0529007 0.345491 0 0.0940456 0.559017 0 0.123435 0.559017 0
+ 0.141068 0.345491 0 0.146946 7.19829e-17 0 0.141068 -0.345491 0 0.123435 -0.559017 0
+ 0.0940456 -0.559017 0 0.0529007 -0.345491 0 0 -1.43966e-16 0 0 0 0
+ 1.10218e-17 7.19829e-17 0 1.95943e-17 1.16471e-16 0 2.57176e-17 1.16471e-16 0 2.93915e-17 7.19829e-17 0
+ 3.06162e-17 1.49976e-32 0 2.93915e-17 -7.19829e-17 0 2.57176e-17 -1.16471e-16 0 1.95943e-17 -1.16471e-16 0
+ 1.10218e-17 -7.19829e-17 0 0 -2.99952e-32 0 0 0 0 -0.0529007 -0.345491 0
+ -0.0940456 -0.559017 0 -0.123435 -0.559017 0 -0.141068 -0.345491 0 -0.146946 -7.19829e-17 0
+ -0.141068 0.345491 0 -0.123435 0.559017 0 -0.0940456 0.559017 0 -0.0529007 0.345491 0
+ 0 1.43966e-16 0 -0 -0 0 -0.0855951 -0.559017 0 -0.152169 -0.904508 0
+ -0.199722 -0.904508 0 -0.228254 -0.559017 0 -0.237764 -1.16471e-16 0 -0.228254 0.559017 0
+ -0.199722 0.904508 0 -0.152169 0.904508 0 -0.0855951 0.559017 0 -0 2.32942e-16 0
+ -0 -0 0 -0.0855951 -0.559017 0 -0.152169 -0.904508 0 -0.199722 -0.904508 0
+ -0.228254 -0.559017 0 -0.237764 -1.16471e-16 0 -0.228254 0.559017 0 -0.199722 0.904508 0
+ -0.152169 0.904508 0 -0.0855951 0.559017 0 -0 2.32942e-16 0 -0 -0 0
+ -0.0529007 -0.345491 0 -0.0940456 -0.559017 0 -0.123435 -0.559017 0 -0.141068 -0.345491 0
+ -0.146946 -7.19829e-17 0 -0.141068 0.345491 0 -0.123435 0.559017 0 -0.0940456 0.559017 0
+ -0.0529007 0.345491 0 -0 1.43966e-16 0 -0 -0 0 -2.20436e-17 -1.43966e-16 0
+ -3.91887e-17 -2.32942e-16 0 -5.14352e-17 -2.32942e-16 0 -5.8783e-17 -1.43966e-16 0 -6.12323e-17 -2.99952e-32 0
+ -5.8783e-17 1.43966e-16 0 -5.14352e-17 2.32942e-16 0 -3.91887e-17 2.32942e-16 0 -2.20436e-17 1.43966e-16 0
+ -0 5.99904e-32 0
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 0.1 0 0 0 0.1 0 0.1 0.1 0
+ 0.2 0 0 0.2 0.1 0 0.3 0 0 0.3 0.1 0
+ 0.4 0 0 0.4 0.1 0 0.5 0 0 0.5 0.1 0
+ 0.6 0 0 0.6 0.1 0 0.7 0 0 0.7 0.1 0
+ 0.8 0 0 0.8 0.1 0 0.9 0 0 0.9 0.1 0
+ 1 0 0 1 0.1 0 0 0.2 0 0.1 0.2 0
+ 0.2 0.2 0 0.3 0.2 0 0.4 0.2 0 0.5 0.2 0
+ 0.6 0.2 0 0.7 0.2 0 0.8 0.2 0 0.9 0.2 0
+ 1 0.2 0 0 0.3 0 0.1 0.3 0 0.2 0.3 0
+ 0.3 0.3 0 0.4 0.3 0 0.5 0.3 0 0.6 0.3 0
+ 0.7 0.3 0 0.8 0.3 0 0.9 0.3 0 1 0.3 0
+ 0 0.4 0 0.1 0.4 0 0.2 0.4 0 0.3 0.4 0
+ 0.4 0.4 0 0.5 0.4 0 0.6 0.4 0 0.7 0.4 0
+ 0.8 0.4 0 0.9 0.4 0 1 0.4 0 0 0.5 0
+ 0.1 0.5 0 0.2 0.5 0 0.3 0.5 0 0.4 0.5 0
+ 0.5 0.5 0 0.6 0.5 0 0.7 0.5 0 0.8 0.5 0
+ 0.9 0.5 0 1 0.5 0 0 0.6 0 0.1 0.6 0
+ 0.2 0.6 0 0.3 0.6 0 0.4 0.6 0 0.5 0.6 0
+ 0.6 0.6 0 0.7 0.6 0 0.8 0.6 0 0.9 0.6 0
+ 1 0.6 0 0 0.7 0 0.1 0.7 0 0.2 0.7 0
+ 0.3 0.7 0 0.4 0.7 0 0.5 0.7 0 0.6 0.7 0
+ 0.7 0.7 0 0.8 0.7 0 0.9 0.7 0 1 0.7 0
+ 0 0.8 0 0.1 0.8 0 0.2 0.8 0 0.3 0.8 0
+ 0.4 0.8 0 0.5 0.8 0 0.6 0.8 0 0.7 0.8 0
+ 0.8 0.8 0 0.9 0.8 0 1 0.8 0 0 0.9 0
+ 0.1 0.9 0 0.2 0.9 0 0.3 0.9 0 0.4 0.9 0
+ 0.5 0.9 0 0.6 0.9 0 0.7 0.9 0 0.8 0.9 0
+ 0.9 0.9 0 1 0.9 0 0 1 0 0.1 1 0
+ 0.2 1 0 0.3 1 0 0.4 1 0 0.5 1 0
+ 0.6 1 0 0.7 1 0 0.8 1 0 0.9 1 0
+ 1 1 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 4 6 7 5
+ 6 8 9 7 8 10 11 9 10 12 13 11
+ 12 14 15 13 14 16 17 15 16 18 19 17
+ 18 20 21 19 2 3 23 22 3 5 24 23
+ 5 7 25 24 7 9 26 25 9 11 27 26
+ 11 13 28 27 13 15 29 28 15 17 30 29
+ 17 19 31 30 19 21 32 31 22 23 34 33
+ 23 24 35 34 24 25 36 35 25 26 37 36
+ 26 27 38 37 27 28 39 38 28 29 40 39
+ 29 30 41 40 30 31 42 41 31 32 43 42
+ 33 34 45 44 34 35 46 45 35 36 47 46
+ 36 37 48 47 37 38 49 48 38 39 50 49
+ 39 40 51 50 40 41 52 51 41 42 53 52
+ 42 43 54 53 44 45 56 55 45 46 57 56
+ 46 47 58 57 47 48 59 58 48 49 60 59
+ 49 50 61 60 50 51 62 61 51 52 63 62
+ 52 53 64 63 53 54 65 64 55 56 67 66
+ 56 57 68 67 57 58 69 68 58 59 70 69
+ 59 60 71 70 60 61 72 71 61 62 73 72
+ 62 63 74 73 63 64 75 74 64 65 76 75
+ 66 67 78 77 67 68 79 78 68 69 80 79
+ 69 70 81 80 70 71 82 81 71 72 83 82
+ 72 73 84 83 73 74 85 84 74 75 86 85
+ 75 76 87 86 77 78 89 88 78 79 90 89
+ 79 80 91 90 80 81 92 91 81 82 93 92
+ 82 83 94 93 83 84 95 94 84 85 96 95
+ 85 86 97 96 86 87 98 97 88 89 100 99
+ 89 90 101 100 90 91 102 101 91 92 103 102
+ 92 93 104 103 93 94 105 104 94 95 106 105
+ 95 96 107 106 96 97 108 107 97 98 109 108
+ 99 100 111 110 100 101 112 111 101 102 113 112
+ 102 103 114 113 103 104 115 114 104 105 116 115
+ 105 106 117 116 106 107 118 117 107 108 119 118
+ 108 109 120 119
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40 44 48
+ 52 56 60 64 68 72 76 80 84 88 92 96
+ 100 104 108 112 116 120 124 128 132 136 140 144
+ 148 152 156 160 164 168 172 176 180 184 188 192
+ 196 200 204 208 212 216 220 224 228 232 236 240
+ 244 248 252 256 260 264 268 272 276 280 284 288
+ 292 296 300 304 308 312 316 320 324 328 332 336
+ 340 344 348 352 356 360 364 368 372 376 380 384
+ 388 392 396 400
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>