From f2e2f221a9b0351be25927a0df9c188de24e00f7 Mon Sep 17 00:00:00 2001
From: Martin Schneider <martin.schneider@iws.uni-stuttgart.de>
Date: Fri, 21 May 2021 11:12:08 +0200
Subject: [PATCH] [md][box] Internal Dirichlet

* Set off-diagonal matrix blocks to zero if internal Dirichlet is set
* use correct problem() function in md framework
---
 dumux/assembly/fvlocalassemblerbase.hh          |  4 ++--
 dumux/multidomain/subdomainboxlocalassembler.hh | 10 ++++++++++
 2 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/dumux/assembly/fvlocalassemblerbase.hh b/dumux/assembly/fvlocalassemblerbase.hh
index 5b77bb6106..35773d0f51 100644
--- a/dumux/assembly/fvlocalassemblerbase.hh
+++ b/dumux/assembly/fvlocalassemblerbase.hh
@@ -221,10 +221,10 @@ public:
         // and set the residual to (privar - dirichletvalue)
         for (const auto& scvI : scvs(this->fvGeometry()))
         {
-            const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scvI);
+            const auto internalDirichletConstraints = asImp_().problem().hasInternalDirichletConstraint(this->element(), scvI);
             if (internalDirichletConstraints.any())
             {
-                const auto dirichletValues = this->problem().internalDirichlet(this->element(), scvI);
+                const auto dirichletValues = asImp_().problem().internalDirichlet(this->element(), scvI);
                 // set the Dirichlet conditions in residual and jacobian
                 for (int eqIdx = 0; eqIdx < internalDirichletConstraints.size(); ++eqIdx)
                     if (internalDirichletConstraints[eqIdx])
diff --git a/dumux/multidomain/subdomainboxlocalassembler.hh b/dumux/multidomain/subdomainboxlocalassembler.hh
index f81dfaca33..9ad5a02366 100644
--- a/dumux/multidomain/subdomainboxlocalassembler.hh
+++ b/dumux/multidomain/subdomainboxlocalassembler.hh
@@ -388,6 +388,7 @@ class SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*
     using FVElementGeometry = typename GridGeometry::LocalView;
     using GridView = typename 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 };
@@ -581,6 +582,15 @@ public:
                             A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
                         else
                             A[scv.dofIndex()][globalJ][eqIdx][pvIdx] += partialDerivs[scv.localDofIndex()][eqIdx];
+
+                        // enforce internal Dirichlet constraints
+                        if constexpr (Problem::enableInternalDirichletConstraints())
+                        {
+                            const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+                            if (internalDirichletConstraints[eqIdx])
+                                A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
+                        }
+
                     }
                 }
 
-- 
GitLab