[disc] Enable optional internal Dirichlet constraints

dumux/assembly/cclocalassembler.hh

@@ -57,11 +57,13 @@ class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Imp
 {
     using ParentType = FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>;
     using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Problem = typename Assembler::Problem;
     using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
     using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
     using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
     using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
     using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
 
 public:
 
@@ -86,6 +88,20 @@ public:
         {
             res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
         }
+
+        auto applyDirichlet = [&] (const auto& scvI,
+                                   const auto& dirichletValues,
+                                   const auto eqIdx,
+                                   const auto pvIdx)
+        {
+            res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+            auto& row = jac[scvI.dofIndex()];
+            for (auto col = row.begin(); col != row.end(); ++col)
+                row[col.index()][eqIdx] = 0.0;
+            row[scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+        };
+
+        this->asImp_().enforceDirichletConstraints(applyDirichlet);
     }
 
     /*!
@@ -96,6 +112,19 @@ public:
     {
         this->asImp_().bindLocalViews();
         this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
+
+        auto applyDirichlet = [&] (const auto& scvI,
+                                   const auto& dirichletValues,
+                                   const auto eqIdx,
+                                   const auto pvIdx)
+        {
+            auto& row = jac[scvI.dofIndex()];
+            for (auto col = row.begin(); col != row.end(); ++col)
+                row[col.index()][eqIdx] = 0.0;
+            row[scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+        };
+
+        this->asImp_().enforceDirichletConstraints(applyDirichlet);
     }
 
     /*!
@@ -106,7 +135,50 @@ public:
         this->asImp_().bindLocalViews();
         const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(this->element());
         res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
+
+        auto applyDirichlet = [&] (const auto& scvI,
+                                   const auto& dirichletValues,
+                                   const auto eqIdx,
+                                   const auto pvIdx)
+        {
+            res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+        };
+
+        this->asImp_().enforceDirichletConstraints(applyDirichlet);
     }
+
+    //! Enforce Dirichlet constraints
+    template<typename ApplyFunction>
+    void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
+    {
+        // currently Dirichlet boundary conditions are weakly enforced for cc-schemes
+        // so here, we only take care of internal Dirichlet constraints (if enabled)
+        this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
+    }
+
+    /*!
+     * \brief Enforces Dirichlet constraints if enabled in the problem
+     */
+    template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<P::enableInternalDirichletConstraints(), int> = 0>
+    void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
+    {
+        // enforce Dirichlet constraints strongly by overwriting partial derivatives with 1 or 0
+        // and set the residual to (privar - dirichletvalue)
+        for (const auto& scvI : scvs(this->fvGeometry()))
+        {
+            if (this->problem().hasInternalDirichletConstraint(this->element(), scvI))
+            {
+                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);
+            }
+        }
+    }
+
+    template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<!P::enableInternalDirichletConstraints(), int> = 0>
+    void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
+    {}
 };
 
 /*!

dumux/common/fvproblem.hh

@@ -246,6 +246,17 @@ public:
                    "a dirichlet() or a dirichletAtPos() method.");
     }
 
+    /*!
+     * \brief If internal Dirichlet contraints are enabled
+     * Enables / disables internal (non-boundary) Dirichlet constraints. If this is overloaded
+     * to return true, the assembler calls problem.hasInternalDirichletConstraint(element, scv)
+     * which has to return a bool signifying whether the dof associated with is contraint.
+     * If true is returned for a dof, the assembler calls problem.internalDiririchlet(element, scv)
+     * which has to return the enforced Dirichlet values for that dof.
+     */
+    static constexpr bool enableInternalDirichletConstraints()
+    { return false; }
+
     /*!
      * \brief Evaluate the boundary conditions for a neumann
      *        boundary segment.

dumux/multidomain/subdomaincclocalassembler.hh

@@ -82,6 +82,8 @@ class SubDomainCCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assem
     using CouplingManager = typename Assembler::CouplingManager;
     using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
 
+    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
+
 public:
     //! export the domain id of this sub-domain
     static constexpr auto domainId = typename Dune::index_constant<id>();
@@ -127,6 +129,20 @@ public:
             if (i != id)
                 this->assembleJacobianCoupling(i, jacRow, res[globalI], gridVariables);
         });
+
+        auto applyDirichlet = [&] (const auto& scvI,
+                                   const auto& dirichletValues,
+                                   const auto eqIdx,
+                                   const auto pvIdx)
+        {
+            res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+            auto& row = jacRow[domainId][scvI.dofIndex()];
+            for (auto col = row.begin(); col != row.end(); ++col)
+                row[col.index()][eqIdx] = 0.0;
+            row[scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+        };
+
+        this->asImp_().enforceDirichletConstraints(applyDirichlet);
     }
 
     /*!
@@ -158,6 +174,16 @@ public:
         this->asImp_().bindLocalViews();
         const auto globalI = this->fvGeometry().fvGridGeometry().elementMapper().index(this->element());
         res[globalI] = this->evalLocalResidual()[0]; // forward to the internal implementation
+
+        auto applyDirichlet = [&] (const auto& scvI,
+                                   const auto& dirichletValues,
+                                   const auto eqIdx,
+                                   const auto pvIdx)
+        {
+            res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+        };
+
+        this->asImp_().enforceDirichletConstraints(applyDirichlet);
     }
 
     /*!
@@ -243,6 +269,39 @@ public:
     const Problem& problem() const
     { return this->assembler().problem(domainId); }
 
+    //! Enforce Dirichlet constraints
+    template<typename ApplyFunction>
+    void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
+    {
+        // currently Dirichlet boundary conditions are weakly enforced for cc-schemes
+        // so here, we only take care of internal Dirichlet constraints (if enabled)
+        this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
+    }
+
+    /*!
+     * \brief Enforces Dirichlet constraints if enabled in the problem
+     */
+    template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<P::enableInternalDirichletConstraints(), int> = 0>
+    void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
+    {
+        // enforce Dirichlet constraints strongly by overwriting partial derivatives with 1 or 0
+        // and set the residual to (privar - dirichletvalue)
+        for (const auto& scvI : scvs(this->fvGeometry()))
+        {
+            if (this->problem().hasInternalDirichletConstraint(this->element(), scvI))
+            {
+                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);
+            }
+        }
+    }
+
+    template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<!P::enableInternalDirichletConstraints(), int> = 0>
+    void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
+    {}
+
     //! return reference to the coupling manager
     CouplingManager& couplingManager()
     { return couplingManager_; }