diff --git a/configure.ac b/configure.ac
index 10f0a7672db3e0365bdef068fda32dd013480701..0a331b511cc1f56e0638aecdda36958971a1ae47 100644
--- a/configure.ac
+++ b/configure.ac
@@ -20,7 +20,8 @@ AC_CONFIG_FILES([dumux.pc
dumux/boxmodels/2p/Makefile
dumux/boxmodels/2p2c/Makefile
dumux/boxmodels/2p2cni/Makefile
- dumux/boxmodels/2pni/Makefile
+ dumux/boxmodels/2pni/Makefile
+ dumux/boxmodels/2pdfm/Makefile
dumux/boxmodels/3p3c/Makefile
dumux/boxmodels/3p3cni/Makefile
dumux/boxmodels/common/Makefile
diff --git a/dumux/boxmodels/2pdfm/2pdfmfluxvariables.hh b/dumux/boxmodels/2pdfm/2pdfmfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..93ca3137d6bba43107bc859f0656bdde8596549c
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmfluxvariables.hh
@@ -0,0 +1,234 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This file contains the data which is required to calculate
+ * all fluxes of fluid phases over a face of a finite volume in the
+ * two phase discrete fracture-matrix model.
+ */
+#ifndef DUMUX_BOXMODELS_2PDFM_FLUX_VARIABLES_HH
+#define DUMUX_BOXMODELS_2PDFM_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/boxmodels/common/boxdarcyfluxvariables.hh>
+
+#include "2pdfmproperties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPDFMBoxModel
+ * \ingroup BoxDFMFluxVariables
+ * \brief Contains the data which is required to calculate the fluxes of
+ * the fluid phases over a face of a finite volume for the two-phase
+ * discrete fracture-matrix model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the intergration point, etc.
+ */
+template <class TypeTag>
+class TwoPDFMFluxVariables : public BoxDarcyFluxVariables<TypeTag>
+{
+ 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,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldVector<Scalar, dimWorld> Vector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry in the box scheme
+ * \param faceIdx 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
+ */
+ TwoPDFMFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ int faceIdx,
+ const ElementVolumeVariables &elemVolVars,
+ const bool onBoundary = false)
+ : BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary),
+ fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
+ {
+ faceSCV_ = &this->face();
+
+ for (int phase = 0; phase < numPhases; ++phase)
+ {
+ potentialGradFracture_[phase] = 0.0;
+ }
+
+ calculateGradientsInFractures_(problem, element, elemVolVars, faceIdx);
+ calculateVelocitiesFracture_(problem, element, elemVolVars, faceIdx);
+ };
+
+public:
+ /*!
+ * \brief Calculates the velocities in the lower-dimenstional fracture.
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param elemVolVars The volume variables of the current element
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ */
+ void calculateVelocitiesFracture_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars,
+ int faceIdx)
+ {
+ isFracture_ = problem.spatialParams().isEdgeFracture(element, faceIdx);
+ fractureWidth_ = problem.spatialParams().fractureWidth(element, faceIdx);
+
+ Scalar KFracture, KFi, KFj; //permeabilities of the fracture
+ if (isFracture_)
+ {
+ KFi = problem.spatialParams().
+ intrinsicPermeabilityFracture(element, this->fvGeometry_, this->face().i);
+ KFj = problem.spatialParams().
+ intrinsicPermeabilityFracture(element, this->fvGeometry_, this->face().j);
+ }
+ else
+ {
+ KFi = 0;
+ KFj = 0;
+ }
+
+ KFracture = Dumux::harmonicMean(KFi, KFj);
+
+ // temporary vector for the Darcy velocity
+ Scalar vDarcyFracture = 0;
+ for (int phase=0; phase < numPhases; phase++)
+ {
+ vDarcyFracture = - (KFracture * potentialGradFracture_[phase]);
+
+ Valgrind::CheckDefined(KFracture);
+ Valgrind::CheckDefined(fractureWidth_);
+ vDarcyFracture_[phase] = (vDarcyFracture * fractureWidth_);
+ Valgrind::CheckDefined(vDarcyFracture_[phase]);
+ }
+
+ // set the upstream and downstream vertices
+ for (int phase = 0; phase < numPhases; ++phase)
+ {
+ upstreamFractureIdx[phase] = faceSCV_->i;
+ downstreamFractureIdx[phase] = faceSCV_->j;
+
+ if (vDarcyFracture_[phase] < 0)
+ {
+ std::swap(upstreamFractureIdx[phase],
+ downstreamFractureIdx[phase]);
+ }
+ }
+ }
+
+ /*!
+ * \brief Return the pressure potential gradient in the lower dimensional fracture.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ const Scalar &potentialGradFracture(int phaseIdx) const
+ {
+ return potentialGradFracture_[phaseIdx];
+ }
+
+
+ PhasesVector vDarcyFracture_;
+
+ int upstreamFractureIdx[numPhases];
+ int downstreamFractureIdx[numPhases];
+protected:
+ // gradients
+ Scalar potentialGradFracture_[numPhases];
+ const FVElementGeometry &fvGeometry_;
+ int faceIdx_;
+ const bool onBoundary_;
+ bool isFracture_;
+ Scalar fractureWidth_;
+ const SCVFace *faceSCV_;
+
+private:
+ /*!
+ * \brief Calculates the gradients in the lower-dimenstional fracture.
+ *
+ * \param problem The problem
+ * \param element The finite element
+ * \param elemVolVars The volume variables of the current element
+ * \param faceIdx The local index of the SCV (sub-control-volume) face
+ */
+ void calculateGradientsInFractures_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars,
+ int faceIdx)
+ {
+ // calculate gradients, loop over adjacent vertices
+ for (int idx = 0; idx < this->fvGeometry_.numFAP; idx++)
+ {
+ int i = this->face().i;
+ int j = this->face().j;
+
+ // compute sum of pressure gradients for each phaseGlobalPosition
+ for (int phase = 0; phase < numPhases; phase++)
+ {
+ const GlobalPosition localIdx_i = element.geometry().corner(i);
+ const GlobalPosition localIdx_j = element.geometry().corner(j);
+
+ isFracture_ = problem.spatialParams().isEdgeFracture(element, faceIdx);
+
+ if (isFracture_)
+ {
+ GlobalPosition diff_ij = localIdx_j;
+ diff_ij -= localIdx_i;
+ potentialGradFracture_[phase] =
+ (elemVolVars[j].pressure(phase) - elemVolVars[i].pressure(phase))
+ / diff_ij.two_norm();
+ }
+ else
+ {
+ potentialGradFracture_[phase] = 0;
+ }
+ }
+ }
+ }
+};
+
+} // end namepace
+
+#endif // DUMUX_BOXMODELS_2PDFM_FLUX_VARIABLES_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmindices.hh b/dumux/boxmodels/2pdfm/2pdfmindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6326b0508fa69ac17ec754eea8b63dffd54481a6
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmindices.hh
@@ -0,0 +1,64 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+ * \file
+ *
+ * \brief Defines the indices required for the finite volume in the
+ * two phase discrete fracture-matrix model.
+ */
+#ifndef DUMUX_BOXMODELS_2PDFM_INDICES_HH
+#define DUMUX_BOXMODELS_2PDFM_INDICES_HH
+
+#include<dumux/boxmodels/2p/2pindices.hh>
+
+namespace Dumux
+{
+// \{
+
+/*!
+ * \ingroup TwoPDFMBoxModel
+ * \ingroup BoxIndices
+ * \brief The common indices for the \f$p_w-S_n\f$ formulation of the
+ * isothermal two-phase discrete fracture-matrix model.
+ *
+ * \tparam TypeTag The problem type tag
+ * \tparam formulation The formulation, either pwSn or pnSw
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag,
+ int formulation = TwoPFormulation::pwSn,
+ int PVOffset = 0>
+struct TwoPDFMIndices : public TwoPIndices <TypeTag, formulation, PVOffset>
+{
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+ static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
+
+ // indices of the primary variables
+ static const int pwIdx = PVOffset + 0; //!< Pressure index of the wetting phase
+ static const int SnIdx = PVOffset + 1; //!< Saturation index of the wetting phase
+
+ // indices of the equations
+ static const int contiWEqIdx = PVOffset + 0; //!< Index of the continuity equation of the wetting phase
+ static const int contiNEqIdx = PVOffset + 1; //!< Index of the continuity equation of the non-wetting phase
+};
+
+// \}
+} // namespace Dumux
+
+#endif // DUMUX_BOXMODELS_2PDFM_INDICES_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmlocalresidual.hh b/dumux/boxmodels/2pdfm/2pdfmlocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3ab534a7b2bbc2b7514cd1fbf283f2c04f45fa04
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmlocalresidual.hh
@@ -0,0 +1,335 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Element-wise calculation of the residual for the finite volume in the
+ * two phase discrete fracture-matrix model.
+ */
+#ifndef DUMUX_BOXMODELS_2PDFM_LOCAL_RESIDUAL_HH
+#define DUMUX_BOXMODELS_2PDFM_LOCAL_RESIDUAL_HH
+
+#include <dumux/boxmodels/2p/2plocalresidual.hh>
+#include "2pdfmproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPDFMBoxModel
+ * \ingroup BoxLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase discrete fracture box model.
+ */
+template<class TypeTag>
+class TwoPDFMLocalResidual : public TwoPLocalResidual<TypeTag>
+{
+protected:
+ typedef TwoPLocalResidual<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum
+ {
+ contiWEqIdx = Indices::contiWEqIdx,
+ contiNEqIdx = Indices::contiNEqIdx,
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef Dune::UGGrid<2> GridType;
+ typedef typename GridType::ctype DT;
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldVector<Scalar, dimWorld> Vector;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief Constructor. Sets the upwind weight.
+ */
+ TwoPDFMLocalResidual()
+ {
+ // 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
+ massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+ };
+
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a finite sub-control volume.
+ *
+ * \param storage The phase mass within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \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.
+ ParentType::computeStorage(storage,scvIdx,usePrevSol);
+ computeStorageFracture(storage,scvIdx,usePrevSol);
+ }
+
+ /*!
+ * \brief Evaluate the storage term of the current solution in a
+ * lower-dimensional fracture.
+ *
+ * \param storage The phase mass within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \param usePrevSol Evaluate function with solution of current or previous time step
+ */
+ void computeStorageFracture(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
+ {
+ const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_()
+ : this->curVolVars_();
+ const VolumeVariables &vertDat = elemDat[scvIdx];
+ const Element &elem = this->element_();
+ bool isFracture = this->problem_().spatialParams().isVertexFracture(elem, scvIdx);
+ /*
+ * Sn = w_F * SnF + w_M * SnM
+ * First simple case before determining the real w_F is to assume that it is 0
+ * and w_M = 1
+ *
+ */
+ ///////////////////////////////////////////////////////////////////////
+ Scalar w_F, w_M; //volumetric fractions of fracture and matrix;
+ Scalar fractureVolume = 0.0;
+ w_F = 0.0;
+ /*
+ * Calculate the fracture volume fraction w_F = 0.5 * Fwidth * 0.5 * Length
+ */
+ Dune::GeometryType gt = elem.geometry().type();
+ const typename Dune::GenericReferenceElementContainer<DT,dim>::value_type&
+ refElem = Dune::GenericReferenceElements<DT,dim>::general(gt);
+
+ Scalar vol; //subcontrol volume
+ FVElementGeometry fvelem = this->fvGeometry_();
+ vol = fvelem.subContVol[scvIdx].volume;
+ for (int faceIdx=0; faceIdx<refElem.size(1); faceIdx++)
+ {
+ SCVFace face = fvelem.subContVolFace[faceIdx];
+ int i=face.i;
+ int j=face.j;
+
+ if (this->problem_().spatialParams().isEdgeFracture(elem, faceIdx)
+ && (i == scvIdx || j == scvIdx))
+ {
+ Scalar fracture_width = this->problem_().spatialParams().fractureWidth(elem, faceIdx);
+
+ const GlobalPosition global_i = elem.geometry().corner(i);
+ const GlobalPosition global_j = elem.geometry().corner(j);
+ GlobalPosition diff_ij = global_j;
+ diff_ij -=global_i;
+ Scalar fracture_length = 0.5*diff_ij.two_norm();
+ //half is taken from the other element
+ fractureVolume += 0.5 * fracture_length * fracture_width;
+ }
+ }
+ w_F = fractureVolume/vol;
+ w_M = 1-w_F;
+ ///////////////////////////////////////////////////////////////////////
+ Scalar storageFracture[numPhases];
+ Scalar storageMatrix [numPhases];
+ storageFracture[wPhaseIdx] = 0.0;
+ storageFracture[nPhaseIdx] = 0.0;
+ storageMatrix[wPhaseIdx] = 0.0;
+ storageMatrix[nPhaseIdx] = 0.0;
+ // const GlobalPosition &globalPos = elem.geometry().corner(scvIdx);
+
+ Scalar dSM_dSF = vertDat.dSM_dSF();
+ if (!this->problem_().useInterfaceCondition())
+ {
+ dSM_dSF = 1.0;
+ }
+
+ if (isFracture)
+ {
+ for (int phaseIdx = 0; phaseIdx<2; phaseIdx++)
+ {
+ storageFracture[phaseIdx] = vertDat.density(phaseIdx)
+ * vertDat.porosityFracture()
+ * w_F
+ * vertDat.saturationFracture(phaseIdx);
+ storageMatrix[phaseIdx] = vertDat.density(phaseIdx)
+ * vertDat.porosity()
+ * w_M
+ * dSM_dSF
+ * vertDat.saturationMatrix(phaseIdx);
+ }
+ }
+ else
+ {
+ for (int phaseIdx = 0; phaseIdx < 2; phaseIdx++)
+ {
+ storageFracture[phaseIdx] = 0.0;
+ storageMatrix[phaseIdx] = vertDat.density(phaseIdx)
+ * vertDat.porosity()
+ * vertDat.saturation(phaseIdx);
+ }
+
+ }
+ // wetting phase mass
+ storage[contiWEqIdx] = storageFracture[wPhaseIdx]
+ + storageMatrix[wPhaseIdx];
+
+ // non-wetting phase mass
+ storage[contiNEqIdx] = storageFracture[nPhaseIdx]
+ + storageMatrix[nPhaseIdx];
+ }
+
+ /*!
+ * \brief Evaluates the mass flux over a face of a sub-control
+ * volume.
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face for each phase
+ * \param faceIdx The index of the SCV face
+ * \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 faceIdx, const bool onBoundary=false) const
+ {
+ FluxVariables fluxVars(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ faceIdx,
+ this->curVolVars_(),
+ onBoundary);
+ flux = 0;
+ asImp_()->computeAdvectiveFlux(flux, fluxVars);
+ asImp_()->computeAdvectiveFluxFracture(flux, fluxVars);
+ asImp_()->computeDiffusiveFlux(flux, fluxVars);
+ asImp_()->computeDiffusiveFluxFracture(flux, fluxVars);
+ }
+
+ /*!
+ * \brief Adds the flux vector in the lower dimensional fracture to the
+ * flux vector over the face of a sub-control volume.
+ *
+ * This method is called by compute flux and is mainly there for
+ * derived models to ease adding equations selectively.
+ *
+ * \param flux The advective flux over the sub-control-volume face for each phase
+ * \param fluxVars The flux variables at the current SCV
+ */
+ void computeAdvectiveFluxFracture(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ {
+ ////////
+ // advective fluxes of all components in all phases
+ ////////
+ Vector tmpVec;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // calculate the flux in direction of the
+ // current fracture
+ //
+ // v = - (K grad p) * n
+ //
+ // (the minus comes from the Darcy law which states that
+ // the flux is from high to low pressure potentials.)
+
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &upFracture =
+ this->curVolVars_(fluxVars.upstreamFractureIdx[phaseIdx]);
+ const VolumeVariables &dnFracture =
+ this->curVolVars_(fluxVars.downstreamFractureIdx[phaseIdx]);
+
+ Scalar fractureFlux =
+ 0.5 * fluxVars.vDarcyFracture_[phaseIdx]
+ * ( massUpwindWeight_ // upstream vertex
+ * (upFracture.density(phaseIdx)
+ * upFracture.mobilityFracture(phaseIdx))
+ + (1 - massUpwindWeight_) // downstream vertex
+ * (dnFracture.density(phaseIdx)
+ * dnFracture.mobilityFracture(phaseIdx)));
+
+ flux[phaseIdx] += fractureFlux;
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive flux of the fracture to the flux vector over
+ * the face of a sub-control volume.
+ *
+ * \param flux The diffusive flux over the sub-control-volume face for each phase
+ * \param fluxData The flux variables at the current SCV
+ *
+ * It doesn't do anything in two-phase model but is used by the
+ * non-isothermal two-phase models to calculate diffusive heat
+ * fluxes
+ */
+ void computeDiffusiveFluxFracture(PrimaryVariables &flux, const FluxVariables &fluxData) const
+ {
+ // diffusive fluxes
+ flux += 0.0;
+ }
+
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param source The source/sink in the SCV for each phase
+ * \param scvIdx The index of the SCV
+ *
+ */
+ void computeSource(PrimaryVariables &source, int scvIdx) const
+ {
+ // retrieve the source term intrinsic to the problem
+ this->problem_().boxSDSource(source,
+ this->element_(),
+ this->fvGeometry_(),
+ scvIdx,
+ this->curVolVars_());
+ }
+
+protected:
+ Implementation *asImp_()
+ {
+ return static_cast<Implementation *> (this);
+ }
+ const Implementation *asImp_() const
+ {
+ return static_cast<const Implementation *> (this);
+ }
+
+private:
+ Scalar massUpwindWeight_;
+};
+
+} // end namespace
+
+#endif // DUMUX_BOXMODELS_2PDFM_LOCAL_RESIDUAL_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmmodel.hh b/dumux/boxmodels/2pdfm/2pdfmmodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..07b9c92fd06861973664b4a72c5781f91a712098
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmmodel.hh
@@ -0,0 +1,390 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+
+/*!
+* \file
+*
+* \brief Adaption of the box scheme to the two-phase flow in discrete fracture-matrix
+* model.
+*/
+
+#ifndef DUMUX_BOXMODELS_2PDFM_MODEL_HH
+#define DUMUX_BOXMODELS_2PDFM_MODEL_HH
+
+#include <dumux/boxmodels/2p/2pmodel.hh>
+
+#include "2pdfmproperties.hh"
+#include "2pdfmlocalresidual.hh"
+#include "2pdfmproblem.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPDFMBoxModel
+ * \brief A two-phase, isothermal flow model using the box scheme.
+ *
+ * This model implements two-phase flow of two immiscible fluids
+ * \f$\alpha \in \{ w, n \}\f$ using a standard multiphase Darcy
+ * approach as the equation for the conservation of momentum, i.e.
+ \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
+ \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right)
+ \f]
+ *
+ * By inserting this into the equation for the conservation of the
+ * phase mass, one gets
+ \f[
+ \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
+ -
+ \text{div} \left\{
+ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ \right\} - q_\alpha = 0 \;,
+ \f]
+ *
+ * These equations are discretized by a fully-coupled vertex centered finite volume
+ * (box) scheme as spatial and the implicit Euler method as time
+ * discretization.
+ *
+ * By using constitutive relations for the capillary pressure \f$p_c =
+ * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
+ * advantage of the fact that \f$S_w + S_n = 1\f$, the number of
+ * unknowns can be reduced to two. Currently the model supports
+ * choosing either \f$p_w\f$ and \f$S_n\f$ or \f$p_n\f$ and \f$S_w\f$
+ * as primary variables. The formulation which ought to be used can be
+ * specified by setting the <tt>Formulation</tt> property to either
+ * <tt>TwoPCommonIndices::pWsN</tt> or <tt>TwoPCommonIndices::pNsW</tt>. By
+ * default, the model uses \f$p_w\f$ and \f$S_n\f$.
+ */
+template<class TypeTag >
+class TwoPDFMModel : public TwoPModel<TypeTag>
+{
+ typedef TwoPModel<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ nPhaseIdx = Indices::nPhaseIdx,
+ wPhaseIdx = Indices::wPhaseIdx,
+ pressureIdx = Indices::pressureIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::ctype CoordScalar;
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief Returns the relative weight of a primary variable for
+ * calculating relative errors.
+ *
+ * \param globalIdx The global index of the vertex
+ * \param pvIdx The index of the primary variable
+ */
+ DUNE_DEPRECATED
+ Scalar primaryVarWeight(int globalIdx, int pvIdx) const
+ {
+ if (pressureIdx == pvIdx)
+ {
+ return std::min(10.0 / this->prevSol()[globalIdx][pvIdx], 1.0);
+ }
+
+ return 1.0;
+ }
+
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ *
+ * \param sol The global solution vector
+ * \param writer The writer for multi-file VTK datasets
+ */
+ template<class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer)
+ {
+ bool velocityOutput = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
+ typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
+
+ // get the number of degrees of freedom and the number of vertices
+ unsigned numDofs = this->numDofs();
+ unsigned numVertices = this->problem_().gridView().size(dim);
+
+ // velocity output currently only works for vertex data
+ if (numDofs != numVertices)
+ {
+ velocityOutput = false;
+ }
+
+ // create the required scalar fields
+ ScalarField *pW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *pN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *pC = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Sw = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Sn = writer.allocateManagedBuffer(numDofs);
+ ScalarField *rhoW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *rhoN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *mobW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *mobN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *poro = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Te = writer.allocateManagedBuffer(numDofs);
+ ScalarField *cellNum =writer.allocateManagedBuffer (numDofs);
+ VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
+
+ if(velocityOutput) // check if velocity output is demanded
+ {
+ // initialize velocity fields
+ for (unsigned int i = 0; i < numVertices; ++i)
+ {
+ (*velocityN)[i] = Scalar(0);
+ (*velocityW)[i] = Scalar(0);
+ (*cellNum)[i] = Scalar(0.0);
+ }
+ }
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField *rank = writer.allocateManagedBuffer(numElements);
+
+ FVElementGeometry fvGeometry;
+ VolumeVariables volVars;
+ ElementVolumeVariables elemVolVars;
+
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt)
+ {
+ if(velocityOutput && !elemIt->geometry().type().isCube()){
+ DUNE_THROW(Dune::InvalidStateException,
+ "Currently, velocity output only works for cubes. "
+ "Please set the EnableVelocityOutput property to false!");
+ }
+ int idx = this->elementMapper().map(*elemIt);
+ (*rank)[idx] = this->gridView_().comm().rank();
+
+ fvGeometry.update(this->gridView_(), *elemIt);
+
+ if (numDofs == numElements) // element data
+ {
+ volVars.update(sol[idx],
+ this->problem_(),
+ *elemIt,
+ fvGeometry,
+ /*scvIdx=*/0,
+ false);
+
+ (*pW)[idx] = volVars.pressure(wPhaseIdx);
+ (*pN)[idx] = volVars.pressure(nPhaseIdx);
+ (*pC)[idx] = volVars.capillaryPressure();
+ (*Sw)[idx] = volVars.saturation(wPhaseIdx);
+ (*Sn)[idx] = volVars.saturation(nPhaseIdx);
+ (*rhoW)[idx] = volVars.density(wPhaseIdx);
+ (*rhoN)[idx] = volVars.density(nPhaseIdx);
+ (*mobW)[idx] = volVars.mobility(wPhaseIdx);
+ (*mobN)[idx] = volVars.mobility(nPhaseIdx);
+ (*poro)[idx] = volVars.porosity();
+ (*Te)[idx] = volVars.temperature();
+ }
+ else // vertex data
+ {
+ int numVerts = elemIt->template count<dim> ();
+ for (int scvIdx = 0; scvIdx < numVerts; ++scvIdx)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvGeometry,
+ scvIdx,
+ false);
+
+ (*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
+ (*pN)[globalIdx] = volVars.pressure(nPhaseIdx);
+ (*pC)[globalIdx] = volVars.capillaryPressure();
+ (*Sw)[globalIdx] = volVars.saturation(wPhaseIdx);
+ (*Sn)[globalIdx] = volVars.saturation(nPhaseIdx);
+ (*rhoW)[globalIdx] = volVars.density(wPhaseIdx);
+ (*rhoN)[globalIdx] = volVars.density(nPhaseIdx);
+ (*mobW)[globalIdx] = volVars.mobility(wPhaseIdx);
+ (*mobN)[globalIdx] = volVars.mobility(nPhaseIdx);
+ (*poro)[globalIdx] = volVars.porosity();
+ (*Te)[globalIdx] = volVars.temperature();
+ if(velocityOutput)
+ {
+ (*cellNum)[globalIdx] += 1;
+ }
+ }
+
+ if(velocityOutput)
+ {
+ // calculate vertex velocities
+ GlobalPosition tmpVelocity[numPhases];
+
+ for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ tmpVelocity[phaseIdx] = Scalar(0.0);
+ }
+
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
+ SCVVelocities scvVelocityW(8), scvVelocityN(8);
+
+ scvVelocityW = 0;
+ scvVelocityN = 0;
+
+ ElementVolumeVariables elemVolVars;
+
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ false /* oldSol? */);
+
+ for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
+ {
+ FluxVariables fluxVars(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ faceIdx,
+ elemVolVars);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // local position of integration point
+ const Dune::FieldVector<Scalar, dim>& localPosIP = fvGeometry.subContVolFace[faceIdx].ipLocal;
+
+ // Transformation of the global normal vector to normal vector in the reference element
+ const typename Element::Geometry::JacobianTransposed& jacobianT1 =
+ elemIt->geometry().jacobianTransposed(localPosIP);
+
+ const GlobalPosition globalNormal = fluxVars.face().normal;
+ GlobalPosition localNormal;
+ jacobianT1.mv(globalNormal, localNormal);
+ // note only works for cubes
+ const Scalar localArea = pow(2,-(dim-1));
+
+ localNormal /= localNormal.two_norm();
+
+ // Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolumeface
+ // in the reference element.
+ PhasesVector q;
+ q[phaseIdx] = fluxVars.volumeFlux(phaseIdx) / localArea;
+
+ // transform the normal Darcy velocity into a vector
+ tmpVelocity[phaseIdx] = localNormal;
+ tmpVelocity[phaseIdx] *= q[phaseIdx];
+
+ if(phaseIdx == wPhaseIdx)
+ {
+ scvVelocityW[fluxVars.face().i] += tmpVelocity[phaseIdx];
+ scvVelocityW[fluxVars.face().j] += tmpVelocity[phaseIdx];
+ }
+ else if(phaseIdx == nPhaseIdx)
+ {
+ scvVelocityN[fluxVars.face().i] += tmpVelocity[phaseIdx];
+ scvVelocityN[fluxVars.face().j] += tmpVelocity[phaseIdx];
+ }
+ }
+ }
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+ const Dune::FieldVector<Scalar, dim> &localPos
+ = ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
+
+ // get the transposed Jacobian of the element mapping
+ const typename Element::Geometry::JacobianTransposed& jacobianT2
+ = elemIt->geometry().jacobianTransposed(localPos);
+
+ // transform vertex velocities from local to global coordinates
+ for (int i = 0; i < numVerts; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ // calculate the subcontrolvolume velocity by the Piola transformation
+ Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
+
+ jacobianT2.mtv(scvVelocityW[i], scvVelocity);
+ scvVelocity /= elemIt->geometry().integrationElement(localPos);
+ // add up the wetting phase subcontrolvolume velocities for each vertex
+ (*velocityW)[globalIdx] += scvVelocity;
+
+ jacobianT2.mtv(scvVelocityN[i], scvVelocity);
+ scvVelocity /= elemIt->geometry().integrationElement(localPos);
+ // add up the nonwetting phase subcontrolvolume velocities for each vertex
+ (*velocityN)[globalIdx] += scvVelocity;
+ }
+ }
+ }
+ }
+
+ if (numDofs == numElements) // element data
+ {
+ writer.attachCellData(*Sn, "Sn");
+ writer.attachCellData(*Sw, "Sw");
+ writer.attachCellData(*pN, "pn");
+ writer.attachCellData(*pW, "pw");
+ writer.attachCellData(*pC, "pc");
+ writer.attachCellData(*rhoW, "rhoW");
+ writer.attachCellData(*rhoN, "rhoN");
+ writer.attachCellData(*mobW, "mobW");
+ writer.attachCellData(*mobN, "mobN");
+ writer.attachCellData(*poro, "porosity");
+ writer.attachCellData(*Te, "temperature");
+ }
+ else // vertex data
+ {
+ writer.attachVertexData(*Sn, "Sn");
+ writer.attachVertexData(*Sw, "Sw");
+ writer.attachVertexData(*pN, "pn");
+ writer.attachVertexData(*pW, "pw");
+ writer.attachVertexData(*pC, "pc");
+ writer.attachVertexData(*rhoW, "rhoW");
+ writer.attachVertexData(*rhoN, "rhoN");
+ writer.attachVertexData(*mobW, "mobW");
+ writer.attachVertexData(*mobN, "mobN");
+ writer.attachVertexData(*poro, "porosity");
+ writer.attachVertexData(*Te, "temperature");
+ if(velocityOutput) // check if velocity output is demanded
+ {
+ // divide the vertex velocities by the number of adjacent scvs i.e. cells
+ for(unsigned int globalIdx = 0; globalIdx < numVertices; ++globalIdx)
+ {
+ (*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
+ (*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
+ }
+ writer.attachVertexData(*velocityW, "velocityW", dim);
+ writer.attachVertexData(*velocityN, "velocityN", dim);
+ }
+ }
+ writer.attachCellData(*rank, "process rank");
+ }
+};
+} // end namespace
+
+#include "2pdfmpropertydefaults.hh"
+
+#endif // DUMUX_BOXMODELS_2PDFM_MODEL_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmproblem.hh b/dumux/boxmodels/2pdfm/2pdfmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fd65e2b68dd2fc9909b991d4f7b4631367888dc1
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmproblem.hh
@@ -0,0 +1,79 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for all problems which use the two-phase DFM box model
+ */
+#ifndef DUMUX_BOXMODELS_2PDFM_PROBLEM_HH
+#define DUMUX_BOXMODELS_2PDFM_PROBLEM_HH
+
+#include <dumux/boxmodels/common/porousmediaboxproblem.hh>
+#include "2pdfmproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup BoxBaseProblems
+ * \ingroup TwoPBoxModel
+ * \brief Base class for all problems which use the two-phase box model
+ */
+template<class TypeTag>
+class TwoPDFMProblem : public PorousMediaBoxProblem<TypeTag>
+{
+ typedef PorousMediaBoxProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ * \param verbose Turn verbosity on or off
+ */
+ DUNE_DEPRECATED_MSG("use PorousMediaBoxProblem instead")
+ TwoPDFMProblem(TimeManager &timeManager,
+ const GridView &gridView,
+ bool verbose = true)
+ : ParentType(timeManager, gridView)
+ {}
+
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ * \param spatialParams The spatial parameters object
+ * \param verbose Turn verbosity on or off
+ */
+ TwoPDFMProblem(TimeManager &timeManager,
+ const GridView &gridView,
+ SpatialParams &spatialParams,
+ bool verbose = true)
+ : ParentType(timeManager, gridView)
+ {
+ this->newSpatialParams_ = false;
+ delete this->spatialParams_;
+ this->spatialParams_ = &spatialParams;
+ }
+};
+} // end namespace
+
+#endif // DUMUX_BOXMODELS_2PDFM_PROBLEM_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmproperties.hh b/dumux/boxmodels/2pdfm/2pdfmproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7d0d2631e056a07bd93b21f6ee3b1d5c4637b3be
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmproperties.hh
@@ -0,0 +1,69 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the two phase DFM BOX model.
+ *
+ * \ingroup Properties
+ * \ingroup BoxProperties
+ * \ingroup TwoPDFMBoxModel
+ */
+
+#ifndef DUMUX_BOXMODELS_2PDFM_PROPERTIES_HH
+#define DUMUX_BOXMODELS_2PDFM_PROPERTIES_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+
+namespace Dumux
+{
+
+////////////////////////////////
+// properties
+////////////////////////////////
+namespace Properties
+{
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for 2pDFM problems
+NEW_TYPE_TAG(BoxTwoPDFM, INHERITS_FROM(BoxModel));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the weight of the upwind direction in the mass conservation equations
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
+NEW_PROP_TAG(Formulation); //!< The formulation of the model
+NEW_PROP_TAG(TwoPDFMIndices); //!< Enumerations for the 2pDFM models
+NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
+NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The context material law (extracted from the spatial parameters)
+NEW_PROP_TAG(WettingPhase); //!< The wetting phase for two-phase models
+NEW_PROP_TAG(NonwettingPhase); //!< The non-wetting phase for two-phase models
+NEW_PROP_TAG(FluidSystem); //!<The fluid systems including the information about the phases
+NEW_PROP_TAG(FluidState); //!<The phases state
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+} // end namespace Properties
+} // end namespace Dumux
+
+#endif // DUMUX_BOXMODELS_2PDFM_PROPERTIES_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmpropertydefaults.hh b/dumux/boxmodels/2pdfm/2pdfmpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..84121c23b0f4ed368d4a5e8b78a18fde2697976c
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmpropertydefaults.hh
@@ -0,0 +1,144 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Defines default values for the properties required by the
+ * 2pDFM box model.
+ *
+ * \ingroup Properties
+ * \ingroup TwoPDFMBoxModel
+ * \ingroup BoxProperties
+ */
+#ifndef DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
+#define DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
+
+#include <dumux/boxmodels/common/boxdarcyfluxvariables.hh>
+#include <dumux/material/components/nullcomponent.hh>
+#include <dumux/material/fluidstates/immisciblefluidstate.hh>
+#include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/spatialparams/boxspatialparams.hh>
+
+#include "2pdfmmodel.hh"
+#include "2pdfmproblem.hh"
+#include "2pdfmindices.hh"
+#include "2pdfmfluxvariables.hh"
+#include "2pdfmvolumevariables.hh"
+#include "2pdfmproperties.hh"
+
+namespace Dumux
+{
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property defaults
+//////////////////////////////////////////////////////////////////
+SET_INT_PROP(BoxTwoPDFM, NumEq, 2); //!< set the number of equations to 2
+SET_INT_PROP(BoxTwoPDFM, NumPhases, 2); //!< The number of phases in the 2p model is 2
+
+//! Set the default formulation to pWsN
+SET_INT_PROP(BoxTwoPDFM,
+ Formulation,
+ TwoPFormulation::pwSn);
+
+//! Use the 2p local jacobian operator for the 2p model
+SET_TYPE_PROP(BoxTwoPDFM,
+ LocalResidual,
+ TwoPDFMLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxTwoPDFM, Model, TwoPDFMModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxTwoPDFM, VolumeVariables, TwoPDFMVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxTwoPDFM, FluxVariables, TwoPDFMFluxVariables<TypeTag>);
+
+//! the upwind weight for the mass conservation equations.
+SET_SCALAR_PROP(BoxTwoPDFM, ImplicitMassUpwindWeight, 1.0);
+
+//! weight for the upwind mobility in the velocity calculation
+SET_SCALAR_PROP(BoxTwoPDFM, ImplicitMobilityUpwindWeight, 1.0);
+
+//! The indices required by the isothermal 2pDFM model
+SET_PROP(BoxTwoPDFM, Indices)
+{ private:
+ enum { Formulation = GET_PROP_VALUE(TypeTag, Formulation) };
+ public:
+ typedef TwoPDFMIndices<TypeTag, Formulation, 0> type;
+};
+
+
+//! The spatial parameters to be employed.
+//! Use BoxSpatialParams by default.
+SET_TYPE_PROP(BoxTwoPDFM, SpatialParams, BoxSpatialParams<TypeTag>);
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_TYPE_PROP(BoxTwoPDFM,
+ MaterialLawParams,
+ typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+SET_PROP(BoxTwoPDFM, WettingPhase)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+};
+
+SET_PROP(BoxTwoPDFM, NonwettingPhase)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+};
+
+SET_PROP(BoxTwoPDFM, FluidSystem)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
+ typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
+
+public:
+ typedef Dumux::FluidSystems::TwoPImmiscible<Scalar,
+ WettingPhase,
+ NonwettingPhase> type;
+};
+
+SET_PROP(BoxTwoPDFM, FluidState)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+public:
+ typedef ImmiscibleFluidState<Scalar, FluidSystem> type;
+};
+
+// disable velocity output by default
+SET_BOOL_PROP(BoxTwoPDFM, VtkAddVelocity, false);
+
+// enable gravity by default
+SET_BOOL_PROP(BoxTwoPDFM, ProblemEnableGravity, true);
+} // end namespace Properties
+} // end namespace Dumux
+
+#endif // DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
diff --git a/dumux/boxmodels/2pdfm/2pdfmvolumevariables.hh b/dumux/boxmodels/2pdfm/2pdfmvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c49cce37b57b620697c56ae13bcdf399bce212ef
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/2pdfmvolumevariables.hh
@@ -0,0 +1,385 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the two-phase discrete fracture-matrix model.
+ */
+#ifndef DUMUX_BOXMODELS_2PDFM_VOLUME_VARIABLES_HH
+#define DUMUX_BOXMODELS_2PDFM_VOLUME_VARIABLES_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dumux/boxmodels/2p/2pvolumevariables.hh>
+
+#include "2pdfmproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPDFMBoxModel
+ * \ingroup BoxVolumeVariables
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the two-phase discrete fracture-matrix model.
+ */
+template <class TypeTag>
+class TwoPDFMVolumeVariables : public TwoPVolumeVariables<TypeTag>
+{
+ typedef TwoPVolumeVariables<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ pwSn = Indices::pwSn,
+ pnSw = Indices::pnSw,
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ formulation = GET_PROP_VALUE(TypeTag, Formulation)
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef Dune::UGGrid<2> GridType;
+ typedef typename GridType::ctype DT;
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+ /*!
+ * \copydoc BoxVolumeVariables::update
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ int scvIdx,
+ bool oldSol)
+ {
+ ParentType::update(priVars,
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ oldSol);
+
+ this->completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidState_);
+
+ const MaterialLawParams &materialParams =
+ problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+
+ mobilityMatrix_[wPhaseIdx] =
+ MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx))
+ / fluidState_.viscosity(wPhaseIdx);
+
+ mobilityMatrix_[nPhaseIdx] =
+ MaterialLaw::krn(materialParams, fluidState_.saturation(wPhaseIdx))
+ / fluidState_.viscosity(nPhaseIdx);
+
+ // porosity
+ porosityMatrix_ = problem.spatialParams().porosity(element, fvGeometry, scvIdx);
+
+ // energy related quantities not belonging to the fluid state
+ asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, oldSol);
+ asImp_().updateFracture(priVars, problem, element, fvGeometry, scvIdx, oldSol);
+ }
+
+ /*!
+ * \brief Construct the volume variables for all fracture vertices.
+ *
+ * \param priVars Primary variables.
+ * \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.
+ */
+ void updateFracture(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ int scvIdx,
+ bool oldSol)
+ {
+ PrimaryVariables varsFracture;
+ const MaterialLawParams &materialParamsMatrix =
+ problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+ Scalar pressure[numPhases];
+ Scalar pMatrix[numPhases];
+ Scalar pFract[numPhases];
+
+ satNMatrix_ = priVars[saturationIdx];
+ satWMatrix_ = 1.0 - satNMatrix_;
+ satN_ = satNMatrix_;
+ satW_ = satWMatrix_;
+
+ pCMatrix_ = MaterialLaw::pC(materialParamsMatrix, satWMatrix_);
+ pC_ = pCMatrix_;
+ //pressures
+ pMatrix[wPhaseIdx] = priVars[pressureIdx];
+ pMatrix[nPhaseIdx] = pMatrix[wPhaseIdx] + pCMatrix_;
+ //Initialize pFract with the same values as the ones in the matrix
+ pFract[wPhaseIdx] = pMatrix[wPhaseIdx];
+ pFract[nPhaseIdx] = satNMatrix_;
+
+ varsFracture[pressureIdx] = pFract[wPhaseIdx];
+ varsFracture[saturationIdx] = pFract[wPhaseIdx];
+
+ this->completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidStateFracture_);
+
+ //Checks if the node is on a fracture
+ isNodeOnFracture_ = problem.spatialParams().isVertexFracture(element, scvIdx);
+
+ ///////////////////////////////////////////////////////////////////////////////
+ if (isNodeOnFracture_)
+ {
+ const MaterialLawParams &materialParamsFracture =
+ problem.spatialParams().materialLawParamsFracture(element, fvGeometry, scvIdx);
+
+ satNFracture_ = priVars[saturationIdx];
+ satWFracture_ = 1 - satNFracture_;
+ pCFracture_ = MaterialLaw::pC(materialParamsFracture, satWFracture_);
+ pFract[wPhaseIdx] = priVars[pressureIdx];
+ pFract[nPhaseIdx] = pFract[wPhaseIdx] + pCFracture_;
+ pEntryMatrix_ = MaterialLaw::pC(materialParamsMatrix, 1);
+
+ //use interface condition - extended capillary pressure inteface condition
+ if (problem.useInterfaceCondition())
+ {
+ interfaceCondition(materialParamsMatrix);
+ }
+ pC_ = pCFracture_;
+ satW_ = satWFracture_; //for plotting we are interested in the saturations of the fracture
+ satN_ = satNFracture_;
+ mobilityFracture_[wPhaseIdx] =
+ MaterialLaw::krw(materialParamsFracture, fluidStateFracture_.saturation(wPhaseIdx))
+ / fluidStateFracture_.viscosity(wPhaseIdx);
+
+ mobilityFracture_[nPhaseIdx] =
+ MaterialLaw::krn(materialParamsFracture, fluidStateFracture_.saturation(wPhaseIdx))
+ / fluidStateFracture_.viscosity(nPhaseIdx);
+
+ // derivative resulted from BrooksCorey pc_Sw formulation
+ dSM_dSF_ = (1 - problem.spatialParams().SwrM_) / (1 - problem.spatialParams().SwrF_)
+ * pow((problem.spatialParams().pdM_/ problem.spatialParams().pdF_),problem.spatialParams().lambdaM_)
+ * (problem.spatialParams().lambdaM_ / problem.spatialParams().lambdaF_)
+ * pow((satWFracture_ - problem.spatialParams().SwrF_ ) / (1 - problem.spatialParams().SwrF_),
+ (problem.spatialParams().lambdaM_ / problem.spatialParams().lambdaF_) - 1);
+ }// end if (node)
+ ///////////////////////////////////////////////////////////////////////////////
+ else
+ {
+ /* the values of pressure and saturation of the fractures in the volumes where
+ there are no fracture are set unphysical*/
+ satNFracture_ = -1;
+ satWFracture_ = -1;
+ pCFracture_ = -1e100;
+ pFract[wPhaseIdx] = -1e100;
+ pFract[nPhaseIdx] = -1e100;
+ pEntryMatrix_ = -1e100;
+ mobilityFracture_[wPhaseIdx] = 0.0;
+ mobilityFracture_[nPhaseIdx] = 0.0;
+ }
+ ///////////////////////////////////////////////////////////////////////////////
+ pressure[wPhaseIdx] = priVars[pressureIdx];
+ pressure[nPhaseIdx] = pressure[wPhaseIdx] + pC_;
+
+ porosityFracture_ = problem.spatialParams().porosityFracture(element,
+ fvGeometry,
+ scvIdx);
+ }
+
+ /*!
+ * \brief Extended capillary pressure saturation interface condition
+ *
+ * \param materialParamsMatrix the material law o calculate the Sw as inverse of capillary pressure function
+ *
+ * This method is called by updateFracture
+ */
+ void interfaceCondition(const MaterialLawParams &materialParamsMatrix)
+ {
+ /*2nd condition Niessner, J., R. Helmig, H. Jakobs, and J.E. Roberts. 2005, eq.10
+ * if the capillary pressure in the fracture is smaller than the entry pressure
+ * in the matrix than in the matrix
+ * */
+ if (pCFracture_ <= pEntryMatrix_)
+ {
+ satWMatrix_ = 1.0;
+ satNMatrix_ = 1 - satWMatrix_;
+ }
+ //3rd condition Niessner, J., R. Helmig, H. Jakobs, and J.E. Roberts. 2005, eq.10
+ else
+ {
+ /*
+ * Inverse capillary pressure function SwM = pcM^(-1)(pcF(SwF))
+ */
+ satWMatrix_ = MaterialLaw::Sw(materialParamsMatrix, pCFracture_);
+ satNMatrix_ = 1 - satWMatrix_;
+ }
+ }
+
+ /*!
+ * \brief Calculates the volume of the fracture inside the SCV
+ */
+ Scalar calculateSCVFractureVolume ( const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ int scvIdx)
+ {
+ Scalar volSCVFracture;
+ Dune::GeometryType gt = element.geometry().type();
+ const typename Dune::GenericReferenceElementContainer<DT,dim>::value_type&
+ refElem = Dune::GenericReferenceElements<DT,dim>::general(gt);
+
+ for (int faceIdx=0; faceIdx<refElem.size(1); faceIdx++)
+ {
+ SCVFace face = fvGeometry.subContVolFace[faceIdx];
+ int i=face.i;
+ int j=face.j;
+
+ if (problem.spatialParams().isEdgeFracture(element, faceIdx)
+ && (i == scvIdx || j == scvIdx))
+ {
+ Scalar fracture_width = problem.spatialParams().fractureWidth();
+
+ const GlobalPosition global_i = element.geometry().corner(i);
+ const GlobalPosition global_j = element.geometry().corner(j);
+ GlobalPosition diff_ij = global_j;
+ diff_ij -=global_i;
+ //fracture length in the subcontrol volume is half d_ij
+ Scalar fracture_length = 0.5*diff_ij.two_norm();
+
+ volSCVFracture += 0.5 * fracture_length * fracture_width;
+ }
+ }
+ return volSCVFracture;
+ }
+
+ /*!
+ * \brief Returns the effective saturation fracture of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturationFracture(int phaseIdx) const
+ {
+ if (phaseIdx == wPhaseIdx)
+ return satWFracture_;
+ else
+ return satNFracture_;
+ }
+ Scalar saturationMatrix(int phaseIdx) const
+ {
+ if (phaseIdx == wPhaseIdx)
+ return satWMatrix_;
+ else
+ return satNMatrix_;
+ }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(int phaseIdx) const
+ { return mobilityMatrix_[phaseIdx]; }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobilityFracture(int phaseIdx) const
+ { return mobilityFracture_[phaseIdx]; }
+
+ /*!
+ * \brief Returns the average porosity within the matrix control volume.
+ */
+ Scalar porosity() const
+ { return porosityMatrix_; }
+
+ /*!
+ * \brief Returns the average porosity within the fracture.
+ */
+ Scalar porosityFracture() const
+ { return porosityFracture_; }
+
+ /*!
+ * \brief Returns the average permeability within the fracture.
+ */
+ Scalar permeabilityFracture() const
+ { return permeabilityFracture_; }
+
+ /*!
+ * \brief Returns the derivative dSM/dSF
+ */
+ Scalar dSM_dSF() const
+ { return dSM_dSF_;}
+
+protected:
+ FluidState fluidState_;
+ FluidState fluidStateFracture_;
+ Scalar porosityMatrix_;
+ Scalar porosityFracture_;
+ Scalar permeability_;
+ Scalar permeabilityFracture_;
+ Scalar mobilityMatrix_[numPhases];
+ Scalar mobilityFracture_[numPhases];
+
+ Scalar satW_;
+ Scalar satWFracture_;
+ Scalar satWMatrix_;
+ Scalar satN_;
+ Scalar satNFracture_;
+ Scalar satNMatrix_;
+
+ Scalar pC_;
+ Scalar pCFracture_;
+ Scalar pCMatrix_;
+ Scalar pEntryMatrix_;
+ Scalar dSM_dSF_;
+
+ bool isNodeOnFracture_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+} // end namespace
+
+#endif // DUMUX_BOXMODELS_2PDFM_VOLUME_VARIABLES_HH
diff --git a/dumux/boxmodels/2pdfm/Makefile.am b/dumux/boxmodels/2pdfm/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..117f16d6c2746e5d5f135c87626218ceadc23ba3
--- /dev/null
+++ b/dumux/boxmodels/2pdfm/Makefile.am
@@ -0,0 +1,4 @@
+2pdfmdir = $(includedir)/dumux/boxmodels/2pdfm
+2pdfm_HEADERS = *.hh
+
+include $(top_srcdir)/am/global-rules
diff --git a/dumux/boxmodels/Makefile.am b/dumux/boxmodels/Makefile.am
index eff7632454eb4a45a3797935a710d2ffefac5e94..0c3d2aaeab649b08a9cec100892c48f7fc0b7b5f 100644
--- a/dumux/boxmodels/Makefile.am
+++ b/dumux/boxmodels/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = common co2 co2ni 1p 1p2c 2p 2p2c 2p2cni 2pni 3p3c 3p3cni mpnc richards
+SUBDIRS = common co2 co2ni 1p 1p2c 2p 2p2c 2p2cni 2pni 2pdfm 3p3c 3p3cni mpnc richards
boxmodelsdir = $(includedir)/dumux/boxmodels