Fix review comments: add MultiGrid module with correctness and hygiene fixes

Bembel/MultiGrid

@@ -0,0 +1,23 @@
+// This file is part of Bembel, the higher order C++ boundary element library.
+// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
+// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
+// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
+// source code is subject to the GNU General Public License version 3 and
+// provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
+// information.
+
+#ifndef BEMBEL_MULTIGRID_MODULE_
+#define BEMBEL_MULTIGRID_MODULE_
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <unsupported/Eigen/KroneckerProduct>
+
+#include "AnsatzSpace"
+#include "Spline"
+
+#include "src/MultiGrid/ProlongationBsplines.hpp"
+#include "src/MultiGrid/Prolongation.hpp"
+#include "src/MultiGrid/MultiGridSolver.hpp"
+
+#endif  // BEMBEL_MULTIGRID_MODULE_

Bembel/src/MultiGrid/MultiGridSolver.hpp

@@ -0,0 +1,75 @@
+// This file is part of Bembel, the higher order C++ boundary element library.
+//
+// Copyright (C) 2022 see <http://www.bembel.eu>
+//
+// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
+// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
+// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
+// source code is subject to the GNU General Public License version 3 and
+// provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
+// information.
+
+#ifndef BEMBEL_SRC_MULTIGRID_MULTIGRIDSOLVER_HPP_
+#define BEMBEL_SRC_MULTIGRID_MULTIGRIDSOLVER_HPP_
+
+namespace Bembel {
+namespace MG {
+
+enum : unsigned int {
+  MULTIGRID_minLvl = 1,
+  MULTIGRID_cycleType = 2,
+  MULTIGRID_preSmooth = 3,
+  MULTIGRID_postSmooth = 3
+};
+
+/**
+ *  \ingroup MultiGrid
+ *  \brief Add Gauss Seidel smoother for use with multigrid.
+ **/
+template <typename Derived, typename Derived2, typename Derived3>
+void GaussSeidelSmoother(const Derived &S, const Eigen::MatrixBase<Derived2> &b,
+                         Eigen::MatrixBase<Derived3> *x,
+                         const unsigned int steps = 3) {
+  for (auto i = 0; i < steps; ++i)
+    x->derived() = S.template triangularView<Eigen::Lower>().solve(
+        b - S.template triangularView<Eigen::StrictlyUpper>() * x->derived());
+}
+
+/**
+ *  \ingroup MultiGrid
+ *  \brief Implements a single step of multiplicative multigrid.
+ **/
+template <typename Derived, typename Derived2, typename Derived3,
+          typename Derived4>
+void multiplicativeMultiGrid(const Derived &Ss,
+                             const Eigen::MatrixBase<Derived2> &f,
+                             Eigen::MatrixBase<Derived3> *x, const Derived4 &Ps,
+                             const unsigned int lvl = MULTIGRID_minLvl) {
+  typedef typename Derived3::Scalar Scalar;
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
+  if (lvl <= MULTIGRID_minLvl) {
+    Matrix A = Matrix(Ss[lvl]);
+    x->derived() = A.fullPivHouseholderQr().solve(f);
+  } else {
+    // Pre-smooth
+    GaussSeidelSmoother(Ss[lvl], f, x, MULTIGRID_preSmooth);
+    Matrix r = f - Ss[lvl] * x->derived();
+    // Restrict the residual
+    r = Ps[lvl - 1].transpose() * r;
+    Matrix y = r;
+    y.setZero();
+    // Recursively calling multiplicativeMultiGrid for computing the error
+    for (auto i = 0; i < MULTIGRID_cycleType; ++i)
+      multiplicativeMultiGrid(Ss, r, &y, Ps, lvl - 1);
+    // Prolongate the error
+    y = Ps[lvl - 1] * y;
+    // Correction
+    x->derived() += y;
+    // Post-smooth
+    GaussSeidelSmoother(Ss[lvl], f, x, MULTIGRID_postSmooth);
+  }
+}
+
+}  // namespace MG
+}  // namespace Bembel
+#endif  // BEMBEL_SRC_MULTIGRID_MULTIGRIDSOLVER_HPP_

Bembel/src/MultiGrid/Prolongation.hpp

