From f795fe77a1f03425c9588da63ff3d4ea53fc1dbc Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Mon, 20 May 2019 16:38:41 +0200
Subject: [PATCH] [discretization] implement jacobian pattern for fem schemes

---
 dumux/assembly/jacobianpattern.hh | 45 +++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/dumux/assembly/jacobianpattern.hh b/dumux/assembly/jacobianpattern.hh
index 3829a66b73..d3be1809c7 100644
--- a/dumux/assembly/jacobianpattern.hh
+++ b/dumux/assembly/jacobianpattern.hh
@@ -159,6 +159,51 @@ auto getJacobianPattern(const GridGeometry& gridGeometry)
     return pattern;
 }
 
+/*!
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for finite element scheme
+ */
+template<class FEBasis>
+Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis& feBasis)
+{
+    const auto numDofs = feBasis.size();
+
+    Dune::MatrixIndexSet pattern;
+    pattern.resize(numDofs, numDofs);
+
+    // matrix pattern for implicit Jacobians
+    for (const auto& element : elements(feBasis.gridView()))
+    {
+        auto localView = feBasis.localView();
+        localView.bind(element);
+
+        const auto& finiteElement = localView.tree().finiteElement();
+        const auto numLocalDofs = finiteElement.localBasis().size();
+        for (size_t i = 0; i < numLocalDofs; i++)
+        {
+            const auto dofIdxI = localView.index(i);
+            for (size_t j = 0; j < numLocalDofs; j++)
+            {
+                const auto dofIdxJ = localView.index(j);
+                pattern.add(dofIdxI, dofIdxJ);
+            }
+        }
+    }
+
+    return pattern;
+}
+
+/*!
+ * \ingroup Assembly
+ * \brief Helper function to generate Jacobian pattern for finite element scheme
+ * \note This interface is for compatibility with the other schemes. The pattern
+ *       in fem is the same independent of the time discretization scheme.
+ */
+template<bool isImplicit, class GridGeometry,
+         typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethod::fem), int> = 0>
+Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
+{ return getFEJacobianPattern(gridGeometry.feBasis()); }
+
 } // namespace Dumux
 
 #endif
-- 
GitLab