diff --git a/dumux/CMakeLists.txt b/dumux/CMakeLists.txt
index a8e7f0f33e13e29735b6d997152bd73f9b74f514..87cb87d8609fb96482771c7c36cd522ef43d221c 100644
--- a/dumux/CMakeLists.txt
+++ b/dumux/CMakeLists.txt
@@ -1,13 +1,11 @@
add_subdirectory("common")
add_subdirectory("discretization")
add_subdirectory("freeflow")
-add_subdirectory("geomechanics")
add_subdirectory("implicit")
add_subdirectory("io")
add_subdirectory("linear")
add_subdirectory("material")
add_subdirectory("mixeddimension")
-add_subdirectory("multidomain")
add_subdirectory("nonlinear")
add_subdirectory("parallel")
add_subdirectory("porousmediumflow")
diff --git a/dumux/geomechanics/CMakeLists.txt b/dumux/geomechanics/CMakeLists.txt
deleted file mode 100644
index 317ac1cf3564d3219e6211b0faf40a852dce5665..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/CMakeLists.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-add_subdirectory("el1p2c")
-add_subdirectory("el2p")
-add_subdirectory("elastic")
-add_subdirectory("implicit")
-add_subdirectory("constitutivelaws")
\ No newline at end of file
diff --git a/dumux/geomechanics/constitutivelaws/CMakeLists.txt b/dumux/geomechanics/constitutivelaws/CMakeLists.txt
deleted file mode 100644
index 513c666aeefecd41e73de8899ee9b040504d7a13..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/constitutivelaws/CMakeLists.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-
-#install headers
-install(FILES
-hokeslaw.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics/constitutivelaws)
diff --git a/dumux/geomechanics/constitutivelaws/hookeslaw.hh b/dumux/geomechanics/constitutivelaws/hookeslaw.hh
deleted file mode 100644
index 8e425a76f5ed0e5861295a9ea75087b088911ee7..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/constitutivelaws/hookeslaw.hh
+++ /dev/null
@@ -1,313 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \brief This file contains the data which is required to calculate
- * the mechanic stresses according to Hooke's law.
- */
-#ifndef DUMUX_GEOMECHANICS_HOOKES_LAW_HH
-#define DUMUX_GEOMECHANICS_HOOKES_LAW_HH
-
-#include
-
-#include
-#include
-
-#include
-
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-}
-
-/*!
- * \ingroup CCTpfaHookesLaw
- * \brief Evaluates the stresses, tractions and compressions on a face according to Hooke's law.
- * Specializations are given for the different discretization methods.
- */
-template
-class HookesLaw
-{};
-
-template
-class HookesLaw::type >
-{
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::IndexSet::IndexType IndexType;
- typedef typename std::vector Stencil;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension} ;
- static constexpr int voigtDim = 0.5*(dim*dim+dim);
-
- typedef Dune::FieldMatrix StiffnessMatrix;
- typedef Dune::FieldVector VoigtVector;
- typedef Dune::FieldMatrix DimMatrix;
- typedef Dune::FieldVector DimVector;
-
- struct FaceData
- {
- Scalar insideLambda, insideMu;
- DimVector insideAlpha, insideN, insideU;
-
- Scalar outsideLambda, outsideMu;
- DimVector outsideAlpha, outsideN, outsideU;
-
- bool valueSet;
-
- FaceData()
- {
- valueSet = false;
- }
- };
-
-public:
-
- static DimVector stressVector(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- DimMatrix sigma = calculateSigma_(problem, scvFace);
-
- // calculate Sigma*n
- DimVector stressVec(0.0);
- sigma.mv(scvFace.unitOuterNormal(), stressVec);
- stressVec *= scvFace.area();
-
- return stressVec;
- }
-
- static DimMatrix stressTensor(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- return calculateSigma_(problem, scvFace);
- }
-
- static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- std::vector stencil;
- if (!scvFace.boundary())
- {
- stencil.push_back(scvFace.insideScvIdx());
- stencil.push_back(scvFace.outsideScvIdx());
- }
- else
- stencil.push_back(scvFace.insideScvIdx());
-
- return stencil;
- }
-
- static DimMatrix calculateInversA(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- FaceData faceData = obtainFaceData_(problem, scvFace);
-
- DimMatrix inversA(0.0);
- addEntriesToMatrix_(faceData.insideLambda, faceData.insideMu, faceData.insideAlpha, faceData.insideN, inversA);
- addEntriesToMatrix_(faceData.outsideLambda, faceData.outsideMu, faceData.outsideAlpha, faceData.outsideN, inversA);
- inversA.invert();
-
- return inversA;
- }
-
- static DimVector interpolateFaceDisplacement(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- FaceData faceData = obtainFaceData_(problem, scvFace);
- return interpolateFaceDisplacement_(problem, scvFace, faceData);
- }
-
-private:
-
- static FaceData obtainFaceData_(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- FaceData container;
-
- const auto insideScvIdx = scvFace.insideScvIdx();
- const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
- const auto& insideVolVars = problem.model().curVolVars(insideScvIdx);
- container.insideU = insideVolVars.displacement();
- container.insideLambda = insideVolVars.lambda();
- container.insideMu = insideVolVars.mu();
- container.insideN = scvFace.unitOuterNormal();
- container.insideAlpha = scvFace.center();
- container.insideAlpha -= insideScv.center();
- container.insideAlpha /= container.insideAlpha.two_norm2();
-
- const auto outsideScvIdx = scvFace.outsideScvIdx();
- const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx);
- container.outsideU = outsideVolVars.displacement();
- container.outsideLambda = outsideVolVars.lambda();
- container.outsideMu = outsideVolVars.mu();
- container.outsideN = scvFace.unitOuterNormal();
- container.outsideN *= -1;
- if (scvFace.boundary())
- container.outsideAlpha = 0.0;
- else
- {
- const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
- container.outsideAlpha = scvFace.center();
- container.outsideAlpha -= outsideScv.center();
- container.outsideAlpha /= container.outsideAlpha.two_norm2();
- }
-
- container.valueSet = true;
-
- return container;
- }
-
- static DimMatrix calculateSigma_(const Problem& problem, const SubControlVolumeFace& scvFace)
- {
- DimMatrix sigma(0.0);
- StiffnessMatrix C(0.0);
- VoigtVector voigtStrain(0.0);
- VoigtVector voigtSigma(0.0);
-
- FaceData faceData = obtainFaceData_(problem, scvFace);
- DimVector faceU = interpolateFaceDisplacement_(problem, scvFace, faceData);
-
- fillStiffnessMatrix_(C, faceData.insideLambda, faceData.insideMu);
- fillStrainVector_(voigtStrain, faceData.insideAlpha, faceData.insideU, faceU);
-
- C.mv(voigtStrain, voigtSigma);
-
- if (dim == 2)
- {
- sigma[0][0] = voigtSigma[0];
- sigma[0][1] = voigtSigma[2];
- sigma[1][0] = voigtSigma[2];
- sigma[1][1] = voigtSigma[1];
- }
- else
- DUNE_THROW(Dune::NotImplemented, "dim = " << dim << " is not implemented yet");
-
- return sigma;
- }
-
- static DimVector interpolateFaceDisplacement_(const Problem& problem, const SubControlVolumeFace& scvFace, const FaceData& faceData, const bool oldSol = false)
- {
- DimVector faceU(0.0);
-
- if (!scvFace.boundary())
- {
- DimMatrix inversA(0.0);
- DimMatrix insideB(0.0);
- DimMatrix outsideB(0.0);
-
- getInversA_(problem, scvFace, faceData, inversA);
- addEntriesToMatrix_(faceData.insideLambda, faceData.insideMu, faceData.insideAlpha, faceData.insideN, insideB);
- addEntriesToMatrix_(faceData.outsideLambda, faceData.outsideMu, faceData.outsideAlpha, faceData.outsideN, outsideB);
-
- DimVector insideTmp(0.0);
- DimVector outsideTmp(0.0);
- insideB.mv(faceData.insideU, insideTmp);
- outsideB.mv(faceData.outsideU, outsideTmp);
-
- insideTmp += outsideTmp;
-
- inversA.mv(insideTmp, faceU);
- }
- else
- {
- if (!oldSol)
- {
- try { return problem.model().curVolVars(scvFace.outsideScvIdx()).displacement(); }
- catch (Dune::Exception& e)
- {
- DUNE_THROW(Dune::InvalidStateException, "Error ocurred during the displacement interpolation on a boundary scv face. Only call this method on inner scv faces or pure Dirichlet boundaries with the volvars bound to the element");
- }
- }
- else
- {
- // TODO
- DUNE_THROW(Dune::NotImplemented, "Reconstruction of the previous boundary vol vars not yet implemented");
- }
- }
-
- return faceU;
- }
-
- template
- static typename std::enable_if::type getInversA_(const Problem& problem,
- const SubControlVolumeFace& scvFace,
- const FaceData& faceData,
- DimMatrix& inversA)
- { inversA = problem.model().fluxVarsCache(scvFace).inversA(); }
-
- template
- static typename std::enable_if::type getInversA_(const Problem& problem,
- const SubControlVolumeFace& scvFace,
- const FaceData& faceData,
- DimMatrix& inversA)
- {
- addEntriesToMatrix_(faceData.insideLambda, faceData.insideMu, faceData.insideAlpha, faceData.insideN, inversA);
- addEntriesToMatrix_(faceData.outsideLambda, faceData.outsideMu, faceData.outsideAlpha, faceData.outsideN, inversA);
- inversA.invert();
- }
-
- static void addEntriesToMatrix_(const Scalar lambda, const Scalar mu, const DimVector& alpha, const DimVector& normal, DimMatrix& matrix)
- {
- if (dim == 2)
- {
- matrix[0][0] += (lambda + 2*mu)*alpha[0]*normal[0] + mu*alpha[1]*normal[1];
- matrix[0][1] += lambda*alpha[1]*normal[0] + mu*alpha[0]*normal[1];
- matrix[1][0] += mu*alpha[1]*normal[0] + lambda*alpha[0]*normal[1];
- matrix[1][1] += mu*alpha[0]*normal[0] + (lambda + 2*mu)*alpha[1]*normal[1];
- }
- else
- DUNE_THROW(Dune::NotImplemented, "dim = " << dim << " is not implemented yet");
- }
-
- static void fillStiffnessMatrix_(StiffnessMatrix& C, const Scalar lambda, const Scalar mu)
- {
- if (dim == 2)
- {
- C[0][0] = lambda + 2*mu;
- C[0][1] = lambda;
- C[0][2] = 0.0;
-
- C[1][0] = lambda;
- C[1][1] = lambda + 2*mu;
- C[1][2] = 0.0;
-
- C[2][0] = 0.0;
- C[2][1] = 0.0;
- C[2][2] = mu;
- }
- }
-
- static void fillStrainVector_(VoigtVector& strain, const DimVector& alpha, const DimVector& insideU, const DimVector& faceU)
- {
- if (dim == 2)
- {
- strain[0] = alpha[0]*(faceU[0] - insideU[0]);
- strain[1] = alpha[1]*(faceU[1] - insideU[1]);
- strain[2] = alpha[1]*(faceU[0] - insideU[0]) + alpha[0]*(faceU[1] - insideU[1]);
- }
- }
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/CMakeLists.txt b/dumux/geomechanics/el1p2c/CMakeLists.txt
deleted file mode 100644
index b1d385c69a89719def138462fc0a42ce5cc3c9ff..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/CMakeLists.txt
+++ /dev/null
@@ -1,13 +0,0 @@
-
-#install headers
-install(FILES
-elementvolumevariables.hh
-fluxvariables.hh
-indices.hh
-localjacobian.hh
-localresidual.hh
-model.hh
-properties.hh
-propertydefaults.hh
-volumevariables.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics/el1p2c)
diff --git a/dumux/geomechanics/el1p2c/elementvolumevariables.hh b/dumux/geomechanics/el1p2c/elementvolumevariables.hh
deleted file mode 100644
index 4e422d0ad11682b99f52ef274b002078eed44e60..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/elementvolumevariables.hh
+++ /dev/null
@@ -1,146 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Volume variables gathered on an element
- */
-#ifndef DUMUX_BOX_EL1P2C_ELEMENT_VOLUME_VARIABLES_HH
-#define DUMUX_BOX_EL1P2C_ELEMENT_VOLUME_VARIABLES_HH
-
-#include
-#include
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ElOnePTwoCBoxModel
- *
- * \brief This class stores an array of VolumeVariables objects, one
- * volume variables object for each of the element's vertices
- */
-template
-class ElOnePTwoCElementVolumeVariables : public BoxElementVolumeVariables
-{
- typedef BoxElementVolumeVariables ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- enum { dim = GridView::dimension };
-
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-
-public:
- /*!
- * \brief The constructor.
- */
- ElOnePTwoCElementVolumeVariables()
- { }
-
- /*!
- * \brief Construct the volume variables for all vertices of an element.
- *
- * \param problem The problem which needs to be simulated.
- * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
- * \param fvGeometry The finite volume geometry of the element
- * \param oldSol Tells whether the model's previous or current solution should be used.
- *
- * This class is required for the update of the effective porosity values at the
- * vertices since it is a function of the divergence of the solid displacement
- * at the integration points
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- bool oldSol)
- {
- ParentType::update(problem, element, fvGeometry, oldSol);
- this->updateEffPorosity(problem, element, fvGeometry);
-
- };
-
- /*!
- * \brief Update the effective porosities and the volumetric strain divU for all vertices of an element.
- *
- * \param problem The problem which needs to be simulated.
- * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
- * \param fvGeometry The finite volume geometry of the element
- *
- * This function is required for the update of the effective porosity / divU values at the
- * vertices.
- *
- * During the partial derivative calculation, changes of the solid displacement
- * at vertex i can affect effective porosities / divU of all element vertices.
- * To correctly update the effective porosities / divU of all element vertices
- * an iteration over all scv faces is required.
- * The remaining volvars are only updated for the vertex whose primary variable
- * is changed for the derivative calculation.
- */
- void updateEffPorosity(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry)
- {
- // we assert that the i-th shape function is
- // associated to the i-th vert of the element.
- int numScv = element.subEntities(dim);
-
- // number of faces which contribute to the porosity value in the sub-control volume
- std::vector numContributingFaces;
- numContributingFaces.resize(numScv);
-
- for (int scvIdx = 0; scvIdx < numScv; scvIdx++) {
- (*this)[scvIdx].effPorosity = 0.0;
- (*this)[scvIdx].divU = 0.0;
- numContributingFaces[scvIdx] = 0.0;
- }
- for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++)
- {
- // evaluate the gradients at the IPs for each subcontrol volume face
- FluxVariables fluxVars;
- fluxVars.update(problem,
- element,
- fvGeometry,
- fIdx,
- *this);
-
- numContributingFaces[fluxVars.face().i] += 1;
- numContributingFaces[fluxVars.face().j] += 1;
-
- // average value for the effective porosity
- (*this)[fluxVars.face().i].effPorosity += fluxVars.effPorosity();
- (*this)[fluxVars.face().j].effPorosity += fluxVars.effPorosity();
- // average value for the volumetric strain
- (*this)[fluxVars.face().i].divU += fluxVars.divU();
- (*this)[fluxVars.face().j].divU += fluxVars.divU();
-
- }
- for (int scvIdx = 0; scvIdx < numScv; scvIdx++) {
- (*this)[scvIdx].effPorosity /= numContributingFaces[scvIdx];
- (*this)[scvIdx].divU /= numContributingFaces[scvIdx];
- }
- };
-
-
-};
-
-} // namespace Dumux
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/fluxvariables.hh b/dumux/geomechanics/el1p2c/fluxvariables.hh
deleted file mode 100644
index 40c5bf0b5316d79acd94dc5d69ecc7cb75460827..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/fluxvariables.hh
+++ /dev/null
@@ -1,327 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file contains the calculation of all the fluxes over the surface of the
- * finite volume that make up the volume, the mass and the momentum balance
- * for the one-phase two-component linear-elastic model.
- *
- * This means pressure, concentration and solid-displacement gradients, phase densities at
- * the integration point, etc.
- *
- * This class inherits from the one-phase two-component model FluxVariables and from the
- * linear elasticity model FluxVariables
- */
-#ifndef DUMUX_ELASTIC1P2C_FLUX_VARIABLES_HH
-#define DUMUX_ELASTIC1P2C_FLUX_VARIABLES_HH
-
-#include
-#include
-
-namespace Dumux
-{
-/*!
- * \ingroup ElOnePTwoCBoxModel
- * \ingroup ImplicitFluxVariables
- * \brief This template class contains the data which is required to
- * calculate the fluxes over the surface of the
- * finite volume that make up the volume, the mass and the momentum balance
- * for the one-phase two-component linear-elastic model.
- *
- * This means pressure, concentration and solid-displacement gradients, phase densities at
- * the integration point, etc.
- *
- */
-template
-class ElOnePTwoCFluxVariables: public ElasticFluxVariablesBase ,
- public OnePTwoCFluxVariables
-{
- friend class ElasticFluxVariablesBase; // be friends with parents
- friend class OnePTwoCFluxVariables; // be friends with parents
-
- typedef ElasticFluxVariablesBase ElasticBase;
- typedef OnePTwoCFluxVariables OnePTwoCBase;
-
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffectiveDiffusivityModel;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
- typedef typename GridView::ctype CoordScalar;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldVector GlobalPosition;
- typedef Dune::FieldVector DimVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
- /*!
- * \brief Compute / update the flux variables
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param fIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- {
- ElasticBase::update(problem, element, fvGeometry, fIdx, elemVolVars);
- OnePTwoCBase::update(problem, element, fvGeometry, fIdx, elemVolVars);
-
- dU_ = 0.0;
- dGradP_ = 0.0;
- porosity_ = 0.0;
- effPorosity_ = 0.0;
- pressure_ = 0.0;
- timeDerivUNormal_ = 0.0;
-
- elOnePTwoCGradients_(problem, element, elemVolVars);
- calculateEffectiveValues_(problem, element, elemVolVars);
- calculateDiffCoeffPM_(problem, element, elemVolVars);
- calculateDDt_(problem, element, elemVolVars);
- }
-
-public:
- /*!
- * \brief Return porosity [-] at the integration point.
- */
- Scalar porosity() const
- {
- return porosity_;
- }
-
- /*!
- * \brief Return effective porosity [-] at the integration point.
- */
- Scalar effPorosity() const
- {
- return effPorosity_;
- }
-
- /*!
- * \brief Return pressure [Pa] at the integration
- * point.
- */
- Scalar pressure() const
- {
- return pressure_;
- }
-
- /*!
- * \brief Return change of pressure gradient with time [Pa/m] at
- * integration point.
- */
- Scalar dGradP(int dimIdx) const
- {
- return dGradP_[dimIdx];
- }
-
- /*!
- * \brief Return gradient of time derivative of pressure [Pa].
- */
- Scalar timeDerivGradPNormal() const
- {
- return timeDerivGradPNormal_;
- }
-
- /*!
- * \brief Return change of u [m] with time at integration point
- * point.
- */
- Scalar dU(int dimIdx) const
- {
- return dU_[dimIdx];
- }
-
- /*!
- * \brief Return time derivative of u [m/s] in normal direction at integration point
- */
- Scalar timeDerivUNormal() const
- {
- return timeDerivUNormal_;
- }
-
- /*!
- * \brief Return porous medium diffusion coefficient [m^2]
- */
- Scalar diffCoeffPM() const
- {
- return diffCoeffPM_;
- }
-
- const SCVFace &face() const
- {
- return ElasticBase::face();
- }
-
-protected:
- // Overload the parent's methods to avoid ambiguous overloads due to multiple inheritance
- // The elastic gradients are already computed in the elastic base class update
- void calculateGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- OnePTwoCBase::calculateGradients_(problem, element, elemVolVars);
- }
-
- /*!
- * \brief Calculation of the solid displacement and pressure gradients.
- *
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemVolVars The parameters stored in the considered element
- */
- void elOnePTwoCGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // calculate gradients
- GlobalPosition tmp(0.0);
- for (int idx = 0; idx < ElasticBase::fvGeometry_().numScv; idx++) // loop over adjacent vertices
- {
- // FE gradient at vertex idx
- const DimVector &feGrad = face().grad[idx];
-
- // the gradient of the temporal pressure change (needed for stabilization term)
- tmp = feGrad;
- tmp *= elemVolVars[idx].dPressure();
- dGradP_ += tmp;
-
- // average the pressure at integration point
- pressure_ += elemVolVars[idx].pressure()
- * face().shapeValue[idx];
- // average temporal displacement change at integration point (for calculation of solid displacement velocity)
- for (int i = 0; i < dim; ++i)
- dU_[i] += elemVolVars[idx].dU(i)
- * face().shapeValue[idx];
- // average porosity at integration point
- porosity_ += elemVolVars[idx].porosity()
- * face().shapeValue[idx];
- }
- }
-
- /*!
- * \brief Calculation of the effective porosity.
- *
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemVolVars The parameters stored in the considered element
- */
- void calculateEffectiveValues_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
-
- // the effective porosity is calculated as a function of solid displacement and initial porosity
- // according to Han & Dusseault (2003)
-
- // calculate effective porosity as a function of solid displacement and initial porosity
- effPorosity_ = (porosity_ + this->divU())
- / (1 + this->divU());
- }
-
- /*!
- * \brief Calculation of the effective porous media diffusion coefficient.
- *
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemVolVars The parameters stored in the considered element
- */
- void calculateDiffCoeffPM_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- const VolumeVariables &volVarsI = elemVolVars[face().i];
- const VolumeVariables &volVarsJ = elemVolVars[face().j];
-
- const Scalar diffCoeffI =
- EffectiveDiffusivityModel::effectiveDiffusivity(volVarsI.porosity(),
- /*sat=*/1.0,
- volVarsI.diffCoeff());
-
- const Scalar diffCoeffJ =
- EffectiveDiffusivityModel::effectiveDiffusivity(volVarsJ.porosity(),
- /*sat=*/1.0,
- volVarsJ.diffCoeff());
-
- diffCoeffPM_ = harmonicMean(diffCoeffI, diffCoeffJ);
- }
-
- /*!
- * \brief Calculation of the time derivative of solid displacement and pressure gradient
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemVolVars The parameters stored in the considered element
- */
- void calculateDDt_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- Scalar dt= problem.timeManager().timeStepSize();
- DimVector tmp(0.0);
-
- //time derivative of solid displacement times normal vector
- for (int i = 0; i < dim; ++i)
- tmp[i] = dU_[i] / dt;
- timeDerivUNormal_ = tmp * face().normal;
- //time derivative of pressure gradient times normal vector
- for (int i = 0; i < dim; ++i)
- tmp[i] = dGradP_[i] / dt;
- timeDerivGradPNormal_ = tmp * face().normal;
- }
-
- //! change of solid displacement with time at integration point
- GlobalPosition dU_;
- //! change of pressure gradient with time at integration point
- GlobalPosition dGradP_;
- //! porosity at integration point
- Scalar porosity_;
- //! effective porosity at integration point
- Scalar effPorosity_;
- //! pressure at integration point
- Scalar pressure_;
- //! time derivative of solid displacement times normal vector at integration point
- Scalar timeDerivUNormal_;
- //! time derivative of pressure gradient times normal vector at integration point
- Scalar timeDerivGradPNormal_;
- //! Parameters
- Scalar diffCoeffPM_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/indices.hh b/dumux/geomechanics/el1p2c/indices.hh
deleted file mode 100644
index 6b5c446286b96ef07e679c79a0720d6ec39ef04e..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/indices.hh
+++ /dev/null
@@ -1,53 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Defines the primary variable and equation indices used by
- * the one-phase two-component linear elasticity model.
- */
-
-#ifndef DUMUX_ELASTIC1P2C_INDICES_HH
-#define DUMUX_ELASTIC1P2C_INDICES_HH
-
-#include
-#include
-
-namespace Dumux
-{
-// \{
-
-/*!
- * \ingroup ElOnePTwoCBoxModel
- * \ingroup ImplicitIndices
- * \brief The indices for the one-phase two-component linear elasticity model.
- *
- * This class inherits from the OnePTwoCIndices and from the ElasticIndices
- */
-template
-// PVOffset is set to 0 for the OnePTwoCIndices and to 2 for the ElasticIndices since
-// the first two primary variables are the primary variables of the one-phase two-component
-// model followed by the primary variables of the elastic model
-class ElOnePTwoCIndices : public OnePTwoCIndices, public ElasticIndices<2>
-{
-};
-
-} // namespace Dumux
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/localjacobian.hh b/dumux/geomechanics/el1p2c/localjacobian.hh
deleted file mode 100644
index 4184d8cd4aa2dc93c0fb74f75241a7952d40e224..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/localjacobian.hh
+++ /dev/null
@@ -1,258 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Calculates the partial derivatives of the local residual for the Jacobian of the
- * one-phase two-component linear elasticity model.
- */
-#ifndef DUMUX_EL1P2C_LOCAL_JACOBIAN_HH
-#define DUMUX_EL1P2C_LOCAL_JACOBIAN_HH
-
-#include
-
-namespace Dumux
-{
-/*!
- * \ingroup ElOnePTwoCBoxModel
- * \brief Calculates the partial derivatives of the local residual for the Jacobian
- *
- * Except for the evalPartialDerivatives function all functions are taken from the
- * base class ImplicitLocalJacobian
- */
-template
-class ElOnePTwoCLocalJacobian : public ImplicitLocalJacobian
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- enum {
- dim = GridView::dimension,
- };
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-
- // copying a local jacobian is not a good idea
- ElOnePTwoCLocalJacobian(const ElOnePTwoCLocalJacobian &);
-
-public:
- ElOnePTwoCLocalJacobian()
- {}
-
- /*!
- * \brief Compute the partial derivatives to a primary variable at
- * an degree of freedom.
- *
- * This method is overwritten here since this model requires a call of the model specific
- * elementvolumevariables which updates the effective porosities correctly.
- *
- * The default implementation of this method uses numeric
- * differentiation, i.e. forward or backward differences (2nd
- * order), or central differences (3rd order). The method used is
- * determined by the "NumericDifferenceMethod" property:
- *
- * - if the value of this property is smaller than 0, backward
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
- * \f]
- *
- * - if the value of this property is 0, central
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
- * \f]
- *
- * - if the value of this property is larger than 0, forward
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
- * \f]
- *
- * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
- * is the value of a sub-control volume's primary variable at the
- * evaluation point and \f$\epsilon\f$ is a small value larger than 0.
- *
- * \param partialDeriv The vector storing the partial derivatives of all
- * equations
- * \param storageDeriv the mass matrix contributions
- * \param col The block column index of the degree of freedom
- * for which the partial derivative is calculated.
- * Box: a sub-control volume index.
- * Cell centered: a neighbor index.
- * \param pvIdx The index of the primary variable
- * for which the partial derivative is calculated
- */
- void evalPartialDerivative_(ElementSolutionVector &partialDeriv,
- PrimaryVariables &storageDeriv,
- const int col,
- const int pvIdx)
- {
- int dofIdxGlobal;
- FVElementGeometry neighborFVGeom;
- auto neighbor = this->element_();
- if (isBox)
- {
- dofIdxGlobal = this->vertexMapper_().subIndex(this->element_(), col, dim);
-
- }
- else
- {
- neighbor = this->fvElemGeom_.neighbors[col];
- neighborFVGeom.updateInner(neighbor);
- dofIdxGlobal = this->problemPtr_->elementMapper().index(neighbor);
-
- }
-
- PrimaryVariables priVars(this->model_().curSol()[dofIdxGlobal]);
- VolumeVariables origVolVars(this->curVolVars_[col]);
-
- this->curVolVars_[col].setEvalPoint(&origVolVars);
- Scalar eps = this->numericEpsilon(col, pvIdx);
- Scalar delta = 0;
-
- if (this->numericDifferenceMethod_ >= 0) {
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- priVars[pvIdx] += eps;
- delta += eps;
-
- // calculate the residual
- if (isBox){
- this->curVolVars_[col].update(priVars,
- this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- col,
- false);
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_);
- }
- else{
- this->curVolVars_[col].update(priVars,
- this->problem_(),
- neighbor,
- neighborFVGeom,
- /*scvIdx=*/0,
- false);
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_);
- }
-
- this->localResidual().eval(this->element_(),
- this->fvElemGeom_,
- this->prevVolVars_,
- this->curVolVars_,
- this->bcTypes_);
-
- // store the residual and the storage term
- partialDeriv = this->localResidual().residual();
- if (isBox || col == 0)
- storageDeriv = this->localResidual().storageTerm()[col];
- }
- else {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv = this->residual_;
- storageDeriv = this->storageTerm_[col];
- }
-
-
- if (this->numericDifferenceMethod_ <= 0) {
- // we are not using forward differences, i.e. we don't
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- priVars[pvIdx] -= delta + eps;
- delta += eps;
-
- // calculate residual again
- if (isBox){
- this->curVolVars_[col].update(priVars,
- this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- col,
- false);
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_);
- }
- else{
- this->curVolVars_[col].update(priVars,
- this->problem_(),
- neighbor,
- neighborFVGeom,
- /*scvIdx=*/0,
- false);
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_);
- }
-
- this->localResidual().eval(this->element_(),
- this->fvElemGeom_,
- this->prevVolVars_,
- this->curVolVars_,
- this->bcTypes_);
- partialDeriv -= this->localResidual().residual();
- if (isBox || col == 0)
- storageDeriv -= this->localResidual().storageTerm()[col];
- }
- else {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv -= this->residual_;
- if (isBox || col == 0)
- storageDeriv -= this->storageTerm_[col];
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- partialDeriv /= delta;
- storageDeriv /= delta;
-
- // restore the original state of the element's volume variables
- this->curVolVars_[col] = origVolVars;
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_);
-
-#if HAVE_VALGRIND
- for (unsigned i = 0; i < partialDeriv.size(); ++i)
- Valgrind::CheckDefined(partialDeriv[i]);
-#endif
- }
-};
-}
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/localresidual.hh b/dumux/geomechanics/el1p2c/localresidual.hh
deleted file mode 100644
index 2e2e7628784362e952c95e618be6c15dc33bdaff..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/localresidual.hh
+++ /dev/null
@@ -1,388 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Element-wise calculation the local Jacobian for the linear elastic,
- * single-phase, two-component model in the fully implicit scheme.
- */
-#ifndef DUMUX_ELASTIC1P2C_LOCAL_RESIDUAL_HH
-#define DUMUX_ELASTIC1P2C_LOCAL_RESIDUAL_HH
-
-#include "properties.hh"
-
-namespace Dumux
-{
- /*!
- * \ingroup ElOnePTwoCModel
- * \ingroup ImplicitLocalResidual
- * \brief Calculate the local Jacobian for a one-phase two-component
- * flow in a linear-elastic porous medium.
- *
- * This class is used to fill the gaps in BoxLocalResidual for the
- * one-phase two-component linear elasticity model.
- */
- template
- class ElOnePTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
- {
- protected:
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- enum { dim = GridView::dimension };
- typedef Dune::FieldMatrix DimMatrix;
- typedef Dune::FieldVector DimVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
- enum {
- //phase index
- phaseIdx = Indices::phaseIdx,
- transportCompIdx = Indices::transportCompIdx
- };
- // indices of the equations
- enum {
- conti0EqIdx = Indices::conti0EqIdx,
- transportEqIdx = Indices::transportEqIdx
- };
-
- //! property that defines whether mole or mass fractions are used
- static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
-
- public:
- /*!
- * \brief Constructor. Sets the upwind weight.
- */
- ElOnePTwoCLocalResidual()
- {
- // retrieve the upwind weight for the mass conservation equations. Use the value
- // specified via the property system as default, and overwrite
- // it by the run-time parameter from the Dune::ParameterTree
- upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
- // retrieve the property which defines if the stabilization terms in the mass balance
- // equations are switched on. Use the value specified via the property system as default,
- // and overwrite it by the run-time parameter from the Dune::ParameterTree
- withStabilization_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, WithStabilization);
- };
-
- /*!
- * \brief Evaluate the amount of all conservation quantities
- * (e.g. phase mass) within a finite volume.
- *
- * \param storage The mass of the component within the sub-control volume
- * \param scvIdx The index of the considered face of the sub-control volume
- * \param usePrevSol Evaluate function with solution of current or previous time step
- */
- void computeStorage(PrimaryVariables &storage, int scvIdx,
- bool usePrevSol) const
- {
- // if flag usePrevSol is set, the solution from the previous
- // time step is used, otherwise the current solution is
- // used. The secondary variables are used accordingly. This
- // is required to compute the derivative of the storage term
- // using the implicit euler method.
- const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
- const VolumeVariables &volVars = elemVolVars[scvIdx];
-
- storage = 0;
- // this model assumes incompressible fluids and solids only the bulk
- // density i.e. the ratio of solid phase and pore fluid can vary
- // storage term of continuity equation
- storage[conti0EqIdx] += volVars.divU;
-
- if(useMoles)
- {
- // storage term of the transport equation - mole fractions
- storage[transportEqIdx] += volVars.moleFraction(transportCompIdx)
- * volVars.effPorosity;
- }
- else
- {
- //storage term of the transport equation - mass fractions
- storage[transportEqIdx] += volVars.massFraction(transportCompIdx)
- * volVars.effPorosity;
- }
- }
-
- /*!
- * \brief Evaluate the mass flux over a face of a sub-control
- * volume.
- *
- * \param flux The flux over the SCV (sub-control-volume) face for each component
- * \param fIdx The index of the considered face of the sub control volume
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- void computeFlux(PrimaryVariables &flux, int fIdx, const bool onBoundary=false) const
- {
- flux = 0;
- FluxVariables fluxVars;
- fluxVars.update(this->problem_(),
- this->element_(),
- this->fvGeometry_(),
- fIdx,
- this->curVolVars_());
-
- this->computeAdvectiveFlux(flux, fluxVars);
- this->computeDiffusiveFlux(flux, fluxVars);
- this->computeStresses(flux, fluxVars, fIdx);
- }
-
- /*!
- * \brief Evaluates the advective mass flux of all phases over
- * a face of a subcontrol volume.
- *
- * \param flux The advective flux over the sub-control-volume face for each component
- * \param fluxVars The flux variables at the current SCV
- */
- void computeAdvectiveFlux(PrimaryVariables &flux,
- const FluxVariables &fluxVars) const
- {
- ////////
- // advective fluxes of all components in all phases
- ////////
-
- // data attached to upstream and the downstream vertices
- // of the current phase
- const VolumeVariables &up =
- this->curVolVars_(fluxVars.upstreamIdx());
- const VolumeVariables &dn =
- this->curVolVars_(fluxVars.downstreamIdx());
-
- // calculate the stabilization term which helps in case of stability problems
- // e.g. observed for small time steps (according to G.Aguilar, F.Gaspar, F.Lisbona
- // and C.Rodrigo (2008))
- Scalar stabilizationTerm(0.0);
- if(withStabilization_){
- // calculate distance h between nodes i and j
- const auto geometry = this->element_().geometry();
- DimVector hVec = geometry.corner(fluxVars.face().j)
- - geometry.corner(fluxVars.face().i);
- Scalar h = hVec.two_norm();
- stabilizationTerm = (h * h) /
- (4 * (fluxVars.lambda()
- + 2 * fluxVars.mu()));
-
- stabilizationTerm *= fluxVars.timeDerivGradPNormal();
- }
-
-
- // flux for mass balance of the solid-fluid mixture
- // KmvpNormal is the Darcy velocity multiplied with the normal vector,
- // calculated in 1p2cfluxvariables.hh
- flux[conti0EqIdx] +=
- fluxVars.KmvpNormal() *
- (( upwindWeight_)/up.viscosity()
- +
- ((1 - upwindWeight_)/dn.viscosity()));
-
-
- // stabilization term
- if(withStabilization_)
- flux[conti0EqIdx] -= stabilizationTerm;
-
- if(useMoles)
- {
- // mass flux of the dissolved second component - massfraction
- // advective flux of the component
- flux[transportEqIdx] +=
- fluxVars.KmvpNormal() *
- (( upwindWeight_)* up.moleFraction(transportCompIdx)/up.viscosity()
- +
- (1 - upwindWeight_)* dn.moleFraction(transportCompIdx)/dn.viscosity());
-
- // flux of the dissolved second component due to solid displacement
- flux[transportEqIdx] +=
- fluxVars.timeDerivUNormal() *
- (( upwindWeight_)* up.moleFraction(transportCompIdx)
- * up.effPorosity
- +
- (1 - upwindWeight_)*dn.moleFraction(transportCompIdx)
- * up.effPorosity);
-
- // stabilization term
- if(withStabilization_)
- flux[transportEqIdx] -=
- stabilizationTerm *
- (( upwindWeight_)* up.moleFraction(transportCompIdx)
- +
- (1 - upwindWeight_)*dn.moleFraction(transportCompIdx));
- }
- else
- {
- // mass flux of the dissolved second component - massfraction
- // advective flux of the component
- flux[transportEqIdx] +=
- fluxVars.KmvpNormal() *
- (( upwindWeight_)* up.massFraction(transportCompIdx)/up.viscosity()
- +
- (1 - upwindWeight_)* dn.massFraction(transportCompIdx)/dn.viscosity());
-
- // flux of the dissolved second component due to solid displacement
- flux[transportEqIdx] +=
- fluxVars.timeDerivUNormal() *
- (( upwindWeight_)* up.massFraction(transportCompIdx)
- * up.effPorosity
- +
- (1 - upwindWeight_)*dn.massFraction(transportCompIdx)
- * up.effPorosity);
-
- // stabilization term
- if(withStabilization_)
- flux[transportEqIdx] -=
- stabilizationTerm *
- (( upwindWeight_)* up.massFraction(transportCompIdx)
- +
- (1 - upwindWeight_)*dn.massFraction(transportCompIdx));
- }
- }
-
- /*!
- * \brief Adds the diffusive mass flux of all components over
- * a face of a sub-control volume.
- *
- * \param flux The diffusive flux over the sub-control-volume face for each component
- * \param fluxVars The flux variables at the current sub-control-volume face
- */
- void computeDiffusiveFlux(PrimaryVariables &flux,
- const FluxVariables &fluxVars) const
- {
- Scalar tmp(0);
-
- // diffusive flux of second component
- if(useMoles)
- {
- // diffusive flux of the second component - mole fraction
- tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
- tmp *= fluxVars.diffCoeffPM();
-
- // dispersive flux of second component - mole fraction
- DimVector normalDisp;
- fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
- tmp -= (normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
-
- flux[transportEqIdx] += tmp;
- }
- else
- {
- // diffusive flux of the second component - mass fraction
- tmp = -(fluxVars.moleFractionGrad(transportCompIdx)*fluxVars.face().normal);
- tmp *= fluxVars.diffCoeffPM();
-
- // dispersive flux of second component - mass fraction
- DimVector normalDisp;
- fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
- tmp -= (normalDisp * fluxVars.moleFractionGrad(transportCompIdx));
-
- // convert it to a mass flux and add it
- flux[transportEqIdx] += tmp * FluidSystem::molarMass(transportCompIdx);
- }
- }
-
- /*!
- * \brief Evaluates the total stress induced by effective stresses and fluid
- * pressure in the solid fluid mixture.
- * \param stress The stress over the sub-control-volume face for each component
- * \param fluxVars The variables at the current sub-control-volume face
- * \param fIdx The index of the current sub-control-volume face
- */
- void computeStresses(PrimaryVariables &stress,
- const FluxVariables &fluxVars, const int fIdx) const
- {
- DimMatrix pressure(0.0), sigma(0.0);
- // the normal vector corresponding to the current sub-control-volume face
- const DimVector &normal(this->fvGeometry_().subContVolFace[fIdx].normal);
-
- // the pressure term of the momentum balance
- for (int i = 0; i < dim; ++i)
- pressure[i][i] += 1.0;
-
- pressure *= fluxVars.pressure();
- // effective stresses
- sigma = fluxVars.sigma();
- // calculate total stresses by subtracting the pressure
- sigma -= pressure;
-
- DimVector tmp(0.0);
- // multiply total stress tensor with normal vector of current face
- sigma.mv(normal, tmp);
-
- // set the stress term equal to the calculated vector
- for (int i = 0; i < dim; ++i)
- stress[Indices::momentum(i)] = tmp[i];
- }
-
- /*!
- * \brief Calculate the source term of the equation
- * \param source The source/sink in the SCV for each component
- * \param scvIdx The index of the vertex of the sub control volume
- *
- */
- void computeSource(PrimaryVariables &source, const int scvIdx)
- {
- source = 0;
-
- const ElementVolumeVariables &elemVolVars = this->curVolVars_();
- const VolumeVariables &volVars = elemVolVars[scvIdx];
-
- DimVector tmp1(0.0), tmp2(0.0);
-
- this->problem_().solDependentSource(source,
- this->element_(),
- this->fvGeometry_(),
- scvIdx,
- this->curVolVars_());
-
- // the gravity term of the momentum balance equation is treated as a source term
- // gravity contribution of solid matrix
- tmp1 = this->problem_().gravity();
- tmp1 *= volVars.rockDensity();
- tmp1 *= (1. - volVars.effPorosity);
-
- // gravity contribution of the fluids
- tmp2 = this->problem_().gravity();
- tmp2 *= volVars.density();
- tmp2 *= volVars.effPorosity;
-
- tmp1 += tmp2;
-
- for (int i = 0; i < dim; ++i)
- source[Indices::momentum(i)] += tmp1[i];
- }
-
- Implementation *asImp_()
- { return static_cast (this); }
- const Implementation *asImp_() const
- { return static_cast (this); }
-
- private:
- Scalar upwindWeight_;
- bool withStabilization_;
- };
-
-}
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/model.hh b/dumux/geomechanics/el1p2c/model.hh
deleted file mode 100644
index a04dbe03a79f5dcb57671af90d24a64239c39d2d..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/model.hh
+++ /dev/null
@@ -1,517 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Base class for all models which use the one-phase two-component linear elasticity model.
- * Adaption of the fully implicit scheme to the one-phase two-component linear elasticity model.
- */
-#ifndef DUMUX_ELASTIC1P2C_MODEL_HH
-#define DUMUX_ELASTIC1P2C_MODEL_HH
-
-#include "properties.hh"
-#include
-
-namespace Dumux {
-/*!
- * \ingroup ElOnePTwoCBoxModel
- * \brief Adaption of the fully implicit scheme to the one-phase two-component linear elasticity model.
- *
- * This model implements a one-phase flow of an incompressible fluid, that consists of two components.
- * The deformation of the solid matrix is described with a quasi-stationary momentum balance equation.
- * The influence of the pore fluid is accounted for through the effective stress concept (Biot 1941).
- * The total stress acting on a rock is partially supported by the rock matrix and partially supported
- * by the pore fluid. The effective stress represents the share of the total stress which is supported
- * by the solid rock matrix and can be determined as a function of the strain according to Hooke's law.
- *
- * As an equation for the conservation of momentum within the fluid phase Darcy's approach is used:
- \f[
- v = - \frac{\textbf K}{\mu}
- \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right)
- \f]
- *
- * Gravity can be enabled or disabled via the property system.
- * By inserting this into the volume balance of the solid-fluid mixture, one gets
- \f[
- \frac{\partial \text{div} \textbf{u}}{\partial t} - \text{div} \left\{
- \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right)\right\} = q \;,
- \f]
- *
- * The transport of the components \f$\kappa \in \{ w, a \}\f$ is described by the following equation:
- \f[
- \frac{ \partial \phi_{eff} X^\kappa}{\partial t}
- - \text{div} \left\lbrace
- X^\kappa \frac{{\textbf K}}{\mu} \left( \textbf{grad}\, p - \varrho_w {\textbf g} \right)
- + D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa
- - \phi_{eff} X^\kappa \frac{\partial \boldsymbol{u}}{\partial t}
- \right\rbrace = q.
- \f]
- *
- * If the model encounters stability problems, a stabilization term can be switched on. The stabilization
- * term is defined in Aguilar et al (2008):
- \f[
- \beta \text{div} \textbf{grad} \frac{\partial p}{\partial t}
- \f]
- with \f$\beta\f$:
- \f[
- \beta = h^2 / 4(\lambda + 2 \mu)
- \f]
- * where \f$h\f$ is the discretization length.
- *
- * The balance equations
- * with the stabilization term are given below:
- \f[
- \frac{\partial \text{div} \textbf{u}}{\partial t} - \text{div} \left\{
- \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho_w {\textbf g} \right)
- + \varrho_w \beta \textbf{grad} \frac{\partial p}{\partial t}
- \right\} = q \;,
- \f]
- *
- * The transport of the components \f$\kappa \in \{ w, a \}\f$ is described by the following equation:
- *
- \f[
- \frac{ \partial \phi_{eff} X^\kappa}{\partial t}
- - \text{div} \left\lbrace
- X^\kappa \frac{{\textbf K}}{\mu} \left( \textbf{grad}\, p - \varrho_w {\textbf g} \right)
- + \varrho_w X^\kappa \beta \textbf{grad} \frac{\partial p}{\partial t}
- + D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa
- - \phi_{eff} X^\kappa \frac{\partial \boldsymbol{u}}{\partial t}
- \right\rbrace = q.
- \f]
- *
- *
- * The quasi-stationary momentum balance equation is:
- \f[
- \text{div}\left( \boldsymbol{\sigma'}- p \boldsymbol{I} \right) + \left( \phi_{eff} \varrho_w + (1 - \phi_{eff}) * \varrho_s \right)
- {\textbf g} = 0 \;,
- \f]
- * with the effective stress:
- \f[
- \boldsymbol{\sigma'} = 2\,G\,\boldsymbol{\epsilon} + \lambda \,\text{tr} (\boldsymbol{\epsilon}) \, \boldsymbol{I}.
- \f]
- *
- * and the strain tensor \f$\boldsymbol{\epsilon}\f$ as a function of the solid displacement gradient \f$\textbf{grad} \boldsymbol{u}\f$:
- \f[
- \boldsymbol{\epsilon} = \frac{1}{2} \, (\textbf{grad} \boldsymbol{u} + \textbf{grad}^T \boldsymbol{u}).
- \f]
- *
- * Here, the rock mechanics sign convention is switch off which means compressive stresses are < 0 and tensile stresses are > 0.
- * The rock mechanics sign convention can be switched on for the vtk output via the property system.
- *
- * The effective porosity is calculated as a function of the solid displacement:
- \f[
- \phi_{eff} = \frac{\phi_{init} + \text{div} \boldsymbol{u}}{1 + \text{div}}
- \f]
- * All equations are discretized using a vertex-centered finite volume (box)
- * or cell-centered finite volume scheme as spatial
- * and the implicit Euler method as time discretization.
- *
- * The primary variables are the pressure \f$p\f$ and the mole or mass fraction of dissolved component \f$x\f$ and the solid
- * displacement vector \f$\boldsymbol{u}\f$.
- */
-
-
-template
-class ElOnePTwoCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
-{
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- enum {
- dim = GridView::dimension
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldVector DimVector;
- typedef Dune::FieldMatrix DimMatrix;
-
-public:
- /*!
- * \brief \copybrief ImplicitModel::addOutputVtkFields
- *
- * Specialization for the ElOnePTwoCBoxModel, add one-phase two-component
- * properties, solid displacement, stresses, effective properties and the
- * process rank to the VTK writer.
- */
- template
- void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer) {
-
- // check whether compressive stresses are defined to be positive
- // (rockMechanicsSignConvention_ == true) or negative
- rockMechanicsSignConvention_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, RockMechanicsSignConvention);
-
- typedef Dune::BlockVector > ScalarField;
- typedef Dune::BlockVector > VectorField;
-
- // create the required scalar and vector fields
- unsigned numVertices = this->gridView_().size(dim);
- unsigned numElements = this->gridView_().size(0);
-
- // create the required fields for vertex data
- ScalarField &pressure = *writer.allocateManagedBuffer(numVertices);
- ScalarField &moleFraction0 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &moleFraction1 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &massFraction0 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &massFraction1 = *writer.allocateManagedBuffer(numVertices);
- VectorField &displacement = *writer.template allocateManagedBuffer(numVertices);
- ScalarField &density = *writer.allocateManagedBuffer(numVertices);
- ScalarField &viscosity = *writer.allocateManagedBuffer(numVertices);
- ScalarField &porosity = *writer.allocateManagedBuffer(numVertices);
- ScalarField &Kx = *writer.allocateManagedBuffer(numVertices);
-
- // create the required fields for element data
- // effective stresses
- VectorField &effStressX = *writer.template allocateManagedBuffer(numElements);
- VectorField &effStressY = *writer.template allocateManagedBuffer(numElements);
- VectorField &effStressZ = *writer.template allocateManagedBuffer(numElements);
- // total stresses
- VectorField &totalStressX = *writer.template allocateManagedBuffer<
- Scalar, dim>(numElements);
- VectorField &totalStressY = *writer.template allocateManagedBuffer<
- Scalar, dim>(numElements);
- VectorField &totalStressZ = *writer.template allocateManagedBuffer<
- Scalar, dim>(numElements);
-
- // principal stresses
- ScalarField &principalStress1 = *writer.allocateManagedBuffer(
- numElements);
- ScalarField &principalStress2 = *writer.allocateManagedBuffer(
- numElements);
- ScalarField &principalStress3 = *writer.allocateManagedBuffer(
- numElements);
-
- ScalarField &effPorosity = *writer.allocateManagedBuffer(numElements);
- ScalarField &cellPorosity = *writer.allocateManagedBuffer(numElements);
- ScalarField &cellKx = *writer.allocateManagedBuffer(numElements);
- ScalarField &cellPressure = *writer.allocateManagedBuffer(numElements);
-
- // initialize cell stresses, cell-wise hydraulic parameters and cell pressure with zero
- for (unsigned int eIdx = 0; eIdx < numElements; ++eIdx) {
- effStressX[eIdx] = Scalar(0.0);
- if (dim >= 2)
- effStressY[eIdx] = Scalar(0.0);
- if (dim >= 3)
- effStressZ[eIdx] = Scalar(0.0);
-
- totalStressX[eIdx] = Scalar(0.0);
- if (dim >= 2)
- totalStressY[eIdx] = Scalar(0.0);
- if (dim >= 3)
- totalStressZ[eIdx] = Scalar(0.0);
-
- principalStress1[eIdx] = Scalar(0.0);
- if (dim >= 2)
- principalStress2[eIdx] = Scalar(0.0);
- if (dim >= 3)
- principalStress3[eIdx] = Scalar(0.0);
-
- effPorosity[eIdx] = Scalar(0.0);
- cellPorosity[eIdx] = Scalar(0.0);
- cellKx[eIdx] = Scalar(0.0);
- cellPressure[eIdx] = Scalar(0.0);
- }
- ScalarField &rank = *writer.allocateManagedBuffer(numElements);
-
-
- FVElementGeometry fvGeometry;
- ElementVolumeVariables elemVolVars;
- ElementBoundaryTypes elemBcTypes;
-
- // initialize start and end of element iterator
- // loop over all elements (cells)
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- unsigned int eIdx = this->problem_().model().elementMapper().index(element);
- rank[eIdx] = this->gridView_().comm().rank();
-
- fvGeometry.update(this->gridView_(), element);
- elemBcTypes.update(this->problem_(), element, fvGeometry);
- elemVolVars.update(this->problem_(), element, fvGeometry, false);
-
- // loop over all local vertices of the cell
- int numScv = element.subEntities(dim);
-
- for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
- {
- unsigned int vIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dim);
-
- pressure[vIdxGlobal] = elemVolVars[scvIdx].pressure();
- moleFraction0[vIdxGlobal] = elemVolVars[scvIdx].moleFraction(0);
- moleFraction1[vIdxGlobal] = elemVolVars[scvIdx].moleFraction(1);
- massFraction0[vIdxGlobal] = elemVolVars[scvIdx].massFraction(0);
- massFraction1[vIdxGlobal] = elemVolVars[scvIdx].massFraction(1);
- // in case of rock mechanics sign convention solid displacement is
- // defined to be negative if it points in positive coordinate direction
- if(rockMechanicsSignConvention_){
- DimVector tmpDispl;
- tmpDispl = Scalar(0);
- tmpDispl -= elemVolVars[scvIdx].displacement();
- displacement[vIdxGlobal] = tmpDispl;
- }
-
- else
- displacement[vIdxGlobal] = elemVolVars[scvIdx].displacement();
-
- density[vIdxGlobal] = elemVolVars[scvIdx].density();
- viscosity[vIdxGlobal] = elemVolVars[scvIdx].viscosity();
- porosity[vIdxGlobal] = elemVolVars[scvIdx].porosity();
- Kx[vIdxGlobal] = this->problem_().spatialParams().intrinsicPermeability(
- element, fvGeometry, scvIdx)[0][0];
- // calculate cell quantities by adding up scv quantities and dividing through numScv
- cellPorosity[eIdx] += elemVolVars[scvIdx].porosity() / numScv;
- cellKx[eIdx] += this->problem_().spatialParams().intrinsicPermeability(
- element, fvGeometry, scvIdx)[0][0] / numScv;
- cellPressure[eIdx] += elemVolVars[scvIdx].pressure() / numScv;
- };
-
- // calculate cell quantities for variables which are defined at the integration point
- Scalar tmpEffPoro;
- DimMatrix tmpEffStress;
- tmpEffStress = Scalar(0);
- tmpEffPoro = Scalar(0);
-
- // loop over all scv-faces of the cell
- for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++) {
-
- //prepare the flux calculations (set up and prepare geometry, FE gradients)
- FluxVariables fluxVars;
- fluxVars.update(this->problem_(),
- element, fvGeometry,
- fIdx,
- elemVolVars);
-
- // divide by number of scv-faces and sum up edge values
- tmpEffPoro = fluxVars.effPorosity() / fvGeometry.numScvf;
- tmpEffStress = fluxVars.sigma();
- tmpEffStress /= fvGeometry.numScvf;
-
- effPorosity[eIdx] += tmpEffPoro;
-
- // in case of rock mechanics sign convention compressive stresses
- // are defined to be positive
- if(rockMechanicsSignConvention_){
- effStressX[eIdx] -= tmpEffStress[0];
- if (dim >= 2) {
- effStressY[eIdx] -= tmpEffStress[1];
- }
- if (dim >= 3) {
- effStressZ[eIdx] -= tmpEffStress[2];
- }
- }
- else{
- effStressX[eIdx] += tmpEffStress[0];
- if (dim >= 2) {
- effStressY[eIdx] += tmpEffStress[1];
- }
- if (dim >= 3) {
- effStressZ[eIdx] += tmpEffStress[2];
- }
- }
- }
-
- // calculate total stresses
- // in case of rock mechanics sign convention compressive stresses
- // are defined to be positive and total stress is calculated by adding the pore pressure
- if(rockMechanicsSignConvention_){
- totalStressX[eIdx][0] = effStressX[eIdx][0] + cellPressure[eIdx];
- totalStressX[eIdx][1] = effStressX[eIdx][1];
- totalStressX[eIdx][2] = effStressX[eIdx][2];
- if (dim >= 2) {
- totalStressY[eIdx][0] = effStressY[eIdx][0];
- totalStressY[eIdx][1] = effStressY[eIdx][1] + cellPressure[eIdx];
- totalStressY[eIdx][2] = effStressY[eIdx][2];
- }
- if (dim >= 3) {
- totalStressZ[eIdx][0] = effStressZ[eIdx][0];
- totalStressZ[eIdx][1] = effStressZ[eIdx][1];
- totalStressZ[eIdx][2] = effStressZ[eIdx][2] + cellPressure[eIdx];
- }
- }
- else{
- totalStressX[eIdx][0] = effStressX[eIdx][0] - cellPressure[eIdx];
- totalStressX[eIdx][1] = effStressX[eIdx][1];
- totalStressX[eIdx][2] = effStressX[eIdx][2];
- if (dim >= 2) {
- totalStressY[eIdx][0] = effStressY[eIdx][0];
- totalStressY[eIdx][1] = effStressY[eIdx][1] - cellPressure[eIdx];
- totalStressY[eIdx][2] = effStressY[eIdx][2];
- }
- if (dim >= 3) {
- totalStressZ[eIdx][0] = effStressZ[eIdx][0];
- totalStressZ[eIdx][1] = effStressZ[eIdx][1];
- totalStressZ[eIdx][2] = effStressZ[eIdx][2] - cellPressure[eIdx];
- }
- }
- }
-
- // calculate principal stresses i.e. the eigenvalues of the total stress tensor
- Scalar a1, a2, a3;
- DimMatrix totalStress;
- DimVector eigenValues;
- a1=Scalar(0);
- a2=Scalar(0);
- a3=Scalar(0);
-
- for (unsigned int eIdx = 0; eIdx < numElements; eIdx++)
- {
- eigenValues = Scalar(0);
- totalStress = Scalar(0);
-
- totalStress[0] = totalStressX[eIdx];
- if (dim >= 2)
- totalStress[1] = totalStressY[eIdx];
- if (dim >= 3)
- totalStress[2] = totalStressZ[eIdx];
-
- calculateEigenValues(eigenValues, totalStress);
-
-
- for (int i = 0; i < dim; i++)
- {
- if (std::isnan(eigenValues[i]))
- eigenValues[i] = 0.0;
- }
-
- // sort principal stresses: principalStress1 >= principalStress2 >= principalStress3
- if (dim == 2) {
- a1 = eigenValues[0];
- a2 = eigenValues[1];
-
- if (a1 >= a2) {
- principalStress1[eIdx] = a1;
- principalStress2[eIdx] = a2;
- } else {
- principalStress1[eIdx] = a2;
- principalStress2[eIdx] = a1;
- }
- }
-
- if (dim == 3) {
- a1 = eigenValues[0];
- a2 = eigenValues[1];
- a3 = eigenValues[2];
-
- if (a1 >= a2) {
- if (a1 >= a3) {
- principalStress1[eIdx] = a1;
- if (a2 >= a3) {
- principalStress2[eIdx] = a2;
- principalStress3[eIdx] = a3;
- }
- else //a3 > a2
- {
- principalStress2[eIdx] = a3;
- principalStress3[eIdx] = a2;
- }
- }
- else // a3 > a1
- {
- principalStress1[eIdx] = a3;
- principalStress2[eIdx] = a1;
- principalStress3[eIdx] = a2;
- }
- } else // a2>a1
- {
- if (a2 >= a3) {
- principalStress1[eIdx] = a2;
- if (a1 >= a3) {
- principalStress2[eIdx] = a1;
- principalStress3[eIdx] = a3;
- }
- else //a3>a1
- {
- principalStress2[eIdx] = a3;
- principalStress3[eIdx] = a1;
- }
- }
- else //a3>a2
- {
- principalStress1[eIdx] = a3;
- principalStress2[eIdx] = a2;
- principalStress3[eIdx] = a1;
- }
- }
- }
-
- }
-
- writer.attachVertexData(pressure, "P");
-
- char nameMoleFraction0[42], nameMoleFraction1[42];
- snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
- snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
- writer.attachVertexData(moleFraction0, nameMoleFraction0);
- writer.attachVertexData(moleFraction1, nameMoleFraction1);
-
- char nameMassFraction0[42], nameMassFraction1[42];
- snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
- snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
- writer.attachVertexData(massFraction0, nameMassFraction0);
- writer.attachVertexData(massFraction1, nameMassFraction1);
-
- writer.attachVertexData(displacement, "u", dim);
- writer.attachVertexData(density, "rho");
- writer.attachVertexData(viscosity, "mu");
- writer.attachVertexData(porosity, "porosity");
- writer.attachVertexData(Kx, "Kx");
- writer.attachCellData(cellPorosity, "porosity");
- writer.attachCellData(cellKx, "Kx");
- writer.attachCellData(effPorosity, "effective porosity");
-
- writer.attachCellData(totalStressX, "total stresses X", dim);
- if (dim >= 2)
- writer.attachCellData(totalStressY, "total stresses Y", dim);
- if (dim >= 3)
- writer.attachCellData(totalStressZ, "total stresses Z", dim);
-
- writer.attachCellData(effStressX, "effective stress changes X", dim);
- if (dim >= 2)
- writer.attachCellData(effStressY, "effective stress changes Y", dim);
- if (dim >= 3)
- writer.attachCellData(effStressZ, "effective stress changes Z", dim);
-
- writer.attachCellData(principalStress1, "principal stress 1");
- if (dim >= 2)
- writer.attachCellData(principalStress2, "principal stress 2");
- if (dim >= 3)
- writer.attachCellData(principalStress3, "principal stress 3");
-
- writer.attachCellData(cellPressure, "P");
-
- writer.attachCellData(rank, "rank");
-
- }
-private:
- bool rockMechanicsSignConvention_;
-
-};
-}
-#include "propertydefaults.hh"
-#endif
diff --git a/dumux/geomechanics/el1p2c/properties.hh b/dumux/geomechanics/el1p2c/properties.hh
deleted file mode 100644
index 66970504124e9932b7cca2a646b6cde3e98de288..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/properties.hh
+++ /dev/null
@@ -1,62 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup ElOnePTwoCBoxModel
- * \file
- *
- * \brief Defines the properties required for the one-phase two-component
- * linear elasticity model.
- *
- * This class inherits from the properties of the one-phase two-component model and
- * from the properties of the linear elasticity model
- */
-
-#ifndef DUMUX_ELASTIC1P2C_PROPERTIES_HH
-#define DUMUX_ELASTIC1P2C_PROPERTIES_HH
-
-#include
-#include
-
-
-namespace Dumux
-{
-// \{
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for the single-phase, two-component linear elasticity problems
-NEW_TYPE_TAG(BoxElasticOnePTwoC, INHERITS_FROM(BoxModel));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-//! Returns whether the stabilization terms are included in the balance equations
-NEW_PROP_TAG(ImplicitWithStabilization);
-//! Returns whether the output should be written according to rock mechanics sign convention (compressive stresses > 0)
-NEW_PROP_TAG(VtkRockMechanicsSignConvention);
-}
-// \}
-}
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/propertydefaults.hh b/dumux/geomechanics/el1p2c/propertydefaults.hh
deleted file mode 100644
index 373a2fe316b3bfa33d6c36ffccb6d2d80f5f5de8..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/propertydefaults.hh
+++ /dev/null
@@ -1,128 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup ElOnePTwoCBoxModel
- * \file
- *
- * \brief Defines the properties required for the one-phase two-component
- * linear-elastic model.
- *
- * This class inherits from the properties of the one-phase two-component model and
- * from the properties of the simple linear-elastic model
- */
-
-#ifndef DUMUX_ELASTIC1P2C_PROPERTY_DEFAULTS_HH
-#define DUMUX_ELASTIC1P2C_PROPERTY_DEFAULTS_HH
-
-#include "properties.hh"
-
-#include "model.hh"
-#include "localresidual.hh"
-#include "localjacobian.hh"
-#include "fluxvariables.hh"
-#include "elementvolumevariables.hh"
-#include "volumevariables.hh"
-#include "indices.hh"
-#include
-#include
-
-
-namespace Dumux
-{
-// \{
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Property defaults
-//////////////////////////////////////////////////////////////////
-//!< set the number of equations to the space dimension of the problem
-SET_PROP(BoxElasticOnePTwoC, NumEq)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- enum{dim = GridView::dimension};
-public:
- static const int value = dim + 2;
-};
-
-SET_INT_PROP(BoxElasticOnePTwoC, NumPhases, 1); //!< The number of phases in the 1p2c model is 1
-SET_INT_PROP(BoxElasticOnePTwoC, NumComponents, 2); //!< The number of components in the 1p2c model is 2
-
-//! Use the linear elasticity local residual function for the elasticity model
-SET_TYPE_PROP(BoxElasticOnePTwoC,
- LocalResidual,
- ElOnePTwoCLocalResidual);
-
-//! Use the linear elasticity local residual function for the elasticity model
-SET_TYPE_PROP(BoxElasticOnePTwoC,
- LocalJacobian,
- ElOnePTwoCLocalJacobian);
-
-//! define the model
-SET_TYPE_PROP(BoxElasticOnePTwoC, Model, ElOnePTwoCModel);
-
-//! define the ElementVolumeVariables
-SET_TYPE_PROP(BoxElasticOnePTwoC, ElementVolumeVariables, ElOnePTwoCElementVolumeVariables);
-
-//! define the VolumeVariables
-SET_TYPE_PROP(BoxElasticOnePTwoC, VolumeVariables, ElOnePTwoCVolumeVariables);
-
-//! define the FluxVariables
-SET_TYPE_PROP(BoxElasticOnePTwoC, FluxVariables, ElOnePTwoCFluxVariables);
-
-//! Set the indices used by the linear elasticity model
-SET_TYPE_PROP(BoxElasticOnePTwoC, Indices, ElOnePTwoCIndices);
-
-//! Set the phaseIndex per default to zero (important for two-phase fluidsystems).
-SET_INT_PROP(BoxElasticOnePTwoC, PhaseIdx, 0);
-
-SET_PROP(BoxElasticOnePTwoC, FluidState){
- private:
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- public:
- typedef CompositionalFluidState type;
-};
-
-//! set default upwind weights to 1.0, i.e. fully upwind
-SET_SCALAR_PROP(BoxElasticOnePTwoC, ImplicitMassUpwindWeight, 1.0);
-SET_SCALAR_PROP(BoxElasticOnePTwoC, ImplicitMobilityUpwindWeight, 1.0);
-
-// enable gravity by default
-SET_BOOL_PROP(BoxElasticOnePTwoC, ProblemEnableGravity, true);
-
-// enable gravity by default
-SET_BOOL_PROP(BoxElasticOnePTwoC, ImplicitWithStabilization, true);
-
-//! The model after Millington (1961) is used for the effective diffusivity
-SET_PROP(BoxElasticOnePTwoC, EffectiveDiffusivityModel)
-{ private :
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- public:
- typedef DiffusivityMillingtonQuirk type;
-};
-
-// write the stress and displacement output according to rock mechanics sign convention (compressive stresses > 0)
-SET_BOOL_PROP(BoxElasticOnePTwoC, VtkRockMechanicsSignConvention, false);
-}
-}
-
-#endif
diff --git a/dumux/geomechanics/el1p2c/volumevariables.hh b/dumux/geomechanics/el1p2c/volumevariables.hh
deleted file mode 100644
index f3773ece1077e80e287c82abb0b639c712e7a36c..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el1p2c/volumevariables.hh
+++ /dev/null
@@ -1,205 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Quantities required by the single-phase, two-component
- * linear elasticity model which are defined on a vertex.
- */
-#ifndef DUMUX_ELASTIC1P2C_VOLUME_VARIABLES_HH
-#define DUMUX_ELASTIC1P2C_VOLUME_VARIABLES_HH
-
-
-#include
-#include
-
-#include "properties.hh"
-
-namespace Dumux {
-/*!
- * \ingroup ElOnePTwoCBoxModel
- * \ingroup ImplicitVolumeVariables
- * \brief Contains the quantities which are constant within a
- * finite volume in the single-phase, two-component, linear elasticity model.
- *
- * This class inherits from the volumevariables of the one-phase
- * two-component model
- */
-template
-class ElOnePTwoCVolumeVariables : public OnePTwoCVolumeVariables{
-
- typedef OnePTwoCVolumeVariables ParentType;
-
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- enum { dim = GridView::dimension };
-
- enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx) };
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldVector DimVector;
-
-public:
- /*!
- * \copydoc ImplicitVolumeVariables::update
- */
- void update(const PrimaryVariables &priVars,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx,
- bool isOldSol)
- {
-
- ParentType::update(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
- int vIdxGlobal = problem.vertexMapper().subIndex(element, scvIdx, dim);
-
- primaryVars_ = priVars;
- prevPrimaryVars_ = problem.model().prevSol()[vIdxGlobal];
-
- ParentType prev1p2cVolVars;
- prev1p2cVolVars.update(problem.model().prevSol()[vIdxGlobal],
- problem,
- element,
- fvGeometry,
- scvIdx,
- true);
-
- dPressure_ = this->pressure() - prev1p2cVolVars.pressure();
-
- for (int i = 0; i < dim; ++i){
- displacement_[i] = primaryVars_[Indices::u(i)];
- prevDisplacement_[i] = prevPrimaryVars_[Indices::u(i)];
- }
-
- dU_ = displacement_ - prevDisplacement_;
-
- const Dune::FieldVector &lameParams =
- problem.spatialParams().lameParams(element, fvGeometry, scvIdx);
-
- lambda_ = lameParams[0];
- mu_ = lameParams[1];
-
- rockDensity_ = problem.spatialParams().rockDensity(element, scvIdx);
- }
-
- /*!
- * \brief Return the vector of primary variables
- */
- const PrimaryVariables &primaryVars() const
- { return primaryVars_; }
-
- /*!
- * \brief Sets the evaluation point used in the by the local jacobian.
- */
- void setEvalPoint(const Implementation *ep)
- { }
-
- /*!
- * \brief Return the Lame parameter lambda \f$\mathrm{[Pa]}\f$ within the control volume.
- */
- Scalar lambda() const
- { return lambda_; }
-
- /*!
- * \brief Return the Lame parameter mu \f$\mathrm{[Pa]}\f$ within the control volume.
- */
- Scalar mu() const
- { return mu_; }
-
- /*!
- * \brief Returns the rock density \f$\mathrm{[kg / m^3]}\f$ within the control volume .
- */
- Scalar rockDensity() const
- { return rockDensity_; }
-
- /*!
- * \brief Returns the change in solid displacement \f$\mathrm{[m]}\f$ between
- * the last and the current timestep for the space direction dimIdx within the control volume.
- */
- Scalar dU(int dimIdx) const
- { return dU_[dimIdx]; }
-
- /*!
- * \brief Returns the change in pressure \f$\mathrm{[Pa]}\f$ between the last and the
- * current timestep within the control volume.
- */
- Scalar dPressure() const
- { return dPressure_; }
-
- /*!
- * \brief Returns the solid displacement \f$\mathrm{[m]}\f$ in space
- * directions dimIdx within the control volume.
- */
- Scalar displacement(int dimIdx) const
- { return displacement_[dimIdx]; }
-
- /*!
- * \brief Returns the solid displacement vector \f$\mathrm{[m]}\f$
- * within the control volume.
- */
- const DimVector &displacement() const
- { return displacement_; }
-
-
- /*!
- * \brief the effective porosity and volumetric strain divU is defined as mutable variable here since it
- * is updated within the elementVolumeVariables.
- */
- mutable Scalar effPorosity;
- mutable Scalar divU;
-
- /*!
- * \brief Returns the mass fraction of a given component in the
- * given fluid phase within the control volume.
- *
- * \param compIdx The component index
- */
- Scalar massFraction(const int compIdx) const
- { return this->fluidState_.massFraction(phaseIdx, compIdx); }
-
- /*!
- * \brief Returns the mole fraction of a given component in the
- * given fluid phase within the control volume.
- *
- * \param compIdx The component index
- */
- Scalar moleFraction(const int compIdx) const
- { return this->fluidState_.moleFraction(phaseIdx, compIdx); }
-
-protected:
- PrimaryVariables primaryVars_, prevPrimaryVars_;
- DimVector displacement_, prevDisplacement_;
- DimVector dU_;
- Scalar dPressure_;
- Scalar lambda_;
- Scalar mu_;
- Scalar rockDensity_;
-};
-
-}
-
-#endif
diff --git a/dumux/geomechanics/el2p/CMakeLists.txt b/dumux/geomechanics/el2p/CMakeLists.txt
deleted file mode 100644
index ab8c45b2f6cef4a4c01bb847084eb43e93e82bd9..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/CMakeLists.txt
+++ /dev/null
@@ -1,18 +0,0 @@
-
-#install headers
-install(FILES
-amgbackend.hh
-assembler.hh
-basemodel.hh
-elementvolumevariables.hh
-fluxvariables.hh
-indices.hh
-localjacobian.hh
-localoperator.hh
-localresidual.hh
-model.hh
-newtoncontroller.hh
-properties.hh
-propertydefaults.hh
-volumevariables.hh
-DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/geomechanics/el2p)
diff --git a/dumux/geomechanics/el2p/amgbackend.hh b/dumux/geomechanics/el2p/amgbackend.hh
deleted file mode 100644
index bef3c01da03396c5aef81cf7d576ac389b96666b..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/amgbackend.hh
+++ /dev/null
@@ -1,235 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Linear
- *
- * \brief Wraps the AMG backend such that it can be used for the el2p model.
- */
-#ifndef DUMUX_EL2P_AMGBACKEND_HH
-#define DUMUX_EL2P_AMGBACKEND_HH
-
-#include
-
-namespace Dumux {
-
-/*!
- * \brief Base class for the ElTwoP AMGBackend.
- */
-template
-class El2PAMGBackendBase : public AMGBackend
-{
- typedef AMGBackend ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-
-public:
- /*!
- * \copydoc AMGBackend::AMGBackend()
- */
- El2PAMGBackendBase(const Problem& problem)
- : ParentType(problem)
- {}
-};
-
-/*!
- * \brief Specialization for the parallel setting.
- */
-template
-class El2PAMGBackendBase : public AMGBackend
-{
- typedef AMGBackend ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
- typedef typename Dune::FieldMatrix MatrixBlock;
- typedef typename Dune::BCRSMatrix BlockMatrix;
- typedef typename Dune::FieldVector VectorBlock;
- typedef typename Dune::BlockVector BlockVector;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- enum { dim = GridView::dimension };
-
-public:
- /*!
- * \copydoc AMGBackend::AMGBackend()
- */
- El2PAMGBackendBase(const Problem& problem)
- : ParentType(problem)
- {
- createBlockMatrixAndVectors_();
- }
-
- /*!
- * \copydoc AMGBackend::solve()
- */
- template
- bool solve(Matrix& A, Vector& x, Vector& b)
- {
- flatToBlocked_(A, x, b);
- int converged = ParentType::solve(*aBlocked_, *xBlocked_, *bBlocked_);
- blockedToFlat_(x, b);
- return converged;
- }
-
-private:
- void createBlockMatrixAndVectors_()
- {
- int numVertices = this->problem().gridView().size(dim);
-
- aBlocked_ = std::make_shared(numVertices, numVertices, BlockMatrix::random);
- xBlocked_ = std::make_shared(numVertices);
- bBlocked_ = std::make_shared(numVertices);
-
- // find out the global indices of the neighboring vertices of
- // each vertex
- typedef std::set NeighborSet;
- std::vector neighbors(numVertices);
- for (const auto& element : elements(this->problem().gridView())) {
-
- // loop over all element vertices
- int n = element.subEntities(dim);
- for (int i = 0; i < n - 1; ++i) {
- int globalI = this->problem().vertexMapper().subIndex(element, i, dim);
- for (int j = i + 1; j < n; ++j) {
- int globalJ = this->problem().vertexMapper().subIndex(element, j, dim);
- // make sure that vertex j is in the neighbor set
- // of vertex i and vice-versa
- neighbors[globalI].insert(globalJ);
- neighbors[globalJ].insert(globalI);
- }
- }
- }
-
- // make vertices neighbors to themselfs
- for (int i = 0; i < numVertices; ++i)
- neighbors[i].insert(i);
-
- // allocate space for the rows of the matrix
- for (int i = 0; i < numVertices; ++i) {
- aBlocked_->setrowsize(i, neighbors[i].size());
- }
- aBlocked_->endrowsizes();
-
- // fill the rows with indices. each vertex talks to all of its
- // neighbors. (it also talks to itself since vertices are
- // sometimes quite egocentric.)
- for (int i = 0; i < numVertices; ++i) {
- auto nIt = neighbors[i].begin();
- const auto& nEndIt = neighbors[i].end();
- for (; nIt != nEndIt; ++nIt) {
- aBlocked_->addindex(i, *nIt);
- }
- }
- aBlocked_->endindices();
- }
-
- template
- void flatToBlocked_(const FlatMatrix& aFlat,
- const FlatVector& xFlat,
- const FlatVector& bFlat)
- {
- unsigned numBlocks = xBlocked_->size();
- static const unsigned numMassEq = numEq - dim;
- for (unsigned rowBlockIdx = 0; rowBlockIdx < numBlocks; ++rowBlockIdx)
- {
- for (unsigned rowEqIdx = 0; rowEqIdx < numEq; ++rowEqIdx)
- {
- unsigned rowFlatIdx;
- if (rowEqIdx < numMassEq)
- rowFlatIdx = rowBlockIdx*numMassEq + rowEqIdx;
- else
- rowFlatIdx = numBlocks*numMassEq + rowBlockIdx*dim + rowEqIdx - numMassEq;
-
- (*xBlocked_)[rowBlockIdx][rowEqIdx] = xFlat[rowFlatIdx];
- (*bBlocked_)[rowBlockIdx][rowEqIdx] = bFlat[rowFlatIdx];
-
- for (auto colBlockIt = (*aBlocked_)[rowBlockIdx].begin();
- colBlockIt != (*aBlocked_)[rowBlockIdx].end(); ++colBlockIt)
- {
- unsigned colBlockIdx = colBlockIt.index();
- auto& aBlock = (*aBlocked_)[rowBlockIdx][colBlockIdx];
-
- for (unsigned colEqIdx = 0; colEqIdx < numEq; ++colEqIdx)
- {
- unsigned colFlatIdx;
- if (colEqIdx < numMassEq)
- colFlatIdx = colBlockIdx*numMassEq + colEqIdx;
- else
- colFlatIdx = numBlocks*numMassEq + colBlockIdx*dim + colEqIdx - numMassEq;
-
- aBlock[rowEqIdx][colEqIdx] = aFlat[rowFlatIdx][colFlatIdx];
- }
- }
- }
- }
- }
-
- template
- void blockedToFlat_(FlatVector& xFlat,
- FlatVector& bFlat)
- {
- unsigned numBlocks = xBlocked_->size();
- static const unsigned numMassEq = numEq - dim;
- for (unsigned rowBlockIdx = 0; rowBlockIdx < numBlocks; ++rowBlockIdx)
- {
- for (unsigned rowEqIdx = 0; rowEqIdx < numEq; ++rowEqIdx)
- {
- unsigned rowFlatIdx;
- if (rowEqIdx < numMassEq)
- rowFlatIdx = rowBlockIdx*numMassEq + rowEqIdx;
- else
- rowFlatIdx = numBlocks*numMassEq + rowBlockIdx*dim + rowEqIdx - numMassEq;
-
- xFlat[rowFlatIdx] = (*xBlocked_)[rowBlockIdx][rowEqIdx];
- bFlat[rowFlatIdx] = (*bBlocked_)[rowBlockIdx][rowEqIdx];
- }
- }
- }
-
- std::shared_ptr aBlocked_;
- std::shared_ptr xBlocked_;
- std::shared_ptr bBlocked_;
-};
-
-/*!
- * \brief Wraps the AMG backend such that it can be used for the el2p model.
- */
-template
-class El2PAMGBackend : public El2PAMGBackendBase<
- TypeTag,
- Dune::Capabilities::canCommunicate::v>
-{
- typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
- enum { dofCodim = Grid::dimension };
- enum { isParallel = Dune::Capabilities::canCommunicate::v };
- typedef El2PAMGBackendBase ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-
-public:
- /*!
- * \copydoc AMGBackend::AMGBackend()
- */
- El2PAMGBackend(const Problem& problem)
- : ParentType(problem)
- {}
-};
-
-} // namespace Dumux
-
-#endif // DUMUX_EL2P_AMGBACKEND_HH
diff --git a/dumux/geomechanics/el2p/assembler.hh b/dumux/geomechanics/el2p/assembler.hh
deleted file mode 100644
index f995b84864d778cc6c742468dc2e5d07ec6d2b8c..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/assembler.hh
+++ /dev/null
@@ -1,605 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file contains an assembler for the Jacobian matrix
- * of the two-phase linear-elastic model based on PDELab.
- */
-#ifndef DUMUX_EL2P_ASSEMBLER_HH
-#define DUMUX_EL2P_ASSEMBLER_HH
-
-#include "properties.hh"
-#include "localoperator.hh"
-
-namespace Dumux {
-
-namespace Properties
-{
-NEW_PROP_TAG(PressureFEM); //!< Finite element space used for pressure, saturation, ...
-NEW_PROP_TAG(DisplacementFEM); //!< Finite element space used for displacement
-NEW_PROP_TAG(PressureGridFunctionSpace); //!< Grid function space used for pressure, saturation, ...
-NEW_PROP_TAG(DisplacementGridFunctionSpace); //!< Grid function space used for displacement
-}
-
-namespace PDELab {
-
-/*!
- * \brief An assembler for the Jacobian matrix
- * of the two-phase linear-elastic model based on PDELab.
- */
-template
-class El2PAssembler
-{
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
-
- typedef typename GET_PROP_TYPE(TypeTag, PressureFEM) PressureFEM;
- typedef typename GET_PROP_TYPE(TypeTag, PressureGridFunctionSpace) PressureGFS;
- typedef typename PressureGFS::template Child<0>::Type PressureScalarGFS;
-
- typedef typename GET_PROP_TYPE(TypeTag, DisplacementFEM) DisplacementFEM;
- typedef typename GET_PROP_TYPE(TypeTag, DisplacementGridFunctionSpace) DisplacementGFS;
- typedef typename DisplacementGFS::template Child<0>::Type DisplacementScalarGFS;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
- typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
- typedef typename GET_PROP_TYPE(TypeTag, ConstraintsTrafo) ConstraintsTrafo;
- typedef typename GET_PROP_TYPE(TypeTag, LocalOperator) LocalOperator;
- typedef typename GET_PROP_TYPE(TypeTag, GridOperator) GridOperator;
-
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
- enum{dim = GridView::dimension};
- typedef typename GridView::template Codim<0>::Entity Element;
-
- typedef typename GridView::template Codim::Entity Vertex;
-
- enum {
- enablePartialReassemble = GET_PROP_VALUE(TypeTag, ImplicitEnablePartialReassemble),
- enableJacobianRecycling = GET_PROP_VALUE(TypeTag, ImplicitEnableJacobianRecycling),
- };
-
- // copying the jacobian assembler is not a good idea
- El2PAssembler(const El2PAssembler &);
-
-public:
- /*!
- * \brief The colors of elements and vertices required for partial
- * Jacobian reassembly.
- */
- enum EntityColor {
- /*!
- * Vertex/element that needs to be reassembled because some
- * relative error is above the tolerance
- */
- Red,
-
- /*!
- * Vertex/element that needs to be reassembled because a
- * neighboring element/vertex is red
- */
- Yellow,
-
- /*!
- * Yellow vertex has only non-green neighbor elements.
- *
- * This means that its relative error is below the tolerance,
- * but its defect can be linearized without any additional
- * cost. This is just an "internal" color which is not used
- * ouside of the jacobian assembler.
- */
- Orange,
-
- /*!
- * Vertex/element that does not need to be reassembled
- */
- Green
- };
-
- El2PAssembler()
- {
- // set reassemble tolerance to 0, so that if partial
- // reassembly of the jacobian matrix is disabled, the
- // reassemble tolerance is always smaller than the current
- // relative tolerance
- reassembleTolerance_ = 0.0;
- }
-
- /*!
- * \brief Initialize the jacobian assembler.
- *
- * At this point we can assume that all objects in the problem and
- * the model have been allocated. We can not assume that they are
- * fully initialized, though.
- *
- * \param problem The problem object
- */
- void init(Problem& problem)
- {
- problemPtr_ = &problem;
-
- constraints_ = std::make_shared();
-
- pressureFEM_ = std::make_shared(problemPtr_->gridView());
- pressureScalarGFS_ = std::make_shared(problemPtr_->gridView(), *pressureFEM_, *constraints_);
- pressureGFS_ = std::make_shared(*pressureScalarGFS_);
-
- displacementFEM_ = std::make_shared(problemPtr_->gridView());
- displacementScalarGFS_ = std::make_shared(problemPtr_->gridView(), *displacementFEM_, *constraints_);
- displacementGFS_ = std::make_shared(*displacementScalarGFS_);
-
- gridFunctionSpace_ = std::make_shared(*pressureGFS_, *displacementGFS_);
-
- constraintsTrafo_ = std::make_shared();
-
- // initialize the grid operator spaces
- localOperator_ = std::make_shared(problemPtr_->model());
- gridOperator_ =
- std::make_shared(*gridFunctionSpace_, *constraintsTrafo_,
- *gridFunctionSpace_, *constraintsTrafo_, *localOperator_);
-
- // allocate raw matrix
- matrix_ = std::make_shared(*gridOperator_);
-
- // initialize the jacobian matrix and the right hand side
- // vector
- *matrix_ = 0;
- reuseMatrix_ = false;
-
- residual_ = std::make_shared(*gridFunctionSpace_);
-
- int numVertices = gridView_().size(dim);
- int numElements = gridView_().size(0);
-
- totalElems_ = gridView_().comm().sum(numElements);
-
- // initialize data needed for partial reassembly
- if (enablePartialReassemble) {
- vertexColor_.resize(numVertices);
- vertexDelta_.resize(numVertices);
- elementColor_.resize(numElements);
- }
- reassembleAll();
- }
-
- /*!
- * \brief Assemble the local jacobian of the problem.
- *
- * The current state of affairs (esp. the previous and the current
- * solutions) is represented by the model object.
- */
- void assemble()
- {
- *matrix_ = 0;
- gridOperator_->jacobian(problemPtr_->model().curSol(), *matrix_);
-
- *residual_ = 0;
- gridOperator_->residual(problemPtr_->model().curSol(), *residual_);
-
- return;
- }
-
- /*!
- * \brief If Jacobian matrix recycling is enabled, this method
- * specifies whether the next call to assemble() just
- * rescales the storage term or does a full reassembly
- *
- * \param yesno If true, only rescale; else do full Jacobian assembly.
- */
- void setMatrixReuseable(bool yesno = true)
- {
- if (enableJacobianRecycling)
- reuseMatrix_ = yesno;
- }
-
- /*!
- * \brief If partial Jacobian matrix reassembly is enabled, this
- * method causes all elements to be reassembled in the next
- * assemble() call.
- */
- void reassembleAll()
- {
- nextReassembleTolerance_ = 0.0;
-
- if (enablePartialReassemble) {
- std::fill(vertexColor_.begin(),
- vertexColor_.end(),
- Red);
- std::fill(elementColor_.begin(),
- elementColor_.end(),
- Red);
- std::fill(vertexDelta_.begin(),
- vertexDelta_.end(),
- 0.0);
- }
- }
-
- /*!
- * \brief Returns the relative error below which a vertex is
- * considered to be "green" if partial Jacobian reassembly
- * is enabled.
- *
- * This returns the _actual_ relative computed seen by
- * computeColors(), not the tolerance which it was given.
- */
- Scalar reassembleTolerance() const
- { return reassembleTolerance_; }
-
- /*!
- * \brief Update the distance where the non-linear system was
- * originally insistently linearized and the point where it
- * will be linerized the next time.
- *
- * This only has an effect if partial reassemble is enabled.
- */
- void updateDiscrepancy(const SolutionVector &u,
- const SolutionVector &uDelta)
- {
- if (!enablePartialReassemble)
- return;
-
- // update the vector with the distances of the current
- // evaluation point used for linearization from the original
- // evaluation point
- using std::abs;
- for (int i = 0; i < vertexDelta_.size(); ++i) {
- PrimaryVariables uCurrent(u[i]);
- PrimaryVariables uNext(uCurrent);
- uNext -= uDelta[i];
-
- // we need to add the distance the solution was moved for
- // this vertex
- Scalar dist = model_().relativeErrorVertex(i,
- uCurrent,
- uNext);
- vertexDelta_[i] += abs(dist);
- }
-
- }
-
- /*!
- * \brief Determine the colors of vertices and elements for partial
- * reassembly given a relative tolerance.
- *
- * The following approach is used:
- *
- * - Set all vertices and elements to 'green'
- * - Mark all vertices as 'red' which exhibit an relative error above
- * the tolerance
- * - Mark all elements which feature a 'red' vetex as 'red'
- * - Mark all vertices which are not 'red' and are part of a
- * 'red' element as 'yellow'
- * - Mark all elements which are not 'red' and contain a
- * 'yellow' vertex as 'yellow'
- *
- * \param relTol The relative error below which a vertex won't be
- * reassembled. Note that this specifies the
- * worst-case relative error between the last
- * linearization point and the current solution and
- * _not_ the delta vector of the Newton iteration!
- */
- void computeColors(Scalar relTol)
- {
- if (!enablePartialReassemble)
- return;
-
- // mark the red vertices and update the tolerance of the
- // linearization which actually will get achieved
- nextReassembleTolerance_ = 0;
- for (int i = 0; i < vertexColor_.size(); ++i) {
- vertexColor_[i] = Green;
- if (vertexDelta_[i] > relTol) {
- // mark vertex as red if discrepancy is larger than
- // the relative tolerance
- vertexColor_[i] = Red;
- }
- using std::max;
- nextReassembleTolerance_ =
- max(nextReassembleTolerance_, vertexDelta_[i]);
- };
-
- // Mark all red elements
- for (const auto& element : elements(gridView_())) {
- // find out whether the current element features a red
- // vertex
- bool isRed = false;
- int numVertices = element.subEntities(dim);
- for (int i=0; i < numVertices; ++i) {
- int globalI = vertexMapper_().subIndex(element, i, dim);
- if (vertexColor_[globalI] == Red) {
- isRed = true;
- break;
- }
- };
-
- // if yes, the element color is also red, else it is not
- // red, i.e. green for the mean time
- int eIdxGlobal = elementMapper_().index(element);
- if (isRed)
- elementColor_[eIdxGlobal] = Red;
- else
- elementColor_[eIdxGlobal] = Green;
- }
-
- // Mark yellow vertices (as orange for the mean time)
- for (const auto& element : elements(gridView_())) {
- int eIdx = this->elementMapper_().index(element);
- if (elementColor_[eIdx] != Red)
- continue; // non-red elements do not tint vertices
- // yellow!
-
- int numVertices = element.subEntities(dim);
- for (int i=0; i < numVertices; ++i) {
- int globalI = vertexMapper_().subIndex(element, i, dim);
- // if a vertex is already red, don't recolor it to
- // yellow!
- if (vertexColor_[globalI] != Red)
- vertexColor_[globalI] = Orange;
- };
- }
-
- // Mark yellow elements
- for (const auto& element : elements(gridView_())) {
- int eIdx = this->elementMapper_().index(element);
- if (elementColor_[eIdx] == Red) {
- continue; // element is red already!
- }
-
- // check whether the element features a yellow
- // (resp. orange at this point) vertex
- bool isYellow = false;
- int numVertices = element.subEntities(dim);
- for (int i=0; i < numVertices; ++i) {
- int globalI = vertexMapper_().subIndex(element, i, dim);
- if (vertexColor_[globalI] == Orange) {
- isYellow = true;
- break;
- }
- };
-
- if (isYellow)
- elementColor_[eIdx] = Yellow;
- }
-
- // Demote orange vertices to yellow ones if it has at least
- // one green element as a neighbor.
- for (const auto& element : elements(gridView_())) {
- int eIdx = this->elementMapper_().index(element);
- if (elementColor_[eIdx] != Green)
- continue; // yellow and red elements do not make
- // orange vertices yellow!
-
- int numVertices = element.subEntities(dim);
- for (int i=0; i < numVertices; ++i) {
- int globalI = vertexMapper_().subIndex(element, i, dim);
- // if a vertex is orange, recolor it to yellow!
- if (vertexColor_[globalI] == Orange)
- vertexColor_[globalI] = Yellow;
- };
- }
-
- // promote the remaining orange vertices to red
- for (int i=0; i < vertexColor_.size(); ++i) {
- // if a vertex is green or yellow don't do anything!
- if (vertexColor_[i] == Green || vertexColor_[i] == Yellow)
- continue;
-
- // make sure the vertex is red (this is a no-op vertices
- // which are already red!)
- vertexColor_[i] = Red;
-
- // set the error of this vertex to 0 because the system
- // will be consistently linearized at this vertex
- vertexDelta_[i] = 0.0;
- };
- };
-
- /*!
- * \brief Returns the reassemble color of a vertex
- *
- * \param element An element which contains the vertex
- * \param vIdx The local index of the vertex in the element.
- */
- int vertexColor(const Element &element, int vIdx) const
- {
- if (!enablePartialReassemble)
- return Red; // reassemble unconditionally!
-
- int vIdxGlobal = vertexMapper_().subIndex(element, vIdx, dim);
-
- return vertexColor_[vIdxGlobal];
- }
-
- /*!
- * \brief Returns the reassemble color of a vertex
- *
- * \param vIdxGlobal The global index of the vertex.
- */
- int vertexColor(int vIdxGlobal) const
- {
- if (!enablePartialReassemble)
- return Red; // reassemble unconditionally!
- return vertexColor_[vIdxGlobal];
- }
-
- /*!
- * \brief Returns the Jacobian reassemble color of an element
- *
- * \param element The Codim-0 DUNE entity
- */
- int elementColor(const Element &element) const
- {
- if (!enablePartialReassemble)
- return Red; // reassemble unconditionally!
-
- int eIdxGlobal = elementMapper_().index(element);
- return elementColor_[eIdxGlobal];
- }
-
- /*!
- * \brief Returns the Jacobian reassemble color of an element
- *
- * \param globalElementIdx The global index of the element.
- */
- int elementColor(int globalElementIdx) const
- {
- if (!enablePartialReassemble)
- return Red; // reassemble unconditionally!
- return elementColor_[globalElementIdx];
- }
-
- /*!
- * \brief Returns a pointer to the PDELab's grid function space.
- */
- const GridFunctionSpace& gridFunctionSpace() const
- {
- return *gridFunctionSpace_;
- }
-
- /*!
- * \brief Returns a pointer to the PDELab's constraints
- * transformation.
- */
- const ConstraintsTrafo& constraintsTrafo() const
- {
- return *constraintsTrafo_;
- }
-
- /*!
- * \brief Return constant reference to global Jacobian matrix.
- */
- const JacobianMatrix& matrix() const
- { return *matrix_; }
-
- /*!
- * \brief Return reference to global Jacobian matrix.
- */
- JacobianMatrix& matrix()
- { return *matrix_; }
-
- /*!
- * \brief Return constant reference to global residual vector.
- */
- const SolutionVector& residual() const
- { return *residual_; }
-
-
- /*!
- * \brief Return reference to global residual vector.
- */
- SolutionVector& residual()
- { return *residual_; }
-
- const GridOperator &gridOperator() const
- { return *gridOperator_;}
-
-private:
- // reset the global linear system of equations. if partial
- // reassemble is enabled, this means that the jacobian matrix must
- // only be erased partially!
- void resetSystem_()
- {
- // always reset the right hand side.
- *residual_ = 0.0;
-
- if (!enablePartialReassemble) {
- // If partial reassembly of the jacobian is not enabled,
- // we can just reset everything!
- (*matrix_) = 0;
- return;
- }
-
- // reset all entries corrosponding to a red vertex
- for (int rowIdx = 0; rowIdx < matrix_->N(); ++rowIdx) {
- if (vertexColor_[rowIdx] == Green)
- continue; // the equations for this control volume are
- // already below the treshold
-
- // set all entries in the row to 0
- typedef typename JacobianMatrix::ColIterator ColIterator;
- ColIterator colIt = (*matrix_)[rowIdx].begin();
- const ColIterator &colEndIt = (*matrix_)[rowIdx].end();
- for (; colIt != colEndIt; ++colIt) {
- (*colIt) = 0.0;
- }
- };
- }
-
- Problem &problem_()
- { return *problemPtr_; }
- const Problem &problem_() const
- { return *problemPtr_; }
- const Model &model_() const
- { return problem_().model(); }
- Model &model_()
- { return problem_().model(); }
- const GridView &gridView_() const
- { return problem_().gridView(); }
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); }
- const ElementMapper &elementMapper_() const
- { return problem_().elementMapper(); }
-
- Problem *problemPtr_;
-
- // the jacobian matrix
- std::shared_ptr matrix_;
- // the right-hand side
- std::shared_ptr residual_;
-
- // attributes required for jacobian matrix recycling
- bool reuseMatrix_;
-
- // attributes required for partial jacobian reassembly
- std::vector vertexColor_;
- std::vector elementColor_;
- std::vector vertexDelta_;
-
- int totalElems_;
- int greenElems_;
-
- Scalar nextReassembleTolerance_;
- Scalar reassembleTolerance_;
-
-
- std::shared_ptr constraints_;
- std::shared_ptr pressureFEM_;
- std::shared_ptr displacementFEM_;
- std::shared_ptr pressureScalarGFS_;
- std::shared_ptr displacementScalarGFS_;
- std::shared_ptr pressureGFS_;
- std::shared_ptr displacementGFS_;
- std::shared_ptr gridFunctionSpace_;
- std::shared_ptr constraintsTrafo_;
- std::shared_ptr localOperator_;
- std::shared_ptr gridOperator_;
-};
-
-} // namespace PDELab
-
-} // namespace Dumux
-
-#endif
diff --git a/dumux/geomechanics/el2p/basemodel.hh b/dumux/geomechanics/el2p/basemodel.hh
deleted file mode 100644
index 35f1a67ff36eb1bdf4d829d4463eea4b5926a9be..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/basemodel.hh
+++ /dev/null
@@ -1,984 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Base class for fully-implicit models
- */
-#ifndef DUMUX_ELASTIC2P_BASE_MODEL_HH
-#define DUMUX_ELASTIC2P_BASE_MODEL_HH
-
-#include
-#include
-
-#include
-#include
-#include
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ElTwoPBoxModel
- * \brief base class for the two-phase geomechanics model
- */
-template
-class ElTwoPBaseModel
-{
- typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, NumEq),
- dim = GridView::dimension
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) LocalJacobian;
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
- typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
-
- typedef typename GridView::ctype CoordScalar;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- typedef typename Dune::ReferenceElements ReferenceElements;
- typedef typename Dune::ReferenceElement ReferenceElement;
-
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dim : 0 };
-
- // copying a model is not a good idea
- ElTwoPBaseModel(const ElTwoPBaseModel &);
-
-public:
- /*!
- * \brief The constructor.
- */
- ElTwoPBaseModel()
- : problemPtr_(0)
- {
- enableHints_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableHints);
- }
-
- /*!
- * \brief Apply the initial conditions to the model.
- *
- * \param problem The object representing the problem which needs to
- * be simulated.
- */
- void init(Problem &problem)
- {
- problemPtr_ = &problem;
-
- updateBoundaryIndices_();
-
- int numDofs = asImp_().numDofs();
- if (isBox)
- boxVolume_.resize(numDofs);
-
- localJacobian_.init(problem_());
- jacAsm_ = std::make_shared();
- jacAsm_->init(problem_());
-
- uCur_ = std::make_shared(jacAsm_->gridFunctionSpace());
- uPrev_ = std::make_shared(jacAsm_->gridFunctionSpace());
-
- asImp_().applyInitialSolution_();
-
- // resize the hint vectors
- if (isBox && enableHints_) {
- int numVertices = gridView_().size(dim);
- curHints_.resize(numVertices);
- prevHints_.resize(numVertices);
- hintsUsable_.resize(numVertices);
- std::fill(hintsUsable_.begin(),
- hintsUsable_.end(),
- false);
- }
-
- // also set the solution of the "previous" time step to the
- // initial solution.
- *uPrev_ = *uCur_;
- }
-
- void setHints(const Element &element,
- ElementVolumeVariables &prevVolVars,
- ElementVolumeVariables &curVolVars) const
- {
- if (!isBox || !enableHints_)
- return;
-
- int n = element.subEntities(dim);
- prevVolVars.resize(n);
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i) {
- int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
-
- if (!hintsUsable_[vIdxGlobal]) {
- curVolVars[i].setHint(NULL);
- prevVolVars[i].setHint(NULL);
- }
- else {
- curVolVars[i].setHint(&curHints_[vIdxGlobal]);
- prevVolVars[i].setHint(&prevHints_[vIdxGlobal]);
- }
- }
- }
-
- void setHints(const Element &element,
- ElementVolumeVariables &curVolVars) const
- {
- if (!isBox || !enableHints_)
- return;
-
- int n = element.subEntities(dim);
- curVolVars.resize(n);
- for (int i = 0; i < n; ++i) {
- int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
-
- if (!hintsUsable_[vIdxGlobal])
- curVolVars[i].setHint(NULL);
- else
- curVolVars[i].setHint(&curHints_[vIdxGlobal]);
- }
- }
-
- void updatePrevHints()
- {
- if (!isBox || !enableHints_)
- return;
-
- prevHints_ = curHints_;
- }
-
- void updateCurHints(const Element &element,
- const ElementVolumeVariables &elemVolVars) const
- {
- if (!isBox || !enableHints_)
- return;
-
- for (unsigned int i = 0; i < elemVolVars.size(); ++i) {
- int vIdxGlobal = vertexMapper().subIndex(element, i, dim);
- curHints_[vIdxGlobal] = elemVolVars[i];
- if (!hintsUsable_[vIdxGlobal])
- prevHints_[vIdxGlobal] = elemVolVars[i];
- hintsUsable_[vIdxGlobal] = true;
- }
- }
-
-
- /*!
- * \brief Compute the global residual for an arbitrary solution
- * vector.
- *
- * \param residual Stores the result
- * \param u The solution for which the residual ought to be calculated
- */
- Scalar globalResidual(SolutionVector &residual,
- const SolutionVector &u)
- {
- jacAsm_->gridOperator().residual(u, residual);
-
- // calculate the square norm of the residual
- Scalar result2 = residual.base().two_norm2();
- if (gridView_().comm().size() > 1)
- result2 = gridView_().comm().sum(result2);
-
- // add up the residuals on the process borders
- if (isBox && gridView_().comm().size() > 1) {
- VertexHandleSum
- sumHandle(residual, vertexMapper());
- gridView_().communicate(sumHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
- }
-
- using std::sqrt;
- return sqrt(result2);
- }
-
- /*!
- * \brief Compute the global residual for the current solution
- * vector.
- *
- * \param residual Stores the result
- */
- Scalar globalResidual(SolutionVector &residual)
- {
- return globalResidual(residual, curSol());
- }
-
- /*!
- * \brief Compute the integral over the domain of the storage
- * terms of all conservation quantities.
- *
- * \param storage Stores the result
- */
- void globalStorage(PrimaryVariables &storage)
- {
- storage = 0;
-
- for (const auto& element : elements(gridView_(), Dune::Partitions::interior))
- {
- localResidual().evalStorage(element);
-
- if (isBox)
- {
- for (int i = 0; i < element.subEntities(dim); ++i)
- storage += localResidual().storageTerm()[i];
- }
- else
- {
- storage += localResidual().storageTerm()[0];
- }
- }
-
- if (gridView_().comm().size() > 1)
- storage = gridView_().comm().sum(storage);
- }
-
- /*!
- * \brief Returns the volume \f$\mathrm{[m^3]}\f$ of a given control volume.
- *
- * \param vIdxGlobal The global index of the control volume's
- * associated vertex
- */
- Scalar boxVolume(const int vIdxGlobal) const
- {
- if (isBox)
- {
- return boxVolume_[vIdxGlobal][0];
- }
- else
- {
- DUNE_THROW(Dune::InvalidStateException,
- "requested box volume for cell-centered model");
- }
- }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- const SolutionVector &curSol() const
- { return *uCur_; }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- SolutionVector &curSol()
- { return *uCur_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- const SolutionVector &prevSol() const
- { return *uPrev_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- SolutionVector &prevSol()
- { return *uPrev_; }
-
- /*!
- * \brief Returns the operator assembler for the global jacobian of
- * the problem.
- */
- JacobianAssembler &jacobianAssembler()
- { return *jacAsm_; }
-
- /*!
- * \copydoc jacobianAssembler()
- */
- const JacobianAssembler &jacobianAssembler() const
- { return *jacAsm_; }
-
- /*!
- * \brief Returns the local jacobian which calculates the local
- * stiffness matrix for an arbitrary element.
- *
- * The local stiffness matrices of the element are used by
- * the jacobian assembler to produce a global linerization of the
- * problem.
- */
- LocalJacobian &localJacobian()
- { return localJacobian_; }
- /*!
- * \copydoc localJacobian()
- */
- const LocalJacobian &localJacobian() const
- { return localJacobian_; }
-
- /*!
- * \brief Returns the local residual function.
- */
- LocalResidual &localResidual()
- { return localJacobian().localResidual(); }
- /*!
- * \copydoc localResidual()
- */
- const LocalResidual &localResidual() const
- { return localJacobian().localResidual(); }
-
- /*!
- * \brief Returns the maximum relative shift between two vectors of
- * primary variables.
- *
- * \param priVars1 The first vector of primary variables
- * \param priVars2 The second vector of primary variables
- */
- Scalar relativeShiftAtDof(const PrimaryVariables &priVars1,
- const PrimaryVariables &priVars2)
- {
- Scalar result = 0.0;
- using std::abs;
- using std::max;
- for (int j = 0; j < numEq; ++j) {
- Scalar eqErr = abs(priVars1[j] - priVars2[j]);
- eqErr /= max(1.0, abs(priVars1[j] + priVars2[j])/2);
-
- result = max(result, eqErr);
- }
- return result;
- }
-
- /*!
- * \brief Try to progress the model to the next timestep.
- *
- * \param solver The non-linear solver
- * \param controller The controller which specifies the behaviour
- * of the non-linear solver
- */
- bool update(NewtonMethod &solver,
- NewtonController &controller)
- {
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().base().size(); ++i)
- Valgrind::CheckDefined(curSol().base()[i]);
-#endif // HAVE_VALGRIND
-
- asImp_().updateBegin();
-
- bool converged = solver.execute(controller);
- if (converged) {
- asImp_().updateSuccessful();
- }
- else
- asImp_().updateFailed();
-
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().base().size(); ++i) {
- Valgrind::CheckDefined(curSol().base()[i]);
- }
-#endif // HAVE_VALGRIND
-
- return converged;
- }
-
- /*!
- * \brief Check the plausibility of the current solution
- *
- * This has to be done by the actual model, it knows
- * best, what (ranges of) variables to check.
- * This is primarily a hook
- * which the actual model can overload.
- */
- void checkPlausibility() const
- { }
-
- /*!
- * \brief Called by the update() method before it tries to
- * apply the newton method. This is primarily a hook
- * which the actual model can overload.
- */
- void updateBegin()
- { }
-
-
- /*!
- * \brief Called by the update() method if it was
- * successful. This is primarily a hook which the actual
- * model can overload.
- */
- void updateSuccessful()
- { }
-
- /*!
- * \brief Called by the update() method if it was
- * unsuccessful. This is primarily a hook which the actual
- * model can overload.
- */
- void updateFailed()
- {
- // Reset the current solution to the one of the
- // previous time step so that we can start the next
- // update at a physically meaningful solution.
- *uCur_ = *uPrev_;
- if (isBox)
- curHints_ = prevHints_;
-
- jacAsm_->reassembleAll();
- }
-
- /*!
- * \brief Called by the problem if a time integration was
- * successful, post processing of the solution is done and
- * the result has been written to disk.
- *
- * This should prepare the model for the next time integration.
- */
- void advanceTimeLevel()
- {
- // make the current solution the previous one.
- *uPrev_ = *uCur_;
- if (isBox)
- prevHints_ = curHints_;
-
- updatePrevHints();
- }
-
- /*!
- * \brief Serializes the current state of the model.
- *
- * \tparam Restarter The type of the serializer class
- *
- * \param res The serializer object
- */
- template
- void serialize(Restarter &res)
- {
- if (isBox)
- res.template serializeEntities(asImp_(), gridView_());
- else
- res.template serializeEntities<0>(asImp_(), gridView_());
- }
-
- /*!
- * \brief Deserializes the state of the model.
- *
- * \tparam Restarter The type of the serializer class
- *
- * \param res The serializer object
- */
- template
- void deserialize(Restarter &res)
- {
- if (isBox)
- res.template deserializeEntities(asImp_(), gridView_());
- else
- res.template deserializeEntities<0>(asImp_(), gridView_());
-
- prevSol() = curSol();
- }
-
- /*!
- * \brief Write the current solution for a vertex to a restart
- * file.
- *
- * \param outstream The stream into which the vertex data should
- * be serialized to
- * \param entity The entity which's data should be
- * serialized, i.e. a vertex for the box method
- * and an element for the cell-centered method
- */
- template
- void serializeEntity(std::ostream &outstream,
- const Entity &entity)
- {
- int dofIdxGlobal = dofMapper().index(entity);
-
- // write phase state
- if (!outstream.good()) {
- DUNE_THROW(Dune::IOError,
- "Could not serialize vertex "
- << dofIdxGlobal);
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- outstream << curSol()[dofIdxGlobal][eqIdx] << " ";
- }
- }
-
- /*!
- * \brief Reads the current solution variables for a vertex from a
- * restart file.
- *
- * \param instream The stream from which the vertex data should
- * be deserialized from
- * \param entity The entity which's data should be
- * serialized, i.e. a vertex for the box method
- * and an element for the cell-centered method
- */
- template
- void deserializeEntity(std::istream &instream,
- const Entity &entity)
- {
- int dofIdxGlobal = dofMapper().index(entity);
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!instream.good())
- DUNE_THROW(Dune::IOError,
- "Could not deserialize vertex "
- << dofIdxGlobal);
- instream >> curSol()[dofIdxGlobal][eqIdx];
- }
- }
-
- /*!
- * \brief Returns the number of global degrees of freedoms (DOFs)
- */
- size_t numDofs() const
- {
- if (isBox)
- return gridView_().size(dim);
- else
- return gridView_().size(0);
- }
-
- /*!
- * \brief Mapper for the entities where degrees of freedoms are
- * defined to indices.
- *
- * Is the box method is used, this means a mapper
- * for vertices, if the cell centered method is used,
- * this means a mapper for elements.
- */
- template
- const typename std::enable_if::type &dofMapper() const
- {
- return problem_().vertexMapper();
- }
- template
- const typename std::enable_if::type &dofMapper() const
- {
- return problem_().elementMapper();
- }
-
- /*!
- * \brief Mapper for vertices to indices.
- */
- const VertexMapper &vertexMapper() const
- { return problem_().vertexMapper(); }
-
- /*!
- * \brief Mapper for elements to indices.
- */
- const ElementMapper &elementMapper() const
- { return problem_().elementMapper(); }
-
- /*!
- * \brief Resets the Jacobian matrix assembler, so that the
- * boundary types can be altered.
- */
- void resetJacobianAssembler ()
- {
- jacAsm_.template reset(0);
- jacAsm_ = std::make_shared();
- jacAsm_->init(problem_());
- }
-
- /*!
- * \brief Update the weights of all primary variables within an
- * element given the complete set of volume variables
- *
- * \param element The DUNE codim 0 entity
- * \param volVars All volume variables for the element
- */
- void updatePVWeights(const Element &element,
- const ElementVolumeVariables &volVars) const
- { }
-
- /*!
- * \brief Add the vector fields for analysing the convergence of
- * the newton method to the a VTK multi writer.
- *
- * \tparam MultiWriter The type of the VTK multi writer
- *
- * \param writer The VTK multi writer object on which the fields should be added.
- * \param u The solution function
- * \param deltaU The delta of the solution function before and after the Newton update
- */
- template
- void addConvergenceVtkFields(MultiWriter &writer,
- const SolutionVector &u,
- const SolutionVector &deltaU)
- {
- typedef Dune::BlockVector > ScalarField;
-
- SolutionVector residual(u);
- asImp_().globalResidual(residual, u);
-
- // create the required scalar fields
- unsigned numDofs = asImp_().numDofs();
-
- // global defect of the two auxiliary equations
- ScalarField* def[numEq];
- ScalarField* delta[numEq];
- ScalarField* x[numEq];
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- x[eqIdx] = writer.allocateManagedBuffer(numDofs);
- delta[eqIdx] = writer.allocateManagedBuffer(numDofs);
- def[eqIdx] = writer.allocateManagedBuffer(numDofs);
- }
-
- for (unsigned int vIdxGlobal = 0; vIdxGlobal < u.base().size(); vIdxGlobal++)
- {
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- (*x[eqIdx])[vIdxGlobal] = u.base()[vIdxGlobal][eqIdx];
- (*delta[eqIdx])[vIdxGlobal] = - deltaU.base()[vIdxGlobal][eqIdx];
- (*def[eqIdx])[vIdxGlobal] = residual.base()[vIdxGlobal][eqIdx];
- }
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- std::ostringstream oss;
- oss.str(""); oss << "x_" << eqIdx;
- if (isBox)
- writer.attachVertexData(*x[eqIdx], oss.str());
- else
- writer.attachCellData(*x[eqIdx], oss.str());
- oss.str(""); oss << "delta_" << eqIdx;
- if (isBox)
- writer.attachVertexData(*delta[eqIdx], oss.str());
- else
- writer.attachCellData(*delta[eqIdx], oss.str());
- oss.str(""); oss << "defect_" << eqIdx;
- if (isBox)
- writer.attachVertexData(*def[eqIdx], oss.str());
- else
- writer.attachCellData(*def[eqIdx], oss.str());
- }
-
- asImp_().addOutputVtkFields(u, writer);
- }
-
- /*!
- * \brief Add the quantities of a time step which ought to be written to disk.
- *
- * This should be overwritten by the actual model if any secondary
- * variables should be written out. Read: This should _always_ be
- * overwritten by well behaved models!
- *
- * \tparam MultiWriter The type of the VTK multi writer
- *
- * \param sol The global vector of primary variable values.
- * \param writer The VTK multi writer where the fields should be added.
- */
- template
- void addOutputVtkFields(const SolutionVector &sol,
- MultiWriter &writer)
- {
- typedef Dune::BlockVector > ScalarField;
-
- // create the required scalar fields
- unsigned numDofs = asImp_().numDofs();
-
- // global defect of the two auxiliary equations
- ScalarField* x[numEq];
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- x[eqIdx] = writer.allocateManagedBuffer(numDofs);
- }
-
- for (int vIdxGlobal = 0; vIdxGlobal < sol.size(); vIdxGlobal++)
- {
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- (*x[eqIdx])[vIdxGlobal] = sol[vIdxGlobal][eqIdx];
- }
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- std::ostringstream oss;
- oss << "primaryVar_" << eqIdx;
- if (isBox)
- writer.attachVertexData(*x[eqIdx], oss.str());
- else
- writer.attachCellData(*x[eqIdx], oss.str());
- }
- }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns true if the entity indicated by 'dofIdxGlobal'
- * is located on / touches the grid's boundary.
- *
- * \param dofIdxGlobal The global index of the entity
- */
- bool onBoundary(const int dofIdxGlobal) const
- { return boundaryIndices_[dofIdxGlobal]; }
-
- /*!
- * \brief Returns true if a vertex is located on the grid's
- * boundary.
- *
- * \param element A DUNE Codim<0> entity which contains the control
- * volume's associated vertex.
- * \param vIdx The local vertex index inside element
- */
- bool onBoundary(const Element &element, const int vIdx) const
- {
- if (isBox)
- return onBoundary(vertexMapper().subIndex(element, vIdx, dim));
- else
- DUNE_THROW(Dune::InvalidStateException,
- "requested for cell-centered model");
- }
-
-
- /*!
- * \brief Returns true if the control volume touches
- * the grid's boundary.
- *
- * \param element A DUNE Codim<0> entity coinciding with the control
- * volume.
- */
- bool onBoundary(const Element &element) const
- {
- if (!isBox)
- return onBoundary(elementMapper().index(element));
- else
- DUNE_THROW(Dune::InvalidStateException,
- "requested for box model");
- }
-
- /*!
- * \brief Fill the fluid state according to the primary variables.
- *
- * Taking the information from the primary variables,
- * the fluid state is filled with every information that is
- * necessary to evaluate the model's local residual.
- *
- * \param priVars The primary variables of the model.
- * \param problem The problem at hand.
- * \param element The current element.
- * \param fvGeometry The finite volume element geometry.
- * \param scvIdx The index of the subcontrol volume.
- * \param fluidState The fluid state to fill.
- */
- template
- static void completeFluidState(const PrimaryVariables& priVars,
- const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const int scvIdx,
- FluidState& fluidState)
- {
- VolumeVariables::completeFluidState(priVars, problem, element,
- fvGeometry, scvIdx, fluidState);
- }
-protected:
- /*!
- * \brief A reference to the problem on which the model is applied.
- */
- Problem &problem_()
- { return *problemPtr_; }
- /*!
- * \copydoc problem_()
- */
- const Problem &problem_() const
- { return *problemPtr_; }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Reference to the local residal object
- */
- LocalResidual &localResidual_()
- { return localJacobian_.localResidual(); }
-
- /*!
- * \brief Applies the initial solution for all vertices of the grid.
- */
- void applyInitialSolution_()
- {
- // first set the whole domain to zero
- *uCur_ = Scalar(0.0);
- boxVolume_ = Scalar(0.0);
-
- FVElementGeometry fvGeometry;
-
- // iterate through leaf grid and evaluate initial
- // condition at the center of each sub control volume
- //
- // the initial condition needs to be unique for
- // each vertex. we should think about the API...
- for (const auto& element : elements(gridView_())) {
- // deal with the current element
- fvGeometry.update(gridView_(), element);
-
- // loop over all element vertices, i.e. sub control volumes
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; scvIdx++)
- {
- // get the global index of the degree of freedom
- int dofIdxGlobal = dofMapper().subIndex(element, scvIdx, dofCodim);
-
- // let the problem do the dirty work of nailing down
- // the initial solution.
- PrimaryVariables initPriVars;
- Valgrind::SetUndefined(initPriVars);
- problem_().initial(initPriVars,
- element,
- fvGeometry,
- scvIdx);
- Valgrind::CheckDefined(initPriVars);
-
- if (isBox)
- {
- // add up the initial values of all sub-control
- // volumes. If the initial values disagree for
- // different sub control volumes, the initial value
- // will be the arithmetic mean.
- initPriVars *= fvGeometry.subContVol[scvIdx].volume;
- boxVolume_[dofIdxGlobal] += fvGeometry.subContVol[scvIdx].volume;
- }
-
- uCur_->base()[dofIdxGlobal] += initPriVars;
- Valgrind::CheckDefined(uCur_->base()[dofIdxGlobal]);
- }
- }
-
- // add up the primary variables and the volumes of the boxes
- // which cross process borders
- if (isBox && gridView_().comm().size() > 1) {
- VertexHandleSum,
- Dune::BlockVector >,
- VertexMapper> sumVolumeHandle(boxVolume_, vertexMapper());
- gridView_().communicate(sumVolumeHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
-
- VertexHandleSum
- sumPVHandle(uCur_->base(), vertexMapper());
- gridView_().communicate(sumPVHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
- }
-
- if (isBox)
- {
- // divide all primary variables by the volume of their boxes
- for (unsigned int i = 0; i < uCur_->base().size(); ++i) {
- uCur_->base()[i] /= boxVolume(i);
- }
- }
- }
-
- /*!
- * \brief Find all indices of boundary vertices (box) / elements (cell centered).
- */
- void updateBoundaryIndices_()
- {
- boundaryIndices_.resize(numDofs());
- std::fill(boundaryIndices_.begin(), boundaryIndices_.end(), false);
-
- for (const auto& element : elements(gridView_())) {
- Dune::GeometryType geomType = element.geometry().type();
- const ReferenceElement &refElement = ReferenceElements::general(geomType);
-
- for (const auto& intersection : intersections(gridView_(), element)) {
- if (intersection.boundary()) {
- if (isBox)
- {
- // add all vertices on the intersection to the set of
- // boundary vertices
- int fIdx = intersection.indexInInside();
- int numFaceVerts = refElement.size(fIdx, 1, dim);
- for (int faceVertexIdx = 0;
- faceVertexIdx < numFaceVerts;
- ++faceVertexIdx)
- {
- int vIdx = refElement.subEntity(fIdx,
- 1,
- faceVertexIdx,
- dim);
- int vIdxGlobal = vertexMapper().subIndex(element, vIdx, dim);
- boundaryIndices_[vIdxGlobal] = true;
- }
- }
- else
- {
- int eIdxGlobal = elementMapper().index(element);
- boundaryIndices_[eIdxGlobal] = true;
- }
- }
- }
- }
- }
-
- // the hint cache for the previous and the current volume
- // variables
- mutable std::vector hintsUsable_;
- mutable std::vector curHints_;
- mutable std::vector prevHints_;
-
- // the problem we want to solve. defines the constitutive
- // relations, matxerial laws, etc.
- Problem *problemPtr_;
-
- // calculates the local jacobian matrix for a given element
- LocalJacobian localJacobian_;
- // Linearizes the problem at the current time step using the
- // local jacobian
- std::shared_ptr jacAsm_;
-
- // the set of all indices of vertices on the boundary
- std::vector boundaryIndices_;
-
- // cur is the current iterative solution, prev the converged
- // solution of the previous time step
- std::shared_ptr uCur_;
- std::shared_ptr uPrev_;
-
- Dune::BlockVector > boxVolume_;
-
-private:
- /*!
- * \brief Returns whether messages should be printed
- */
- bool verbose_() const
- { return gridView_().comm().rank() == 0; }
-
- Implementation &asImp_()
- { return *static_cast(this); }
- const Implementation &asImp_() const
- { return *static_cast(this); }
-
- bool enableHints_;
-};
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/geomechanics/el2p/elementvolumevariables.hh b/dumux/geomechanics/el2p/elementvolumevariables.hh
deleted file mode 100644
index 0b2f3690023a15b469358f594e55dbc32aeae797..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/elementvolumevariables.hh
+++ /dev/null
@@ -1,303 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Volume variables gathered on an element
- */
-#ifndef DUMUX_BOX_EL2P_ELEMENT_VOLUME_VARIABLES_HH
-#define DUMUX_BOX_EL2P_ELEMENT_VOLUME_VARIABLES_HH
-
-#include
-
-#include
-#include "properties.hh"
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ElTwoPBoxModel
- *
- * \brief This class stores an array of VolumeVariables objects, one
- * volume variables object for each of the element's vertices
- */
-template
-class ElTwoPElementVolumeVariables : public std::vector
-{
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename Element::Geometry::JacobianInverseTransposed JacobianInverseTransposed;
-
- enum { dim = GridView::dimension };
- enum { dimWorld = GridView::dimensionworld };
- typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
-
- typedef Dune::PDELab::LocalFunctionSpace LocalFunctionSpace;
-
- typedef Dune::FieldVector GlobalPosition;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- enum {
- pressureIdx = Indices::pressureIdx,
- saturationIdx = Indices::saturationIdx
- };
-
-public:
- /*!
- * \brief The constructor.
- */
- ElTwoPElementVolumeVariables()
- { }
-
- /*!
- * \brief Construct the volume variables for all vertices of an element.
- *
- * \param problem The problem which needs to be simulated.
- * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
- * \param fvGeometry The finite volume geometry of the element
- * \param isOldSol Tells whether the model's previous or current solution should be used.
- *
- * This class is required for the update of the effective porosity values at the
- * vertices since it is a function of the divergence of the solid displacement
- * at the integration points
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- bool isOldSol)
- {
- // retrieve the current or the previous solution vector and write the values into globalSol
- const SolutionVector &globalSol =
- isOldSol?
- problem.model().prevSol():
- problem.model().curSol();
-
- const GridFunctionSpace& gridFunctionSpace = problem.model().jacobianAssembler().gridFunctionSpace();
- const typename GridFunctionSpace::Ordering& ordering = gridFunctionSpace.ordering();
- // copy the values of the globalSol vector to the localFunctionSpace values of the current element
- LocalFunctionSpace localFunctionSpace(gridFunctionSpace);
- localFunctionSpace.bind(element);
- std::vector values(localFunctionSpace.size());
- for (typename LocalFunctionSpace::Traits::IndexContainer::size_type k=0; k::Type PressSatLFS;
- const PressSatLFS& pressSatLFS = localFunctionSpace.template child<0>();
- // local function space for pressure
- typedef typename PressSatLFS::template Child<0>::Type PressLFS;
- const PressLFS& pressLFS = pressSatLFS.template child<0>();
- // local function space for saturation
- typedef typename PressSatLFS::template Child<1>::Type SatLFS;
- const SatLFS& satLFS = pressSatLFS.template child<1>();
- // local function space for solid displacement
- typedef typename LocalFunctionSpace::template Child<1>::Type DisplacementLFS;
- const DisplacementLFS& displacementLFS = localFunctionSpace.template child<1>();
- typedef typename DisplacementLFS::template Child<0>::Type ScalarDispLFS;
-
- int numScv = element.subEntities(dim);
- this->resize(numScv);
-
- for (int scvIdx = 0; scvIdx < numScv; scvIdx++)
- {
- // solution vector solI for each vertex
- PrimaryVariables solI;
- // pressure and saturation values
- solI[pressureIdx] = values[pressLFS.localIndex(scvIdx)];
- solI[saturationIdx] = values[satLFS.localIndex(scvIdx)];
- // solid displacement values for each coordinate direction
- for (int coordDir = 0; coordDir < dim; coordDir++)
- {
- const ScalarDispLFS& scalarDispLFS = displacementLFS.child(coordDir);
- solI[Indices::u(coordDir)] = values[scalarDispLFS.localIndex(scvIdx)];
- }
- // reset evaluation point to zero
- (*this)[scvIdx].setEvalPoint(0);
-
- (*this)[scvIdx].update(solI,
- problem,
- element,
- fvGeometry,
- scvIdx,
- isOldSol);
-
- Valgrind::CheckDefined((*this)[scvIdx]);
- }
- this->updateEffPorosity(problem, element, fvGeometry, isOldSol);
-
- if (isOldSol)
- prevValues_ = values;
- else
- dofValues_ = values;
- }
-
- /*!
- * \brief Update the effective porosities for all vertices of an element.
- *
- * \param problem The problem which needs to be simulated.
- * \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
- * \param fvGeometry The finite volume geometry of the element
- * \param isOldSol Specifies whether this is the previous solution or the current one
- *
- * This function is required for the update of the effective porosity values at the
- * vertices.
- *
- * During the partial derivative calculation, changes of the solid displacement
- * at vertex i can affect effective porosities of all element vertices.
- * To correctly update the effective porosities of all element vertices
- * an iteration over all scv faces is required.
- * The remaining volvars are only updated for the vertex whose primary variable
- * is changed for the derivative calculation.
- */
- void updateEffPorosity(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- bool isOldSol)
- {
- int numScv = element.subEntities(dim);
-
- // retrieve the current or the previous solution vector and write the values into globalSol
- const SolutionVector &globalSol =
- isOldSol?
- problem.model().prevSol():
- problem.model().curSol();
-
- // copy the values of the globalSol vector to the localFunctionSpace values of the current element
- const GridFunctionSpace& gridFunctionSpace = problem.model().jacobianAssembler().gridFunctionSpace();
- const typename GridFunctionSpace::Ordering& ordering = gridFunctionSpace.ordering();
- LocalFunctionSpace localFunctionSpace(gridFunctionSpace);
- localFunctionSpace.bind(element);
- std::vector values(localFunctionSpace.size());
- for (typename LocalFunctionSpace::Traits::IndexContainer::size_type k=0; k::Type DisplacementLFS;
- const DisplacementLFS& displacementLFS = localFunctionSpace.template child<1>();
- const unsigned int dispSize = displacementLFS.child(0).size();
- typedef typename DisplacementLFS::template Child<0>::Type ScalarDispLFS;
- // further types required for gradient calculations
- typedef typename ScalarDispLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
- typedef typename ScalarDispLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::RangeFieldType RF;
- typedef Dune::FieldMatrix Tensor;
-
- for (int scvIdx = 0; scvIdx < numScv; scvIdx++)
- (*this)[scvIdx].effPorosity = 0.0;
-
- for (int scvIdx = 0; scvIdx < numScv; scvIdx++)
- {
- GlobalPosition scvCenter = fvGeometry.subContVol[scvIdx].localCenter;
-
- // evaluate gradient of displacement shape functions at the center of
- // the sub control volume in the reference element
- std::vector vRefShapeGradient(dispSize);
- displacementLFS.child(0).finiteElement().localBasis().evaluateJacobian(scvCenter, vRefShapeGradient);
-
- // transform gradient to element in global coordinates
- const JacobianInverseTransposed jacInvT = element.geometry().jacobianInverseTransposed(scvCenter);
- std::vector > vShapeGradient(dispSize);
-
- // loop over element vertices
- for (size_t i = 0; i < dispSize; i++)
- {
- vShapeGradient[i] = 0.0;
- jacInvT.umv(vRefShapeGradient[i][0],vShapeGradient[i]);
- }
-
- // calculate gradient of current displacement
- // (gradient of a vector is a tensor)
- Tensor uGradient(0.0);
- // loop over coordinate directions
- for(int coordDir = 0; coordDir < dim; ++coordDir)
- {
- const ScalarDispLFS & scalarDispLFS = displacementLFS.child(coordDir);
- // loop over element vertices
- for (size_t i = 0; i < scalarDispLFS.size(); i++)
- uGradient[coordDir].axpy((*this)[i].displacement(coordDir), vShapeGradient[i]);
- }
-
- // calculate the divergence of u
- (*this)[scvIdx].divU = 0.0;
-
- for (int coordDir = 0; coordDir < dim; coordDir++)
- (*this)[scvIdx].divU += uGradient[coordDir][coordDir];
-
- // calculate the effective porosity
- if(problem.coupled() == true)
- {
- if ((*this)[scvIdx].divU < -(*this)[scvIdx].porosity())
- {
- (*this)[scvIdx].effPorosity = (*this)[scvIdx].porosity();
- std::cout<<"volume change too large"<& dofValues() const
- {
- return dofValues_;
- }
-
- Scalar& dofValues(int k)
- {
- return dofValues_[k];
- }
-
- const std::vector& prevValues() const
- {
- return prevValues_;
- }
-
-private:
- std::vector dofValues_;
- std::vector prevValues_;
-};
-
-} // namespace Dumux
-
-#endif
diff --git a/dumux/geomechanics/el2p/fluxvariables.hh b/dumux/geomechanics/el2p/fluxvariables.hh
deleted file mode 100644
index 157e010c6f244b15a345599d8799fbe57c227b80..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/fluxvariables.hh
+++ /dev/null
@@ -1,245 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file contains the calculation of all the fluxes over the surface of the
- * finite volume that make up the volume, the mass and the momentum balance
- * for the two-phase linear-elastic model.
- *
- * This means pressure, concentration and solid-displacement gradients, phase densities at
- * the integration point, etc.
- *
- * This class inherits from the two-phase model FluxVariables
- */
-#ifndef DUMUX_EL2P_FLUX_VARIABLES_HH
-#define DUMUX_EL2P_FLUX_VARIABLES_HH
-
-#include
-#include
-#include "properties.hh"
-
-namespace Dumux
-{
-
-namespace Properties
-{
-// forward declaration of properties
-NEW_PROP_TAG(SpatialParams);
-}
-/*!
- * \ingroup ElTwoPBoxModel
- * \ingroup ImplicitFluxVariables
- * \brief This template class contains the data which is required to
- * calculate the fluxes over the surface of the
- * finite volume that make up the volume, the mass and the momentum balance
- * for the two-phase linear-elastic model.
- *
- * This means pressure, concentration and solid-displacement gradients, phase densities at
- * the integration point, etc.
- *
- */
-template
-class ElTwoPFluxVariables: public ImplicitDarcyFluxVariables
-{
- friend class ImplicitDarcyFluxVariables; // be friends with parent
- typedef ImplicitDarcyFluxVariables ParentType;
-
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- enum
- {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
- typedef typename GridView::ctype CoordScalar;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldVector GlobalPosition;
- typedef Dune::FieldVector DimVector;
- typedef Dune::FieldMatrix DimMatrix;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
- enum {numEq = GET_PROP_VALUE(TypeTag, NumEq)};
-
-public:
- /*!
- * \brief Compute / update the flux variables
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param fIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- {
- ParentType::update(problem, element, fvGeometry, fIdx, elemVolVars);
-
- dU_ = 0.0;
- timeDerivUNormal_ = 0.0;
-
- elTwoPGradients_(problem, element, elemVolVars);
- calculateDDt_(problem, element, elemVolVars);
- }
-
-public:
- /*!
- * \brief Return change of u [m] with time at integration point
- * point.
- */
- Scalar dU(int dimIdx) const
- {
- return dU_[dimIdx];
- }
-
- /*!
- * \brief Return time derivative of u [m/s] in normal
- * direction at integration point
- */
- Scalar timeDerivUNormal() const
- {
- return timeDerivUNormal_;
- }
-
- /*
- * \brief Return the intrinsic permeability.
- */
- const DimMatrix &intrinsicPermeability() const
- {
- return K_;
- }
-
- /*
- * \brief Return the gradient of the potential for each phase.
- */
- DimVector potentialGrad(int phaseIdx) const
- {
- return this->potentialGrad_[phaseIdx];
- }
-
- const SCVFace &face() const
- {
- return this->fvGeometry_().subContVolFace[this->faceIdx_];
- }
-
-protected:
- /*!
- * \brief Calculation of the solid displacement gradients.
- *
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemVolVars The parameters stored in the considered element
- */
- void elTwoPGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
- typedef Dune::PDELab::LocalFunctionSpace LocalFunctionSpace;
- const GridFunctionSpace& gridFunctionSpace = problem.model().jacobianAssembler().gridFunctionSpace();
- const typename GridFunctionSpace::Ordering& ordering = gridFunctionSpace.ordering();
- LocalFunctionSpace localFunctionSpace(gridFunctionSpace);
- localFunctionSpace.bind(element);
- // copy values of previous solution into prevSolutionValues Vector
- std::vector prevSolutionValues(localFunctionSpace.size());
- // copy values of current solution into curSolutionValues Vector
- std::vector curSolutionValues(localFunctionSpace.size());
- for (typename LocalFunctionSpace::Traits::IndexContainer::size_type k=0; k::Type DisplacementLFS;
- const DisplacementLFS& displacementLFS = localFunctionSpace.template child<1>();
- // number of degrees of freedom for each displacement value (here number of element vertices)
- const unsigned int dispSize = displacementLFS.child(0).size();
- // type of function space of solid displacement value (one for each coordinate direction)
- typedef typename DisplacementLFS::template Child<0>::Type ScalarDispLFS;
- typedef typename ScalarDispLFS::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType RT_V;
-
- for(int coordDir = 0; coordDir < dim; ++coordDir) {
- // get displacement function space for coordinate direction coordDir
- const ScalarDispLFS & scalarDispLFS = displacementLFS.child(coordDir);
- std::vector vShape(dispSize);
- // evaluate shape functions of all element vertices for current integration point and write it into vector vShape
- scalarDispLFS.finiteElement().localBasis().evaluateFunction(face().ipLocal, vShape);
-
- dU_[coordDir] = 0;
- // subtract previous displacement value from current displacement value for each coordinate direction
- // coordDir and for each node i and interpolate values at integration point via the shape function vShape.
- // TODO: Check if evaluation of prevVolVars is possible
- for (size_t i = 0; i < dispSize; i++){
- dU_[coordDir] += (elemVolVars[i].primaryVars()[(numEq - dim)+coordDir]
- - prevSolutionValues[scalarDispLFS.localIndex(i)])*vShape[i];
- }
- }
- }
-
- /*!
- * \brief Calculation of the time derivative of solid displacement
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemVolVars The parameters stored in the considered element
- */
- void calculateDDt_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- Scalar dt = problem.timeManager().timeStepSize();
-
- DimVector tmp(0.0);
- // calculate time derivative of solid displacement vector
- for (int coordDir = 0; coordDir < dim; ++coordDir)
- tmp[coordDir] = dU(coordDir) / dt;
-
- // multiply time derivative of solid displacement vector with
- // normal vector of current scv-face
- timeDerivUNormal_ = tmp * face().normal;
- }
-
- //! time derivative of solid displacement times normal vector at integration point
- Scalar timeDerivUNormal_;
- //! change of solid displacement with time at integration point
- GlobalPosition dU_;
- // intrinsic permeability
- DimMatrix K_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/geomechanics/el2p/indices.hh b/dumux/geomechanics/el2p/indices.hh
deleted file mode 100644
index c0675daf9bf5070c13c576310187c8fe54a57984..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/indices.hh
+++ /dev/null
@@ -1,58 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Defines the primary variable and equation indices used by
- * the two-phase linear elasticity model.
- */
-#ifndef DUMUX_ELASTIC2P_INDICES_HH
-#define DUMUX_ELASTIC2P_INDICES_HH
-
-#include
-#include
-
-namespace Dumux
-{
-// \{
-
-namespace Properties
-{
-
-/*!
- * \ingroup ElTwoPBoxModel
- * \ingroup ImplicitIndices
- * \brief The indices for the two-phase linear elasticity model.
- *
- * This class inherits from the TwoPIndices and from the ElasticIndices
- */
-
-// PVOffset is set to 0 for the TwoPIndices and to 2 for the ElasticIndices since
-// the first two primary variables are the primary variables of the two-phase
-// model followed by the primary variables of the elastic model
-template
-class ElTwoPIndices : public ElasticIndices, public TwoPIndices
-{};
-
-}
-}
-
-#endif
diff --git a/dumux/geomechanics/el2p/localjacobian.hh b/dumux/geomechanics/el2p/localjacobian.hh
deleted file mode 100644
index 67a960a31e5f04369f833bcd35d3956b6ba03397..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/localjacobian.hh
+++ /dev/null
@@ -1,257 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Calculates the partial derivatives of the local residual for the Jacobian of the
- * two-phase linear elasticity model.
- */
-#ifndef DUMUX_EL2P_LOCAL_JACOBIAN_HH
-#define DUMUX_EL2P_LOCAL_JACOBIAN_HH
-
-#include
-#include
-#include "properties.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup ElTwoPBoxModel
- * \brief Calculates the partial derivatives of the local residual for the Jacobian
- *
- * Except for the evalPartialDerivatives function all functions are taken from the
- * base class ImplicitLocalJacobian
- */
-template
-class ElTwoPLocalJacobian : public ImplicitLocalJacobian
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- enum {
- dim = GridView::dimension,
- };
- typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- enum {
- pressureIdx = Indices::pressureIdx,
- saturationIdx = Indices::saturationIdx
- };
- // copying a local jacobian is not a good idea
- ElTwoPLocalJacobian(const ElTwoPLocalJacobian &);
-
-public:
- ElTwoPLocalJacobian() {}
-
- /*!
- * \brief Compute the partial derivatives to a primary variable at
- * an degree of freedom.
- *
- * This method is overwritten here since this model requires a call of the model specific
- * elementvolumevariables which updates the effective porosities correctly.
- *
- * The default implementation of this method uses numeric
- * differentiation, i.e. forward or backward differences (2nd
- * order), or central differences (3rd order). The method used is
- * determined by the "NumericDifferenceMethod" property:
- *
- * - if the value of this property is smaller than 0, backward
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
- * \f]
- *
- * - if the value of this property is 0, central
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
- * \f]
- *
- * - if the value of this property is larger than 0, forward
- * differences are used, i.e.:
- * \f[
- \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
- * \f]
- *
- * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
- * is the value of a sub-control volume's primary variable at the
- * evaluation point and \f$\epsilon\f$ is a small value larger than 0.
- *
- * \param partialDeriv The vector storing the partial derivatives of all
- * equations
- * \param storageDeriv the mass matrix contributions
- * \param col The block column index of the degree of freedom
- * for which the partial derivative is calculated.
- * Box: a sub-control volume index.
- * Cell centered: a neighbor index.
- * \param pvIdx The index of the primary variable
- * for which the partial derivative is calculated
- */
- void evalPartialDerivative_(ElementSolutionVector &partialDeriv,
- PrimaryVariables &storageDeriv,
- int col,
- int pvIdx)
- {
- typedef typename GET_PROP_TYPE(TypeTag, GridFunctionSpace) GridFunctionSpace;
- typedef Dune::PDELab::LocalFunctionSpace LocalFunctionSpace;
-
- // copy the values of the globalSol vector to the localFunctionSpace values of the current element
- const GridFunctionSpace& gridFunctionSpace = this->problemPtr_->model().jacobianAssembler().gridFunctionSpace();
- const typename GridFunctionSpace::Ordering& ordering = gridFunctionSpace.ordering();
- LocalFunctionSpace localFunctionSpace(gridFunctionSpace);
- localFunctionSpace.bind(this->element_());
- std::vector elementValues(localFunctionSpace.size());
- for (typename LocalFunctionSpace::Traits::IndexContainer::size_type k=0; kproblemPtr_->model().curSol()[ci];
- }
- // pressure and saturation local function space (mass balance equations)
- typedef typename LocalFunctionSpace::template Child<0>::Type PressSatLFS;
- const PressSatLFS& pressSatLFS = localFunctionSpace.template child<0>();
- // local function space for pressure
- typedef typename PressSatLFS::template Child<0>::Type PressLFS;
- const PressLFS& pressLFS = pressSatLFS.template child<0>();
- // local function space for saturation
- typedef typename PressSatLFS::template Child<1>::Type SatLFS;
- const SatLFS& satLFS = pressSatLFS.template child<1>();
- // local function space for solid displacement
- typedef typename LocalFunctionSpace::template Child<1>::Type DisplacementLFS;
- const DisplacementLFS& displacementLFS = localFunctionSpace.template child<1>();
- typedef typename DisplacementLFS::template Child<0>::Type ScalarDispLFS;
-
- //primary variable vector priVars for each vertex
- PrimaryVariables priVars;
- priVars[pressureIdx] = elementValues[pressLFS.localIndex(col)];
- priVars[saturationIdx] = elementValues[satLFS.localIndex(col)];
- for (int coordDir = 0; coordDir < dim; coordDir++)
- {
- const ScalarDispLFS& scalarDispLFS = displacementLFS.child(coordDir);
- priVars[Indices::u(coordDir)] = elementValues[scalarDispLFS.localIndex(col)];
- }
-
- VolumeVariables origVolVars(this->curVolVars_[col]);
- this->curVolVars_[col].setEvalPoint(&origVolVars);
- Scalar eps = this->numericEpsilon(col, pvIdx);
- Scalar delta = 0;
-
- if (this->numericDifferenceMethod_ >= 0) {
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- priVars[pvIdx] += eps;
- delta += eps;
-
- // calculate the residual
- this->curVolVars_[col].update(priVars,
- this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- col,
- false);
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- false);
-
- this->localResidual().eval(this->element_(),
- this->fvElemGeom_,
- this->prevVolVars_,
- this->curVolVars_,
- this->bcTypes_);
-
- // store the residual
- partialDeriv = this->localResidual().residual();
- storageDeriv = this->localResidual().storageTerm()[col];
- }
- else {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv = this->residual_;
- storageDeriv = this->storageTerm_[col];
- }
-
-
- if (this->numericDifferenceMethod_ <= 0) {
- // we are not using forward differences, i.e. we don't
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- priVars[pvIdx] -= delta + eps;
- delta += eps;
-
- // calculate residual again
- this->curVolVars_[col].update(priVars,
- this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- col,
- false);
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- false);
- this->localResidual().eval(this->element_(),
- this->fvElemGeom_,
- this->prevVolVars_,
- this->curVolVars_,
- this->bcTypes_);
- partialDeriv -= this->localResidual().residual();
- storageDeriv -= this->localResidual().storageTerm()[col];
-
- }
- else {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv -= this->residual_;
- storageDeriv -= this->storageTerm_[col];
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
-// if(partialDeriv[col][pvIdx] == -350045)
-
- partialDeriv /= delta;
- storageDeriv /= delta;
- // restore the orignal state of the element's volume variables
- this->curVolVars_[col] = origVolVars;
- // update the effective porosities
- this->curVolVars_.updateEffPorosity(this->problem_(),
- this->element_(),
- this->fvElemGeom_,
- false);
-
-#if HAVE_VALGRIND
- for (unsigned i = 0; i < partialDeriv.size(); ++i)
- Valgrind::CheckDefined(partialDeriv[i]);
-#endif
- }
-};
-}
-
-#endif
diff --git a/dumux/geomechanics/el2p/localoperator.hh b/dumux/geomechanics/el2p/localoperator.hh
deleted file mode 100644
index 3720aff1ec3f4ac880657575108b0b979556b93c..0000000000000000000000000000000000000000
--- a/dumux/geomechanics/el2p/localoperator.hh
+++ /dev/null
@@ -1,636 +0,0 @@
-// -*- 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 . *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file contains a local operator for PDELab which
- * wraps the contributions from
- * el2plocalresidual (box discretized mass balances)
- * and alphaMomentum (FE discretized momentum balance).
- */
-#ifndef DUMUX_EL2P_LOCAL_OPERATOR_HH
-#define DUMUX_EL2P_LOCAL_OPERATOR_HH
-
-#include
-#include
-
-#include
-#include
-#include
-#include
-#include
-#include "properties.hh"
-
-namespace Dumux {
-
-namespace PDELab {
-
-/*!
- * \brief A local operator for PDELab which wraps the contributions from
- * el2plocalresidual (box discretized mass balances)
- * and alphaMomentum (FE discretized momentum balance).
- */
-template
-class El2PLocalOperator
- :
- public Dune::PDELab::FullVolumePattern,
- public Dune::PDELab::LocalOperatorDefaultFlags
-{
- // copying the local operator for PDELab is not a good idea
- El2PLocalOperator(const El2PLocalOperator &);
-
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity::Geometry::JacobianInverseTransposed JacobianInverseTransposed;
- typedef typename GridView::Intersection Intersection;
- typedef typename Dune::PDELab::IntersectionGeometry::ctype DT;
-
- enum{numEq = GET_PROP_VALUE(TypeTag, NumEq)};
- enum{dim = GridView::dimension};
- enum{dimWorld = GridView::dimensionworld};
- typedef Dune::FieldVector GlobalPosition;
- typedef Dune::FieldVector DimVector;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-
- enum {
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx
- };
-
-public:
- // pattern assembly flags
- enum { doPatternVolume = true };
-
- // residual assembly flags
- enum { doAlphaVolume = true };
-
- /*!
- * \param model The physical model for the box scheme.
- */
- El2PLocalOperator(Model &model)
- : model_(model)
- {}
-
- /*!
- * \brief Volume integral depending on test and ansatz functions
- *
- * \tparam EG The entity geometry type from PDELab
- * \tparam LFSU The type of the local function space of the ansatz functions
- * \tparam X The type of the container for the coefficients for the ansatz functions
- * \tparam LFSV The type of the local function space of the test functions
- * \tparam R The range type (usually FieldVector)
- *
- * \param eg The entity geometry object
- * \param lfsu The local function space object of the ansatz functions
- * \param x The object of the container for the coefficients for the ansatz functions
- * \param lfsv The local function space object of the test functions
- * \param r The object storing the volume integral
- */
- template
- void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
- const LFSV& lfsv, R& r) const
- {
- typedef typename LFSU::Traits::SizeType size_type;
-
- // evaluate the local residual of the box mass balance equation for the current element
- model_.localResidual().eval(eg.entity());
-
- // pressure and saturation local function space (mass balance equations)
- typedef typename LFSU::template Child<0>::Type PressSatLFS;
- // local function space for pressure
- typedef typename PressSatLFS::template Child<0>::Type PressLFS;
- const PressSatLFS& pressSatLFS = lfsu.template child<0>();
- const PressLFS& pressLFS = pressSatLFS.template child<0>();
- // local function space for saturation
- typedef typename PressSatLFS::template Child<1>::Type SatLFS;
- const SatLFS& satLFS = pressSatLFS.template child<1>();
-
- unsigned int numScv = eg.entity().subEntities(dim);
-
- for (size_type i = 0; i < (numEq-dim) * numScv; i++)
- {
- // retrieve the local residual value for vertex=i%Vertices and equation i/numScv (here 0 or 1)
- Scalar tmp = model_.localResidual().residual(i%numScv)[i/numScv];
- // get residual for brine phase mass balance equation
- if(i < numScv)
- r.rawAccumulate(pressLFS, i, tmp);
- // get residual for CO2 phase mass balance equation
- else
- r.rawAccumulate(satLFS,i-numScv, tmp);
- }
- // get residual for momentum balance equation
- alphaMomentum(eg, lfsu, x, lfsv, r);
- }
-
-
- /*!
- * \brief Calculate the local residual of the momentum balance equation
- * with the finite element method. This requires numerical
- * integration which is done via a quadrature rule.
- *
- * \tparam EG The entity geometry type from PDELab
- * \tparam LFSU The type of the local function space of the ansatz functions
- * \tparam X The type of the container for the coefficients for the ansatz functions
- * \tparam LFSV The type of the local function space of the test functions
- * \tparam R The range type (usually FieldVector)
- *
- * \param eg The entity geometry object
- * \param lfsu The local function space object of the ansatz functions
- * \param x The object of the container for the coefficients for the ansatz functions
- * \param lfsv The local function space object of the test functions
- * \param r The object storing the volume integral
- *
- *
- */
- template
- void alphaMomentum (const EG& eg, const LFSU& lfsu, const X& x,
- const LFSV& lfsv, R& r) const
- {
- FVElementGeometry fvGeometry;
- fvGeometry.update(model_.problem().gridView(), eg.entity());
- // retrieve lame parameters for calculation of effective stresses
- const Dune::FieldVector lameParams = model_.problem().spatialParams().lameParams(eg.entity(), fvGeometry, 0);
- Scalar lambda = lameParams[0];
- Scalar mu = lameParams[1];
- // retrieve materialParams for calculate of capillary pressure
- const MaterialLawParams& materialParams =
- model_.problem().spatialParams().materialLawParams(eg.entity(), fvGeometry, 0);
- // retrieve initial porosity
- Scalar porosity = model_.problem().spatialParams().porosity(eg.entity(), fvGeometry, 0);
-
- // order of quadrature rule
- const int qorder = 3;
-
- // extract local function spaces
- // pressure and saturation local function space (mass balance equations)
- typedef typename LFSU::template Child<0>::Type PressSatLFS;
- const PressSatLFS& pressSatLFS = lfsu.template child<0>();
- // local function space for pressure
- typedef typename PressSatLFS::template Child<0>::Type PressLFS;
- const PressLFS& pressLFS = pressSatLFS.template child<0>();
- const unsigned int pressSize = pressLFS.size();
- // local function space for saturation
- typedef typename PressSatLFS::template Child<1>::Type SatLFS;
- const SatLFS& satLFS = pressSatLFS.template child<1>();
- // local function space for solid displacement
- typedef typename LFSU::template Child<1>::Type DisplacementLFS;
- typedef typename DisplacementLFS::template Child<0>::Type DisplacementScalarLFS;
- const DisplacementLFS& displacementLFS = lfsu.template child<1>();
- const DisplacementScalarLFS& uScalarLFS = displacementLFS.template child<0>();
- const unsigned int dispSize = displacementLFS.template child<0>().size();
-
- // domain and range field type
- typedef typename DisplacementScalarLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::RangeFieldType RF;
- typedef typename DisplacementScalarLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::RangeType RT_V;
- typedef typename DisplacementScalarLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
- typedef typename PressLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::DomainFieldType DF;
- typedef typename PressLFS::Traits::FiniteElementType::
- Traits::LocalBasisType::Traits::RangeType RT_P;
-
- // select quadrature rule for the element geometry type and with the order=qorder
- const auto geometry = eg.geometry();
- Dune::GeometryType geomType = geometry.type();
- const Dune::QuadratureRule& rule = Dune::QuadratureRules::rule(geomType,qorder);
-
- // loop over quadrature points
- for (typename Dune::QuadratureRule::const_iterator it=rule.begin(); it!=rule.end(); ++it)
- {
- // evaluate reference element gradients of shape functions at quadrature point
- // (we assume Galerkin method lfsu=lfsv)
- std::vector vGradRef(dispSize);
- uScalarLFS.finiteElement().localBasis().evaluateJacobian(it->position(),vGradRef);
-
-
- // get inverse transposed jacobian for quadrature point
- const JacobianInverseTransposed jacobian = geometry.jacobianInverseTransposed(it->position());
-
- // calculate shape function gradients at the quadrature point in global coordinates. This is done
- // by multiplying the reference element shape functions with the inverse transposed jacobian
- std::vector > vGrad(dispSize);
- for (size_t i = 0; i < dispSize; i++)
- {
- vGrad[i] = 0.0;
- jacobian.umv(vGradRef[i][0],vGrad[i]);
- }
-
- // calculate the gradient of the solid displacement vector uGrad
- // x(uLFS,i) is the solid displacement entry of the solution vector
- // for element vertex i and coordinate direction coordDir
-
- Dune::FieldMatrix uGrad(0.0);
- for(int coordDir = 0; coordDir < dim; ++coordDir) {
- const DisplacementScalarLFS& uLFS = displacementLFS.child(coordDir);
- // compute gradient of u
- for (size_t i = 0; i < dispSize; i++)
- uGrad[coordDir].axpy(x(uLFS,i),vGrad[i]);
- }
- // calculate the strain tensor epsilon
- Dune::FieldMatrix epsilon;
- for(int i = 0; i < dim; ++i)
- for(int j = 0; j < dim; ++j)
- epsilon[i][j] = 0.5*(uGrad[i][j] + uGrad[j][i]);
-
- RF traceEpsilon = 0;
- for(int i = 0; i < dim; ++i)
- traceEpsilon += epsilon[i][i];
-
- // calculate the effective stress tensor effStress
- Dune::FieldMatrix effStress(0.0);
- for(int i = 0; i < dim; ++i)
- {
- effStress[i][i] = lambda*traceEpsilon;
- for(int j = 0; j < dim; ++j)
- effStress[i][j] += 2.0*mu*epsilon[i][j];
- }
-
- // retrieve the shape functions for interpolating the primary variables at the
- // current quadrature point
- std::vector q(pressSize);
- pressLFS.finiteElement().localBasis().evaluateFunction(it->position(),q);
-
- RT_P pw(0.0);
- RT_P sn(0.0);
- RT_P ux(0.0);
- RT_P uy(0.0);
- RT_P uz(0.0);
-
- // interpolate primary variables at current quadrature point
- for (size_t i = 0; i < pressLFS.size(); i++)
- {
- pw += x(pressLFS,i) * q[i];
- sn += x(satLFS,i) * q[i];
- ux += x(displacementLFS.child(0),i) * q[i];
- if (dim > 1)
- uy += x(displacementLFS.child(1),i) * q[i];
- if (dim > 2)
- uz += x(displacementLFS.child(2),i) * q[i];
- }
- RT_P sw = 1.0 - sn;
- RT_P pn = pw + MaterialLaw::pc(materialParams, sw);
- RT_P pEff;
-
- const GlobalPosition& globalPos = geometry.global(it->position());
-
- // calculate change in effective pressure with respect to initial conditions pInit (pInit is negativ)
- pEff = pw*sw + pn*sn + model_.problem().pInit(globalPos, it->position(), eg.entity());
- RF uDiv = traceEpsilon;
- RF porosityEff;
-
- // assume deformation induced porosity changes
- if(model_.problem().coupled() == true){
- if (porosity + uDiv < 1e-3*porosity){
- DUNE_THROW(NumericalProblem, "volume change too large");
- }
- else
- // this equation would be correct if the bulk volume could change (Vol_new = Vol_init * (1+div u)), however, we
- // have a constant bulk volume therefore we should apply phi_eff = phi_init + div u
- // but this causes convergence problems. Since div u is very small here the chosen relation is
- // assumed to be a good approximation
- porosityEff = (porosity + uDiv)/(1.0 + uDiv);
- }
- // neglect deformation induced porosity changes
- else
- porosityEff = porosity;
-
- // fill primary variable vector for current quadrature point
- PrimaryVariables primVars;
-
- primVars[wPhaseIdx] = pw;
- primVars[nPhaseIdx] = sn;
- primVars[Indices::uxIdx] = ux;
- if (dim > 1)
- primVars[Indices::uyIdx] = uy;
- if (dim > 2)
- primVars[Indices::uzIdx] = uz;
-
- VolumeVariables volVars;
- // evaluate volume variables for this quadrature point
- // NOTE: this overwrites the entries of the volumevariables of node 0
- // and can cause errors
- volVars.update(primVars, model_.problem(), eg.entity(), fvGeometry, 0, false);
-
- // calculate the density difference for the gravity term
- RF rhoDiff = volVars.density(nPhaseIdx) - volVars.density(wPhaseIdx);
-
- // geometric weight need for quadrature rule evaluation (numerical integration)
- RF qWeight = it->weight() * geometry.integrationElement(it->position());
-
- // evaluate basis functions
- std::vector vBasis(dispSize);
- displacementLFS.child(0).finiteElement().localBasis().evaluateFunction(it->position(), vBasis);
-
- for(int coordDir = 0; coordDir < dim; ++coordDir) {
- const DisplacementScalarLFS& uLFS = displacementLFS.child(coordDir);
- // assemble momentum balance equation
- for (size_t i = 0; i < dispSize; i++){
- // multiply effective stress with gradient of weighting function and geometric weight of quadrature rule
- Scalar tmp = (effStress[coordDir] * vGrad[i]) * qWeight;
- r.rawAccumulate(uLFS,i,tmp);
-
- // subtract effective pressure change contribution multiplied with gradient of weighting function
- // and geometric weight of quadrature rule (soil mechanics sign conventions, compressive stresses are negative)
- tmp = -(pEff * vGrad[i][coordDir]) * qWeight;
- r.rawAccumulate(uLFS,i,tmp);
-
- // evaluate gravity term (soil mechanics sign conventions, compressive stresses are negative)
- // multiplied with weighting function and geometric weight of quadrature rule.
- // This assumes that the solid phase density remains constant, that the changes in porosity are very small,
- // and that the density of the brine phase remains constant
- tmp = sn*porosityEff*rhoDiff*model_.problem().gravity()[coordDir]*vBasis[i]* qWeight;
- r.rawAccumulate(uLFS,i,tmp);
- }
- }
- }
- // include boundary conditions
- // iterate over element intersections of codim dim-1
- for (const auto& intersection : intersections(model_.problem().gridView(), eg.entity()))
- {
- // handle only faces on the boundary
- if (!intersection.boundary())
- continue;
-
- // select quadrature rule for intersection faces (dim-1)
- Dune::GeometryType gtface = intersection.geometryInInside().type();
- const Dune::QuadratureRule& faceRule = Dune::QuadratureRules::rule(gtface,qorder);
-
- // get face index of this intersection
- int fIdx = intersection.indexInInside();
- // get dimension of face
- const int dimIs = Dune::PDELab::IntersectionGeometry::Entity::Geometry::mydimension;
-
- // get reference element for intersection geometry (reference element for face if dim = 3)
- const Dune::ReferenceElement& refElement = Dune::ReferenceElements::general(geomType);
- // get reference element for edges of intersection geometry (reference element for edge if dim = 3), needed for Dirichlet BC
- const Dune::ReferenceElement &face_refElement =
- Dune::ReferenceElements::general(intersection.geometryInInside().type());
-
- // Treat Neumann boundary conditions
- // loop over quadrature points and integrate normal stress changes (traction changes)
- for (typename Dune::QuadratureRule::const_iterator it=faceRule.begin(); it!=faceRule.end(); ++it)
- {
- // position of quadrature point in local coordinates of element
- DimVector local = intersection.geometryInInside().global(it->position());
-
- GlobalPosition globalPos = geometry.global(local);
-
- // evaluate boundary condition type
- BoundaryTypes boundaryTypes;
- model_.problem().boundaryTypesAtPos(boundaryTypes, globalPos);
-
- // skip rest if we are on Dirichlet boundary
- if (!boundaryTypes.hasNeumann())
- continue;
-
- // evaluate basis functions of all all element vertices for quadrature point location "local"
- std::vector