@@ -0,0 +1,65 @@
+// This file is part of Bembel, the higher order C++ boundary element library.
+//
+// Copyright (C) 2022 see <http://www.bembel.eu>
+//
+// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
+// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
+// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
+// source code is subject to the GNU General Public License version 3 and
+// provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
+// information.
+
+#ifndef BEMBEL_SRC_MULTIGRID_PROLONGATION_HPP_
+#define BEMBEL_SRC_MULTIGRID_PROLONGATION_HPP_
+
+namespace Bembel {
+namespace MG {
+
+/**
+ *    \ingroup MultiGrid
+ *    \brief Return the prolongation matrix from level l-1 to level l.
+ **/
+template <typename Derived>
+Eigen::SparseMatrix<double> prolongationMatrix(
+    const AnsatzSpace<Derived> &ansatz_space, unsigned int level) {
+  assert(level > 0 &&
+         "Level must be 1 or larger for build prolongation matrix from level-1 "
+         "to level");
+  unsigned int knot_repetition = ansatz_space.get_knot_repetition();
+  unsigned int polynomial_degree = ansatz_space.get_polynomial_degree();
+  unsigned int num_patches = ansatz_space.get_number_of_patches();
+
+  SuperSpace<Derived> coarse_super_space;
+  coarse_super_space.init_SuperSpace(
+      Geometry(ansatz_space.get_superspace().get_geometry()), level - 1,
+      polynomial_degree);
+  Projector<Derived> coarse_proj(coarse_super_space, knot_repetition);
+  Glue<Derived> coarse_glue(coarse_super_space, coarse_proj);
+
+  SuperSpace<Derived> fine_super_space;
+  fine_super_space.init_SuperSpace(
+      Geometry(ansatz_space.get_superspace().get_geometry()), level,
+      polynomial_degree);
+  Projector<Derived> fine_proj(fine_super_space, knot_repetition);
+  Glue<Derived> fine_glue(fine_super_space, fine_proj);
+
+  const Eigen::SparseMatrix<double> fine_glue_mat =
+      fine_glue.get_glue_matrix();
+  const Eigen::SparseMatrix<double> coarse_glue_mat =
+      coarse_glue.get_glue_matrix();
+
+  Eigen::VectorXd degree_vector =
+      fine_glue_mat.transpose() *
+      Eigen::VectorXd::Ones(fine_glue_mat.rows());
+  // the reason we multiply 0.5 in the begining is because we use the scaled
+  // B-spline basis function in the Bembel
+  return 0.5 * degree_vector.cwiseInverse().asDiagonal() *
+         fine_glue_mat.transpose() *
+         MG::prolongateBsplinesMultiPatches2D(polynomial_degree, level,
+                                              num_patches, knot_repetition) *
+         coarse_glue_mat;
+}
+
+}  // namespace MG
+}  // namespace Bembel
+#endif  // BEMBEL_SRC_MULTIGRID_PROLONGATION_HPP_

Bembel/src/MultiGrid/ProlongationBsplines.hpp

