diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 004bb7e8f09c4ce0c77f1afd1dad6c11d7e3f100..37e51257d379532d0b23debe41745d4af85c922d 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -63,8 +63,6 @@ class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Imp
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
- static_assert(!Assembler::Problem::enableInternalDirichletConstraints(), "Internal Dirichlet constraints are currently not implemented for cc-methods!");
-
public:
using ParentType::ParentType;
@@ -141,6 +139,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/tru
using FVElementGeometry = typename GridGeometry::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
enum { dim = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension };
@@ -275,15 +274,54 @@ public:
// add the current partial derivatives to the global jacobian matrix
// no special treatment is needed if globalJ is a ghost because then derivatives have been assembled to 0 above
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ if constexpr (Problem::enableInternalDirichletConstraints())
{
- // the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraintsOwnElement = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsOwnElement[eqIdx])
+ {
+ origResiduals[0][eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ }
// off-diagonal entries
j = 1;
for (const auto& dataJ : connectivityMap[globalI])
- A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ {
+ const auto& neighborElement = neighborElements[j-1];
+ const auto& neighborScv = fvGeometry.scv(dataJ.globalJ);
+ const auto internalDirichletConstraintsNeighbor = this->problem().hasInternalDirichletConstraint(neighborElement, neighborScv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsNeighbor[eqIdx])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] = 0.0;
+ else
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j][eqIdx];
+ }
+
+ ++j;
+ }
+ }
+ else // no internal Dirichlet constraints specified
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // the diagonal entries
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+
+ // off-diagonal entries
+ j = 1;
+ for (const auto& dataJ : connectivityMap[globalI])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ }
}
// restore the original state of the scv's volume variables
@@ -325,6 +363,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/fal
using Element = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::template Codim<0>::Entity;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
@@ -343,7 +382,7 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
// assemble the undeflected residual
- const auto residual = this->evalLocalResidual()[0];
+ auto residual = this->evalLocalResidual()[0];
const auto storageResidual = this->evalLocalStorageResidual()[0];
//////////////////////////////////////////////////////////////////////////////////////////////////
@@ -405,8 +444,28 @@ public:
// add the current partial derivatives to the global jacobian matrix
// only diagonal entries for explicit jacobians
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
+ }
+ else
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
// restore the original state of the scv's volume variables
curVolVars = origVolVars;
@@ -435,6 +494,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/tr
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
@@ -461,7 +521,7 @@ public:
}
// assemble the undeflected residual
- const auto residual = this->evalLocalResidual()[0];
+ auto residual = this->evalLocalResidual()[0];
// get some aliases for convenience
const auto& problem = this->problem();
@@ -507,6 +567,30 @@ public:
}
}
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+
+ // inner faces
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (!scvf.boundary())
+ A[globalI][fvGeometry.scv(scvf.outsideScvIdx()).dofIndex()][eqIdx][pvIdx] = 0.0;
+ }
+ }
+ }
+ }
+
// return element residual
return residual;
}
@@ -527,6 +611,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/fa
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
@@ -563,6 +648,25 @@ public:
// add hand-code derivative of storage term
this->localResidual().addStorageDerivatives(A[globalI][globalI], this->problem(), this->element(), this->fvGeometry(), volVars, scv);
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ }
+ }
+ }
+
// return the original residual
return residual;
}
diff --git a/dumux/assembly/fvlocalassemblerbase.hh b/dumux/assembly/fvlocalassemblerbase.hh
index 7d894ca276fce2b65e621189b2ad167908ed3bbd..5b77bb610671ad2d7cdb5be6e00f448e8b68fb81 100644
--- a/dumux/assembly/fvlocalassemblerbase.hh
+++ b/dumux/assembly/fvlocalassemblerbase.hh
@@ -221,12 +221,14 @@ public:
// and set the residual to (privar - dirichletvalue)
for (const auto& scvI : scvs(this->fvGeometry()))
{
- if (this->problem().hasInternalDirichletConstraint(this->element(), scvI))
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scvI);
+ if (internalDirichletConstraints.any())
{
const auto dirichletValues = this->problem().internalDirichlet(this->element(), scvI);
// set the Dirichlet conditions in residual and jacobian
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- applyDirichlet(scvI, dirichletValues, eqIdx, eqIdx);
+ for (int eqIdx = 0; eqIdx < internalDirichletConstraints.size(); ++eqIdx)
+ if (internalDirichletConstraints[eqIdx])
+ applyDirichlet(scvI, dirichletValues, eqIdx, eqIdx);
}
}
}
diff --git a/dumux/multidomain/newtonsolver.hh b/dumux/multidomain/newtonsolver.hh
index 5d555124041799c3fcadaad06b72f36dba872df1..dba8c58a41b89b18c546275b1fb0fe29dda0235e 100644
--- a/dumux/multidomain/newtonsolver.hh
+++ b/dumux/multidomain/newtonsolver.hh
@@ -167,7 +167,7 @@ private:
const auto& problem = this->assembler().problem(id);
const auto& gridGeometry = this->assembler().gridGeometry(id);
auto& gridVariables = this->assembler().gridVariables(id);
- priVarSwitch.updateBoundary(problem, gridGeometry, gridVariables, sol[id]);
+ priVarSwitch.updateDirichletConstraints(problem, gridGeometry, gridVariables, sol[id]);
}
/*!
diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh
index 60eb4e8ac4d718045b873f727b1ce93201b127db..a4657ac225a24e6c0027d334a28d44ba2d3095de 100644
--- a/dumux/multidomain/subdomaincclocalassembler.hh
+++ b/dumux/multidomain/subdomaincclocalassembler.hh
@@ -84,8 +84,6 @@ class SubDomainCCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assem
using CouplingManager = typename Assembler::CouplingManager;
using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
- static_assert(!Problem::enableInternalDirichletConstraints(), "Internal Dirichlet constraints are currently not implemented for cc-methods!");
-
public:
//! export the domain id of this sub-domain
static constexpr auto domainId = typename Dune::index_constant<id>();
@@ -282,6 +280,7 @@ class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*i
{
using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>;
using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/true>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
@@ -406,15 +405,54 @@ public:
eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
// add the current partial derivatives to the global jacobian matrix
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ if constexpr (Problem::enableInternalDirichletConstraints())
{
- // the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraintsOwnElement = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsOwnElement[eqIdx])
+ {
+ origResiduals[0][eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ }
// off-diagonal entries
j = 1;
for (const auto& dataJ : connectivityMap[globalI])
- A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ {
+ const auto& neighborElement = neighborElements[j-1];
+ const auto& neighborScv = fvGeometry.scv(dataJ.globalJ);
+ const auto internalDirichletConstraintsNeighbor = this->problem().hasInternalDirichletConstraint(neighborElement, neighborScv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsNeighbor[eqIdx])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] = 0.0;
+ else
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j][eqIdx];
+ }
+
+ ++j;
+ }
+ }
+ else // no internal Dirichlet constraints specified
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // the diagonal entries
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+
+ // off-diagonal entries
+ j = 1;
+ for (const auto& dataJ : connectivityMap[globalI])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ }
}
// restore the original state of the scv's volume variables
@@ -526,8 +564,31 @@ public:
epsCoupl(priVarsJ[pvIdx], pvIdx), numDiffMethod);
// add the current partial derivatives to the global jacobian matrix
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ const auto& scv = fvGeometry.scv(globalI);
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ if (internalDirichletConstraints.none())
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
+ else
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ A[globalI][globalJ][eqIdx][pvIdx] = 0.0;
+ else
+ A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
+ }
+ }
+ else
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
// restore the current element solution
priVarsJ[pvIdx] = origPriVarsJ[pvIdx];
@@ -560,6 +621,7 @@ class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*i
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
static constexpr auto domainI = Dune::index_constant<id>();
@@ -645,8 +707,28 @@ public:
// add the current partial derivatives to the global jacobian matrix
// only diagonal entries for explicit jacobians
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
+ }
+ else
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
// restore the original state of the scv's volume variables
curVolVars = origVolVars;
@@ -687,6 +769,7 @@ class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*
using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
enum { dim = GridView::dimension };
@@ -749,8 +832,35 @@ public:
}
}
- // return element residual
- return this->evalLocalResidual()[0];
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ auto residual = this->evalLocalResidual()[0];
+
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+
+ // inner faces
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (!scvf.boundary())
+ A[globalI][fvGeometry.scv(scvf.outsideScvIdx()).dofIndex()][eqIdx][pvIdx] = 0.0;
+ }
+ }
+ }
+ return residual;
+ }
+ else
+ // return element residual
+ return this->evalLocalResidual()[0];
}
/*!
@@ -782,6 +892,20 @@ public:
{
const auto& elementJ = this->assembler().gridGeometry(domainJ).element(globalJ);
this->couplingManager().addCouplingDerivatives(A[globalI][globalJ], domainI, element, fvGeometry, curElemVolVars, domainJ, elementJ);
+
+ // add the current partial derivatives to the global jacobian matrix
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ const auto& scv = fvGeometry.scv(globalI);
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ if (internalDirichletConstraints.any())
+ {
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (internalDirichletConstraints[eqIdx])
+ A[globalI][globalJ][eqIdx][pvIdx] = 0.0;
+ }
+ }
}
}
};
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 261d2434d61efef7792b0ed674ac1783dfd49e52..a5e34418f45e7692e4b1c29d1b803a05cbc78074 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -808,7 +808,7 @@ protected:
const auto& problem = this->assembler().problem();
const auto& gridGeometry = this->assembler().gridGeometry();
auto& gridVariables = this->assembler().gridVariables();
- priVarSwitch_->updateBoundary(problem, gridGeometry, gridVariables, sol);
+ priVarSwitch_->updateDirichletConstraints(problem, gridGeometry, gridVariables, sol);
}
/*!
diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
index 49cca272818914e00f886721773ec0ec4f3f9464..d6a02cda41f21fe6bf4d2946c7c3d34c7cb3ea43 100644
--- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
@@ -203,15 +203,30 @@ public:
/*!
* \brief Updates the the primary variables state at the boundary.
*
- * Required when a Dirichlet BC differes from the initial conditon (only for box method).
+ * Required when a Dirichlet BC differs from the initial condition (only for box method).
*/
template<class Problem, class GridVariables, class SolutionVector>
+ [[deprecated("Use updateDirichletConstraints() instead. Will be remove after 3.3")]]
void updateBoundary(const Problem& problem,
const typename GridVariables::GridGeometry& gridGeometry,
GridVariables& gridVariables,
SolutionVector& sol)
{
- if constexpr (GridVariables::GridGeometry::discMethod == DiscretizationMethod::box)
+ updateDirichletConstraints(problem, gridGeometry, gridVariables, sol);
+ }
+
+ /*!
+ * \brief Updates the the primary variables state at the boundary.
+ *
+ * Required when a Dirichlet constraint (at a boundary or internal) differs from the initial condition.
+ */
+ template<class Problem, class GridVariables, class SolutionVector>
+ void updateDirichletConstraints(const Problem& problem,
+ const typename GridVariables::GridGeometry& gridGeometry,
+ GridVariables& gridVariables,
+ SolutionVector& sol)
+ {
+ if constexpr (GridVariables::GridGeometry::discMethod == DiscretizationMethod::box || Problem::enableInternalDirichletConstraints())
{
std::vector<bool> stateChanged(sol.size(), false);
std::size_t countChanged = 0;
@@ -221,8 +236,8 @@ public:
auto fvGeometry = localView(gridGeometry);
fvGeometry.bindElement(element);
- // skip if the element is not at a boundary
- if (!fvGeometry.hasBoundaryScvf())
+ // skip if the element is not at a boundary or if no internal Dirichlet constraints are set
+ if (!Problem::enableInternalDirichletConstraints() && !fvGeometry.hasBoundaryScvf())
continue;
auto elemVolVars = localView(gridVariables.curGridVolVars());
@@ -230,40 +245,32 @@ public:
for (const auto& scv : scvs(fvGeometry))
{
- // this implies that state is set equal for all scvs associated with the dof
const auto dofIdx = scv.dofIndex();
- if (!gridGeometry.dofOnBoundary(dofIdx) || stateChanged[dofIdx])
+
+ if (stateChanged[dofIdx])
continue;
- const auto bcTypes = problem.boundaryTypes(element, scv);
- if (bcTypes.hasDirichlet())
+ if constexpr (GridVariables::GridGeometry::discMethod == DiscretizationMethod::box)
{
- const auto dirichletValues = problem.dirichlet(element, scv);
-
- if (sol[dofIdx].state() != dirichletValues.state())
+ if (gridGeometry.dofOnBoundary(dofIdx))
{
- if (verbosity() > 1)
- std::cout << "Changing primary variable state at boundary (" << sol[dofIdx].state()
- << ") to the one given by the Dirichlet condition (" << dirichletValues.state() << ") at dof " << dofIdx
- << ", coordinates: " << scv.dofPosition()
- << std::endl;
-
- // make sure the solution vector has the right state (given by the Dirichlet BC)
- sol[dofIdx].setState(dirichletValues.state());
- stateChanged[dofIdx] = true;
- ++countChanged;
-
- // overwrite initial with Dirichlet values
- for (int eqIdx = 0; eqIdx < SolutionVector::block_type::dimension; ++eqIdx)
+ if (handleDirichletBoundaryCondition_(problem, element, scv, sol))
{
- if (bcTypes.isDirichlet(eqIdx))
- {
- const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- sol[dofIdx][pvIdx] = dirichletValues[pvIdx];
- }
+ stateChanged[dofIdx] = true;
+ ++countChanged;
+ continue;
}
}
}
+
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ if (handleInternalDirichletConstraint_(problem, element, scv, sol))
+ {
+ stateChanged[dofIdx] = true;
+ ++countChanged;
+ }
+ }
}
// update the volVars if caching is enabled
@@ -285,7 +292,7 @@ public:
}
if (verbosity_ > 0 && countChanged > 0)
- std::cout << "Changed primary variable states and solution values at boundary to Dirichlet states and values at " << countChanged << " dof locations on processor "
+ std::cout << "Changed primary variable states and solution values to Dirichlet states and values at " << countChanged << " dof locations on processor "
<< gridGeometry.gridView().comm().rank() << "." << std::endl;
}
}
@@ -338,10 +345,27 @@ protected:
const typename Geometry::SubControlVolume& scv,
const Problem& problem)
{
- if constexpr (Geometry::GridGeometry::discMethod != DiscretizationMethod::box)
+ // Dofs can be only constrained when using the Box method or when imposing internal Dirichlet constraints
+ if constexpr (Geometry::GridGeometry::discMethod != DiscretizationMethod::box && !Problem::enableInternalDirichletConstraints())
return false;
- else
+ // check for internally constrained Dofs
+ const auto isInternallyConstrainedDof = [&]()
+ {
+ if constexpr (!Problem::enableInternalDirichletConstraints())
+ return false;
+ else
+ {
+ const auto internalDirichletConstraints = problem.hasInternalDirichletConstraint(element, scv);
+ return internalDirichletConstraints.any();
+ }
+ };
+
+ if (isInternallyConstrainedDof())
+ return true;
+
+ // check for a Dirichlet BC when using the Box method
+ if constexpr (Geometry::GridGeometry::discMethod == DiscretizationMethod::box)
{
if (!fvGeometry.hasBoundaryScvf())
return false;
@@ -353,9 +377,87 @@ protected:
const auto bcTypes = problem.boundaryTypes(element, scv);
if (bcTypes.hasDirichlet())
return true;
+ }
- return false;
+ return false;
+ }
+
+ template<class Problem, class Element, class SubControlVolume, class SolutionVector>
+ bool handleDirichletBoundaryCondition_(const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ SolutionVector& sol)
+ {
+ bool changed = false;
+ const auto bcTypes = problem.boundaryTypes(element, scv);
+ if (bcTypes.hasDirichlet())
+ {
+ const auto dirichletValues = problem.dirichlet(element, scv);
+ const auto dofIdx = scv.dofIndex();
+
+ if (sol[dofIdx].state() != dirichletValues.state())
+ {
+ if (verbosity() > 1)
+ std::cout << "Changing primary variable state at boundary (" << sol[dofIdx].state()
+ << ") to the one given by the Dirichlet condition (" << dirichletValues.state() << ") at dof " << dofIdx
+ << ", coordinates: " << scv.dofPosition()
+ << std::endl;
+
+ // make sure the solution vector has the right state (given by the Dirichlet BC)
+ sol[dofIdx].setState(dirichletValues.state());
+ changed = true;
+
+ // overwrite initial with Dirichlet values
+ for (int eqIdx = 0; eqIdx < SolutionVector::block_type::dimension; ++eqIdx)
+ {
+ if (bcTypes.isDirichlet(eqIdx))
+ {
+ const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+ sol[dofIdx][pvIdx] = dirichletValues[pvIdx];
+ }
+ }
+ }
}
+
+ return changed;
+ }
+
+ template<class Problem, class Element, class SubControlVolume, class SolutionVector>
+ bool handleInternalDirichletConstraint_(const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ SolutionVector& sol)
+ {
+ bool changed = false;
+
+ const auto internalDirichletConstraints = problem.hasInternalDirichletConstraint(element, scv);
+ if (internalDirichletConstraints.none())
+ return changed;
+
+ const auto dirichletValues = problem.internalDirichlet(element, scv);
+ const auto dofIdx = scv.dofIndex();
+
+ if (sol[dofIdx].state() != dirichletValues.state())
+ {
+ if (verbosity() > 1)
+ std::cout << "Changing primary variable state at internal DOF (" << sol[dofIdx].state()
+ << ") to the one given by the internal Dirichlet constraint (" << dirichletValues.state() << ") at dof " << dofIdx
+ << ", coordinates: " << scv.dofPosition()
+ << std::endl;
+
+ // make sure the solution vector has the right state (given by the Dirichlet constraint)
+ sol[dofIdx].setState(dirichletValues.state());
+ changed = true;
+
+ // overwrite initial with Dirichlet values
+ for (int pvIdx = 0; pvIdx < SolutionVector::block_type::dimension; ++pvIdx)
+ {
+ if (internalDirichletConstraints[pvIdx])
+ sol[dofIdx][pvIdx] = dirichletValues[pvIdx];
+ }
+ }
+
+ return changed;
}
std::vector<bool> wasSwitched_;
diff --git a/test/porousmediumflow/1p/implicit/internaldirichlet/CMakeLists.txt b/test/porousmediumflow/1p/implicit/internaldirichlet/CMakeLists.txt
index 1b5956b7acf5a1690686c32ad1eff2eb1df8f23f..a206bda3ea01d99e2651f8af9d287e2f6fb69a27 100644
--- a/test/porousmediumflow/1p/implicit/internaldirichlet/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/internaldirichlet/CMakeLists.txt
@@ -1,20 +1,14 @@
-# This feature (internal Dirichlet contraints) is currently not implemented for cc models
-# The compilation triggers a copmiler error (static_assert)
-# Once its implemented replace OnePInternalDirichletBox by OnePInternalDirichletTpfa (and uncomment OnePInternalDirichletTpfa in the problem.hh)
-# and remove the EXPECT_COMPILE_FAIL flag
-
-# dumux_add_test(NAME test_1p_internaldirichlet_tpfa
-# EXPECT_COMPILE_FAIL
-# SOURCES ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/main.cc
-# LABELS porousmediumflow 1p
-# COMPILE_DEFINITIONS TYPETAG=OnePInternalDirichletTpfa
-# COMPILE_DEFINITIONS NUMDIFFMETHOD=DiffMethod::analytic
-# COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-# CMD_ARGS --script fuzzy
-# --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_internaldirichlet_tpfa-reference.vtu
-# ${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_tpfa-00001.vtu
-# --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_tpfa ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/params.input -Problem.Name test_1p_internaldirichlet_tpfa -Problem.EnableGravity false")
+dumux_add_test(NAME test_1p_internaldirichlet_tpfa
+ SOURCES ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/main.cc
+ LABELS porousmediumflow 1p
+ COMPILE_DEFINITIONS TYPETAG=OnePInternalDirichletTpfa
+ COMPILE_DEFINITIONS NUMDIFFMETHOD=DiffMethod::analytic
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_internaldirichlet_tpfa-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_tpfa-00001.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_tpfa ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/params.input -Problem.Name test_1p_internaldirichlet_tpfa -Problem.EnableGravity false")
dumux_add_test(NAME test_1p_internaldirichlet_box
SOURCES ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/main.cc
diff --git a/test/porousmediumflow/1p/implicit/internaldirichlet/problem.hh b/test/porousmediumflow/1p/implicit/internaldirichlet/problem.hh
index e87a9c4f1db916fa2846013cecb9a080c90d7d56..864c32b0e0c7db3c875d59e1840a72be71800066 100644
--- a/test/porousmediumflow/1p/implicit/internaldirichlet/problem.hh
+++ b/test/porousmediumflow/1p/implicit/internaldirichlet/problem.hh
@@ -25,6 +25,7 @@
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_INTERNAL_DIRICHLET_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_INTERNAL_DIRICHLET_HH
+#include <bitset>
#include <dumux/common/boundarytypes.hh>
#include <test/porousmediumflow/1p/implicit/incompressible/problem.hh>
@@ -36,8 +37,7 @@ namespace Properties {
// Create new type tags
namespace TTag {
struct OnePInternalDirichlet {};
-// internal Dirichlet BC are currently not implemented for cc-models
-//struct OnePInternalDirichletTpfa { using InheritsFrom = std::tuple<OnePInternalDirichlet, OnePIncompressibleTpfa>; };
+struct OnePInternalDirichletTpfa { using InheritsFrom = std::tuple<OnePInternalDirichlet, OnePIncompressibleTpfa>; };
struct OnePInternalDirichletBox { using InheritsFrom = std::tuple<OnePInternalDirichlet, OnePIncompressibleBox>; };
} // end namespace TTag
@@ -61,11 +61,15 @@ class OnePTestProblemInternalDirichlet : public OnePTestProblem<TypeTag>
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NeumannValues = GetPropType<TypeTag, Properties::NumEqVector>;
- using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using BoundaryTypes = Dumux::BoundaryTypes<ModelTraits::numEq()>;
+ using Indices = typename ModelTraits::Indices;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ static constexpr auto numEq = ModelTraits::numEq();
+
public:
OnePTestProblemInternalDirichlet(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
@@ -118,11 +122,15 @@ public:
* \param element The finite element
* \param scv The sub-control volume
*/
- bool hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
+ std::bitset<numEq> hasInternalDirichletConstraint(const Element& element,
+ const SubControlVolume& scv) const
{
// the pure Neumann problem is only defined up to a constant
// we create a well-posed problem by fixing the pressure at one dof in the middle of the domain
- return (scv.dofIndex() == static_cast<std::size_t>(this->gridGeometry().numDofs()/2));
+ std::bitset<numEq> values;
+ if (scv.dofIndex() == static_cast<std::size_t>(this->gridGeometry().numDofs()/2))
+ values.set(Indices::pressureIdx);
+ return values;
}
/*!