[md] Add periodic constraints to multidomain fvassembler

dumux/multidomain/fvassembler.hh

@@ -203,6 +203,9 @@ public:
             auto& jacRow = (*jacobian_)[domainId];
             auto& subRes = (*residual_)[domainId];
             this->assembleJacobianAndResidual_(domainId, jacRow, subRes, curSol);
+
+            const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
+            enforcePeriodicConstraints_(domainId, jacRow, subRes, *gridGeometry, curSol[domainId]);
         });
     }
 
@@ -551,6 +554,48 @@ private:
                                                         domainJ, gridGeometry(domainJ));
     }
 
+    // build periodic contraints into the system matrix
+    template<std::size_t i, class JacRow, class Sol, class GG>
+    void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Sol& res, const GG& gridGeometry, const Sol& curSol)
+    {
+        if constexpr (GG::discMethod == DiscretizationMethod::box || GG::discMethod == DiscretizationMethod::fcstaggered)
+        {
+            for (const auto& m : gridGeometry.periodicVertexMap())
+            {
+                if (m.first < m.second)
+                {
+                    auto& jac = jacRow[domainI];
+
+                    // add the second row to the first
+                    res[m.first] += res[m.second];
+
+                    // enforce the solution of the first periodic DOF to the second one
+                    res[m.second] = curSol[m.second] - curSol[m.first];
+
+                    const auto end = jac[m.second].end();
+                    for (auto it = jac[m.second].begin(); it != end; ++it)
+                        jac[m.first][it.index()] += (*it);
+
+                    // enforce constraint in second row
+                    for (auto it = jac[m.second].begin(); it != end; ++it)
+                        (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
+
+                    using namespace Dune::Hybrid;
+                    forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
+                    {
+                        auto& jacCoupling = jacRow[couplingDomainId];
+
+                        for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
+                            jacCoupling[m.first][it.index()] += (*it);
+
+                        for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
+                            (*it) = 0.0;
+                    });
+                }
+            }
+        }
+    }
+
     //! pointer to the problem to be solved
     ProblemTuple problemTuple_;