@@ -0,0 +1,144 @@
+// This file is part of Bembel, the higher order C++ boundary element library.
+//
+// Copyright (C) 2022 see <http://www.bembel.eu>
+//
+// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
+// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
+// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
+// source code is subject to the GNU General Public License version 3 and
+// provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
+// information.
+
+#ifndef BEMBEL_SRC_MULTIGRID_PROLONGATIONBSPLINES_HPP_
+#define BEMBEL_SRC_MULTIGRID_PROLONGATIONBSPLINES_HPP_
+
+#include <tuple>
+
+namespace Bembel {
+namespace MG {
+
+template <typename Scalar>
+std::vector<Eigen::Triplet<Scalar>> toEigenTriplets(
+    Eigen::SparseMatrix<Scalar> &M) {
+  std::vector<Eigen::Triplet<Scalar>> v;
+  for (int i = 0; i < M.outerSize(); i++)
+    for (typename Eigen::SparseMatrix<Scalar>::InnerIterator it(M, i); it; ++it)
+      v.emplace_back(it.row(), it.col(), it.value());
+  return v;
+}
+
+/**
+ *  \ingroup MultiGrid
+ *  \brief Return the nonzero entries info of the row id of the prolongation
+ * matrix using Oslo-Algorithm
+ * https://www.uio.no/studier/emner/matnat/ifi/nedlagte-emner/INF-MAT5340/v07/undervisningsmateriale/kap4.pdf
+ */
+std::tuple<Eigen::MatrixXd, unsigned int> nonZerosInfoRow(
+    const std::vector<double> &knots,
+    const std::vector<double> &uniform_refined_knots,
+    unsigned int polynomial_degree, unsigned int row_id,
+    unsigned int guess = 0) {
+  unsigned int d = guess;
+  const unsigned int knots_size = static_cast<unsigned int>(knots.size());
+  while (d + 1 < knots_size &&
+         uniform_refined_knots[row_id] >= knots[d]) ++d;
+  if (d > 0) --d;
+  Eigen::MatrixXd non_zero_alpha =
+      Eigen::MatrixXd::Ones(1, polynomial_degree + 1);
+  for (auto i = 0; i < polynomial_degree; ++i) {
+    Eigen::MatrixXd R = Eigen::MatrixXd::Zero(i + 1, i + 2);
+    for (auto j = 0; j < i + 1; ++j) {
+      R(j, j) = (knots[d + j + 1] - uniform_refined_knots[row_id + i + 1]) /
+                (knots[d + j + 1] - knots[d + j - i]);
+      R(j, j + 1) = 1.0 - R(j, j);
+    }
+    non_zero_alpha.block(0, 0, 1, i + 2) =
+        non_zero_alpha.block(0, 0, 1, i + 1) * R;
+  }
+
+  return std::make_tuple(non_zero_alpha, d);
+}
+
+/**
+ *  \ingroup MultiGrid
+ *  \brief Return prolongation matrix from  level-1 to level in 1D.
+ */
+Eigen::SparseMatrix<double> prolongateBsplines1D(
+    unsigned int polynomial_degree, unsigned int level,
+    unsigned int knotrepetition = 1) {
+  assert(level > 0 &&
+         "Level must be 1 or larger for build prolongation matrix from level-1 "
+         "to level");
+  std::vector<double> knots = Bembel::Spl::MakeUniformKnotVector(
+      polynomial_degree + 1, ((1 << (level - 1)) - 1), knotrepetition);
+  std::vector<double> uniform_refined_knots = Spl::MakeUniformKnotVector(
+      polynomial_degree + 1, ((1 << level) - 1), knotrepetition);
+  unsigned int num_basis = knots.size() - polynomial_degree - 1;
+  unsigned int num_fine_basis =
+      uniform_refined_knots.size() - polynomial_degree - 1;
+  typedef typename Eigen::Triplet<double> Triplet;
+  std::vector<Triplet> tripletList;
+  tripletList.reserve(num_fine_basis * (polynomial_degree + 1));
+  int d = 0;
+  for (int i = 0; i < num_fine_basis; ++i) {
+    Eigen::MatrixXd non_zero_alpha;
+    std::tie(non_zero_alpha, d) = nonZerosInfoRow(
+        knots, uniform_refined_knots, polynomial_degree, i, d);
+    for (int j = 0; j < polynomial_degree + 1; ++j) {
+      tripletList.push_back(
+          Triplet(i, d - polynomial_degree + j, non_zero_alpha(0, j)));
+    }
+  }
+  Eigen::SparseMatrix<double> P;
+  P.resize(num_fine_basis, num_basis);
+  P.setFromTriplets(tripletList.begin(), tripletList.end());
+  return P;
+}
+
+/**
+ *  \ingroup MultiGrid
+ *  \brief Return prolongation matrix from  level-1 to level in 2D.
+ */
+Eigen::SparseMatrix<double> prolongateBsplines2D(
+    unsigned int polynomial_degree, unsigned int level,
+    unsigned int knotrepetition = 1) {
+  assert(level > 0 &&
+         "Level must be 1 or larger for build prolongation matrix from level-1 "
+         "to level");
+  Eigen::SparseMatrix<double> P_ =
+      prolongateBsplines1D(polynomial_degree, level, knotrepetition);
+  Eigen::SparseMatrix<double> P = Eigen::kroneckerProduct(P_, P_);
+  return P;
+}
+
+/**
+ *  \ingroup MultiGrid
+ *  \brief Return prolongation matrix from  level-1 to level in 2D in  case of
+ * multiple patches.
+ */
+Eigen::SparseMatrix<double> prolongateBsplinesMultiPatches2D(
+    unsigned int polynomial_degree, unsigned int level,
+    unsigned int num_patches, unsigned int knotrepetition = 1) {
+  assert(level > 0 &&
+         "Level must be 1 or larger for build prolongation matrix from level-1 "
+         "to level");
+  Eigen::SparseMatrix<double> P_ =
+      prolongateBsplines2D(polynomial_degree, level, knotrepetition);
+  std::vector<Eigen::Triplet<double>> triplets_ = toEigenTriplets(P_);
+  std::vector<Eigen::Triplet<double>> triplets;
+  Eigen::SparseMatrix<double> P(num_patches * P_.rows(),
+                                num_patches * P_.cols());
+  for (auto i = 0; i < num_patches; ++i) {
+    for (auto entry : triplets_) {
+      triplets.push_back(Eigen::Triplet<double>(entry.row() + i * P_.rows(),
+                                                entry.col() + i * P_.cols(),
+                                                entry.value()));
+    }
+  }
+  P.setFromTriplets(triplets.begin(), triplets.end());
+  return P;
+}
+
+}  // namespace MG
+}  // namespace Bembel
+#endif  // BEMBEL_SRC_MULTIGRID_PROLONGATIONBSPLINES_HPP_

examples/CMakeLists.txt

@@ -14,6 +14,7 @@ set(CIFILES
     Geometry
     MassMatrix
     LaplaceBeltrami
+    LaplaceBeltramiMultiGrid
     LaplaceSingleLayerFull
     LaplaceSingleLayerH2
     LaplaceHypersingularH2