From 1dc10fbb7e1548f008b51a3bed71dcacfbd88439 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 24 Jul 2018 15:32:55 +0200
Subject: [PATCH] [assembly][box] treat ghost elements explicitly

This recovers results from 2.12. However, these results also are wrong.
---
 dumux/assembly/boxlocalassembler.hh | 32 ++++++++++++++++++++++++++++-
 1 file changed, 31 insertions(+), 1 deletion(-)

diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index ad9baffa0b..c8ebb3d836 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -88,12 +88,42 @@ public:
             for (const auto& scv : scvs(this->fvGeometry()))
                 res[scv.dofIndex()] += residual[scv.localDofIndex()];
         }
-        else
+        else if (!this->elementIsGhost())
         {
             const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables, partialReassembler); // forward to the internal implementation
             for (const auto& scv : scvs(this->fvGeometry()))
                 res[scv.dofIndex()] += residual[scv.localDofIndex()];
         }
+        else
+        {
+            using GridGeometry = typename GridVariables::GridGeometry;
+            using GridView = typename GridGeometry::GridView;
+            static constexpr auto dim = GridView::dimension;
+
+            int numVerticesLocal = this->element().subEntities(dim);
+
+            for (int i = 0; i < numVerticesLocal; ++i) {
+                auto vertex = this->element().template subEntity<dim>(i);
+
+                if (vertex.partitionType() == Dune::InteriorEntity ||
+                    vertex.partitionType() == Dune::BorderEntity)
+                {
+                    // do not change the non-ghost vertices
+                    continue;
+                }
+
+                // set main diagonal entries for the vertex
+                int vIdx = this->assembler().fvGridGeometry().vertexMapper().index(vertex);
+
+                typedef typename JacobianMatrix::block_type BlockType;
+                BlockType &J = jac[vIdx][vIdx];
+                for (int j = 0; j < BlockType::rows; ++j)
+                    J[j][j] = 1.0;
+
+                // set residual for the vertex
+                res[vIdx] = 0;
+            }
+        }
 
         auto applyDirichlet = [&] (const auto& scvI,
                                    const auto& dirichletValues,
-- 
GitLab