Commit 45de14f1 authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[common,boundarytypes,multidomain] Add new multidomainboundarytypes.

Small cleanup, removal of unused code and Mortar.
parent 0e56fefe
......@@ -23,8 +23,6 @@
#ifndef BOUNDARY_TYPES_HH
#define BOUNDARY_TYPES_HH
#include <dumux/common/valgrind.hh>
namespace Dumux
......@@ -45,8 +43,8 @@ public:
* \brief Reset the boundary types.
*
* After this method no equations will be disabled and neither
* neumann nor dirichlet conditions will be evaluated. This
* corrosponds to a neumann zero boundary.
* Neumann nor Dirichlet conditions will be evaluated. This
* corrosponds to a Neumann zero boundary.
*/
void reset()
{
......@@ -90,65 +88,42 @@ public:
};
/*!
* \brief Set all boundary conditions to neuman.
* \brief Set all boundary conditions to Neumann.
*/
void setAllNeumann()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
boundaryInfo_[eqIdx].visited = 1;
boundaryInfo_[eqIdx].isDirichlet = 0;
boundaryInfo_[eqIdx].isNeumann = 1;
boundaryInfo_[eqIdx].isOutflow = 0;
boundaryInfo_[eqIdx].isCouplingInflow = 0;
boundaryInfo_[eqIdx].isCouplingOutflow = 0;
boundaryInfo_[eqIdx].isMortarCoupling = 0;
Valgrind::SetDefined(boundaryInfo_[eqIdx]);
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
setNeumann(eqIdx);
}
}
/*!
* \brief Set all boundary conditions to dirichlet.
* \brief Set all boundary conditions to Dirichlet.
*/
void setAllDirichlet()
{
for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
boundaryInfo_[eqIdx].visited = 1;
boundaryInfo_[eqIdx].isDirichlet = 1;
boundaryInfo_[eqIdx].isNeumann = 0;
boundaryInfo_[eqIdx].isOutflow = 0;
boundaryInfo_[eqIdx].isCouplingInflow = 0;
boundaryInfo_[eqIdx].isCouplingOutflow = 0;
boundaryInfo_[eqIdx].isMortarCoupling = 0;
eq2pvIdx_[eqIdx] = eqIdx;
pv2eqIdx_[eqIdx] = eqIdx;
Valgrind::SetDefined(boundaryInfo_[eqIdx]);
for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
{
setDirichlet(eqIdx);
}
}
/*!
* \brief Set all boundary conditions to neuman.
* \brief Set all boundary conditions to Neumann.
*/
void setAllOutflow()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
boundaryInfo_[eqIdx].visited = 1;
boundaryInfo_[eqIdx].isDirichlet = 0;
boundaryInfo_[eqIdx].isNeumann = 0;
boundaryInfo_[eqIdx].isOutflow = 1;
boundaryInfo_[eqIdx].isCouplingInflow = 0;
boundaryInfo_[eqIdx].isCouplingOutflow = 0;
boundaryInfo_[eqIdx].isMortarCoupling = 0;
Valgrind::SetDefined(boundaryInfo_[eqIdx]);
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
setOutflow(eqIdx);
}
}
/*!
* \brief Set all boundary conditions to coupling inflow.
*/
DUNE_DEPRECATED_MSG("setAllCouplingInflow() is deprecated")
void setAllCouplingInflow()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
......@@ -167,6 +142,7 @@ public:
/*!
* \brief Set all boundary conditions to coupling outflow.
*/
DUNE_DEPRECATED_MSG("setAllCouplingOutflow() is deprecated")
void setAllCouplingOutflow()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
......@@ -185,6 +161,7 @@ public:
/*!
* \brief Set all boundary conditions to mortar coupling.
*/
DUNE_DEPRECATED_MSG("setAllMortarCoupling() is deprecated")
void setAllMortarCoupling()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
......@@ -201,7 +178,7 @@ public:
}
/*!
* \brief Set a neumann boundary condition for a single a single
* \brief Set a Neumann boundary condition for a single a single
* equation.
*
* \param eqIdx The index of the equation
......@@ -220,7 +197,7 @@ public:
}
/*!
* \brief Set a dirichlet boundary condition for a single primary
* \brief Set a Dirichlet boundary condition for a single primary
* variable
*
* \param pvIdx The index of the primary variable for which the
......@@ -246,7 +223,7 @@ public:
}
/*!
* \brief Set a neumann boundary condition for a single a single
* \brief Set a Neumann boundary condition for a single a single
* equation.
*
* \param eqIdx The index of the equation on which the outflow
......@@ -268,6 +245,7 @@ public:
/*!
* \brief Set a boundary condition for a single equation to coupling inflow.
*/
DUNE_DEPRECATED_MSG("setCouplingInflow() is deprecated")
void setCouplingInflow(int eqIdx)
{
boundaryInfo_[eqIdx].visited = 1;
......@@ -284,6 +262,7 @@ public:
/*!
* \brief Set a boundary condition for a single equation to coupling outflow.
*/
DUNE_DEPRECATED_MSG("setCouplingOutflow() is deprecated")
void setCouplingOutflow(int eqIdx)
{
boundaryInfo_[eqIdx].visited = 1;
......@@ -300,6 +279,7 @@ public:
/*!
* \brief Set a boundary condition for a single equation to mortar coupling.
*/
DUNE_DEPRECATED_MSG("setMortarCoupling() is deprecated")
void setMortarCoupling(int eqIdx)
{
boundaryInfo_[eqIdx].visited = 1;
......@@ -314,7 +294,7 @@ public:
}
/*!
* \brief Set a dirichlet boundary condition for a single primary
* \brief Set a Dirichlet boundary condition for a single primary
* variable.
*
* Depending on the discretization, setting the Dirichlet condition
......@@ -330,7 +310,7 @@ public:
/*!
* \brief Returns true if an equation is used to specify a
* dirichlet condition.
* Dirichlet condition.
*
* \param eqIdx The index of the equation
*/
......@@ -339,7 +319,7 @@ public:
/*!
* \brief Returns true if some equation is used to specify a
* dirichlet condition.
* Dirichlet condition.
*/
bool hasDirichlet() const
{
......@@ -351,7 +331,7 @@ public:
/*!
* \brief Returns true if an equation is used to specify a
* neumann condition.
* Neumann condition.
*
* \param eqIdx The index of the equation
*/
......@@ -360,7 +340,7 @@ public:
/*!
* \brief Returns true if some equation is used to specify a
* neumann condition.
* Neumann condition.
*/
bool hasNeumann() const
{
......@@ -397,6 +377,7 @@ public:
*
* \param eqIdx The index of the equation
*/
DUNE_DEPRECATED_MSG("isCouplingInflow() is deprecated")
bool isCouplingInflow(unsigned eqIdx) const
{ return boundaryInfo_[eqIdx].isCouplingInflow; }
......@@ -404,6 +385,7 @@ public:
* \brief Returns true if some equation is used to specify an
* inflow coupling condition.
*/
DUNE_DEPRECATED_MSG("hasCouplingInflow() is deprecated")
bool hasCouplingInflow() const
{
for (int i = 0; i < numEq; ++i)
......@@ -418,6 +400,7 @@ public:
*
* \param eqIdx The index of the equation
*/
DUNE_DEPRECATED_MSG("isCouplingOutflow() is deprecated")
bool isCouplingOutflow(unsigned eqIdx) const
{ return boundaryInfo_[eqIdx].isCouplingOutflow; }
......@@ -425,6 +408,7 @@ public:
* \brief Returns true if some equation is used to specify an
* outflow coupling condition.
*/
DUNE_DEPRECATED_MSG("hasCouplingOutflow() is deprecated")
bool hasCouplingOutflow() const
{
for (int i = 0; i < numEq; ++i)
......@@ -439,6 +423,7 @@ public:
*
* \param eqIdx The index of the equation
*/
DUNE_DEPRECATED_MSG("isMortarCoupling() is deprecated")
bool isMortarCoupling(unsigned eqIdx) const
{
return boundaryInfo_[eqIdx].isMortarCoupling;
......@@ -448,6 +433,7 @@ public:
* \brief Returns true if some equation is used to specify a
* mortar coupling condition.
*/
DUNE_DEPRECATED_MSG("hasMortarCoupling() is deprecated")
bool hasMortarCoupling() const
{
for (int i = 0; i < numEq; ++i)
......@@ -456,6 +442,28 @@ public:
return false;
}
/*!
* \brief Returns true if an equation is used to specify a
* coupling condition.
*
* \param eqIdx The index of the equation
*
* This only used for correcting the pressure in the Stokes equation.
*/
bool isCoupling(unsigned eqIdx) const
{
return false;
}
/*!
* \brief Returns true if some equation is used to specify a
* coupling condition.
*/
bool hasCoupling() const
{
return false;
}
/*!
* \brief Returns the index of the equation which should be used
* for the Dirichlet condition of the pvIdx's primary
......
......@@ -624,7 +624,7 @@ protected:
if (dim == 3)
DUNE_THROW(Dune::NotImplemented, "The function interpolateCornerPoints_() is only implemented for 2D.");
if (bcTypes.isCouplingInflow(massBalanceIdx) || bcTypes.isCouplingOutflow(massBalanceIdx))
if (bcTypes.isCoupling(massBalanceIdx))
{
if (scvIdx == 0 || scvIdx == 3)
this->residual_[scvIdx][massBalanceIdx] =
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Class to specify the type of a boundary for multidomain problems.
*/
#ifndef MULTIDOMAIN_BOUNDARY_TYPES_HH
#define MULTIDOMAIN_BOUNDARY_TYPES_HH
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/valgrind.hh>
namespace Dumux
{
/*!
* \ingroup BC
* \brief Class to specify the type of a boundary for multidomain problems.
*/
template <int numEq>
class MultidomainBoundaryTypes : public BoundaryTypes<numEq>
{
typedef BoundaryTypes<numEq> ParentType;
public:
MultidomainBoundaryTypes()
{ reset(); }
/*!
* \brief Reset the boundary types.
*
* After this method no equations will be disabled and neither
* neumann nor dirichlet conditions will be evaluated. This
* corrosponds to a neumann zero boundary.
*/
void reset()
{
ParentType::reset();
for (int i=0; i < numEq; ++i) {
boundaryCouplingInfo_[i].visited = 0;
boundaryCouplingInfo_[i].isCouplingDirichlet = 0;
boundaryCouplingInfo_[i].isCouplingNeumann = 0;
}
}
/*!
* \brief Set all boundary conditions to coupling Dirichlet.
*/
void setAllCouplingDirichlet()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
setCouplingDirichlet(eqIdx);
}
}
/*!
* \brief Set all boundary conditions to coupling Neumann.
*/
void setAllCouplingNeumann()
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
setCouplingNeumann(eqIdx);
}
}
/*!
* \brief Set a boundary condition for a single equation to coupling inflow.
*/
void setCouplingDirichlet(int eqIdx)
{
boundaryCouplingInfo_[eqIdx].visited = 1;
boundaryCouplingInfo_[eqIdx].isCouplingDirichlet = 1;
boundaryCouplingInfo_[eqIdx].isCouplingNeumann = 0;
Valgrind::SetDefined(boundaryCouplingInfo_[eqIdx]);
}
/*!
* \brief Set a boundary condition for a single equation to coupling outflow.
*/
void setCouplingNeumann(int eqIdx)
{
boundaryCouplingInfo_[eqIdx].visited = 1;
boundaryCouplingInfo_[eqIdx].isCouplingDirichlet = 0;
boundaryCouplingInfo_[eqIdx].isCouplingNeumann = 1;
Valgrind::SetDefined(boundaryCouplingInfo_[eqIdx]);
}
/*!
* \brief Returns true if an equation is used to specify an
* Dirichlet coupling condition.
*
* \param eqIdx The index of the equation
*/
bool isCouplingDirichlet(unsigned eqIdx) const
{ return boundaryCouplingInfo_[eqIdx].isCouplingDirichlet; }
/*!
* \brief Returns true if some equation is used to specify an
* Dirichlet coupling condition.
*/
bool hasCouplingDirichlet() const
{
for (int i = 0; i < numEq; ++i)
if (boundaryCouplingInfo_[i].isCouplingDirichlet)
return true;
return false;
}
/*!
* \brief Returns true if an equation is used to specify an
* Neumann coupling condition.
*
* \param eqIdx The index of the equation
*/
bool isCouplingNeumann(unsigned eqIdx) const
{ return boundaryCouplingInfo_[eqIdx].isCouplingNeumann; }
/*!
* \brief Returns true if some equation is used to specify an
* Neumann coupling condition.
*/
bool hasCouplingNeumann() const
{
for (int i = 0; i < numEq; ++i)
if (boundaryCouplingInfo_[i].isCouplingNeumann)
return true;
return false;
}
/*!
* \brief Returns true if an equation is used to specify a
* Mortar coupling condition.
*
* \param eqIdx The index of the equation
*/
bool isCouplingMortar(unsigned eqIdx) const
{
return false;
}
/*!
* \brief Returns true if some equation is used to specify an
* Mortar coupling condition.
*/
bool hasCouplingMortar() const
{
return false;
}
/*!
* \brief Returns true if an equation is used to specify a
* coupling condition.
*
* \param eqIdx The index of the equation
*/
bool isCoupling(unsigned eqIdx) const
{
return boundaryCouplingInfo_[eqIdx].isCouplingDirichlet
|| boundaryCouplingInfo_[eqIdx].isCouplingNeumann;
}
/*!
* \brief Returns true if some equation is used to specify an
* coupling condition.
*/
bool hasCoupling() const
{
for (int i = 0; i < numEq; ++i)
if (isCoupling(i))
return true;
return false;
}
private:
// this is a bitfield structure!
struct __attribute__((__packed__)) {
unsigned char visited : 1;
unsigned char isCouplingDirichlet : 1;
unsigned char isCouplingNeumann : 1;
} boundaryCouplingInfo_[numEq];
};
}
#endif
......@@ -47,8 +47,8 @@
#include "multidomainnewtoncontroller.hh"
#include "splitandmerge.hh"
#include <dumux/nonlinear/newtonmethod.hh>
#include <dumux/common/timemanager.hh>
#include <dumux/nonlinear/newtonmethod.hh>
namespace Dumux
{
......@@ -215,7 +215,8 @@ public:
};
SET_PROP(MultiDomain, NumEq1)
{ private:
{
private:
typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) TypeTag1;
enum {numEq = GET_PROP_VALUE(TypeTag1, NumEq)};
public:
......@@ -223,7 +224,8 @@ public:
};
SET_PROP(MultiDomain, NumEq2)
{ private:
{
private:
typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) TypeTag2;
enum {numEq = GET_PROP_VALUE(TypeTag2, NumEq)};
public:
......
......@@ -33,6 +33,7 @@
#include <dune/pdelab/constraints/conforming.hh>
#include "subdomainproperties.hh"
#include "multidomainboundarytypes.hh"
#include "multidomainproperties.hh"
#include "multidomainlocaloperator.hh"
#include <dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh>
......@@ -132,6 +133,10 @@ public:
{ return ParameterTree::unusedNewRunTimeParams(); }
};
//! Boundary types at a single degree of freedom
SET_TYPE_PROP(SubDomain, BoundaryTypes,
MultidomainBoundaryTypes<GET_PROP_VALUE(TypeTag, NumEq)>);
} // namespace Properties
} // namespace Dumux
......
......@@ -107,14 +107,13 @@ public:
// loop over the single vertices on the current face
for (int faceVertexIdx = 0; faceVertexIdx < numFaceVertices; ++faceVertexIdx)
{
const int boundaryFaceIdx = this->fvGeometry_().boundaryFaceIndex(fIdx, faceVertexIdx);
const int vIdx = refElement.subEntity(fIdx, 1, faceVertexIdx, dim);
// only evaluate, if we consider the same face vertex as in the outer
// loop over the element vertices
if (vIdx != idx)
continue;
if (boundaryHasCoupling_(this->bcTypes_(idx)))
if (this->bcTypes_(idx).hasCoupling())
evalCouplingVertex_(idx);
}
}
......@@ -261,6 +260,7 @@ public:
/*!
* \brief Check if one of the boundary conditions is coupling.
*/
DUNE_DEPRECATED_MSG("boundaryHasCoupling_() is unused in dumux and therefore deprecated")
bool boundaryHasCoupling_(const BoundaryTypes& bcTypes) const
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
......@@ -272,6 +272,7 @@ public:
/*!
* \brief Check if one of the boundary conditions is mortar coupling.
*/
DUNE_DEPRECATED_MSG("boundaryHasMortarCoupling_() is unused in dumux and therefore deprecated")
bool boundaryHasMortarCoupling_(const BoundaryTypes& bcTypes) const
{
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
......
......@@ -174,41 +174,8 @@ public:
asImp_()->evalNeumannSegment_(&intersection, idx, boundaryFaceIdx, boundaryVars);
asImp_()->evalOutflowSegment_(&intersection, idx, boundaryFaceIdx, boundaryVars);
//for the corner points, the boundary flux across the vertical non-coupling boundary faces
//has to be calculated to fulfill the mass balance
//convert suddomain intersection into multidomain intersection and check whether it is an outer boundary
if(!GridView::Grid::multiDomainIntersection(intersection).neighbor()
&& this->boundaryHasMortarCoupling_(this->bcTypes_(idx)))
{
const GlobalPosition& globalPos = this->fvGeometry_().subContVol[idx].global;
//problem specific function, in problem orientation of interface is known
if(this->problem_().isInterfaceCornerPoint(globalPos))
{
PrimaryVariables priVars(0.0);
// DimVector faceCoord = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
// std::cout<<faceCoord<<std::endl;
const int numVertices = refElement.size(dim);
bool evalBoundaryFlux = false;
for(int equationIdx = 0; equationIdx < numEq; ++equationIdx)
{
for(int i= 0; i < numVertices; i++)
{
//if vertex is on boundary and not the coupling vertex: check whether an outflow condition is set
if(this->model_().onBoundary(this->element_(), i) && i!=idx)
if (!this->bcTypes_(i).isOutflow(equationIdx))
evalBoundaryFlux = true;
}
//calculate the actual boundary fluxes and add to residual (only for momentum and transport equation, mass balance already has outflow)
asImp_()->computeFlux(priVars, boundaryFaceIdx, true/*on boundary*/);
if(evalBoundaryFlux)
this->residual_[idx][equationIdx] += priVars[equationIdx];
}
}
}
// Beavers-Joseph condition at the coupling boundary/interface
if(boundaryHasCoupling_(bcTypes) || boundaryHasMortarCoupling_(bcTypes))
if(boundaryHasCoupling_(bcTypes) || bcTypes.hasCouplingMortar())
{
evalBeaversJoseph_(&intersection, idx, boundaryFaceIdx, boundaryVars);