docs

refs

wording

additions from review

wording

some work

Author: henrij22

docs/website/01_framework/dirichletBCs.md

@@ -24,23 +24,25 @@ The interface of `#!cpp Ikarus::DirichletValues` is represented by the following
 Ikarus::DirichletValues dirichletValues2(basis); // (1)!
 void fixBoundaryDOFs(f); // (2)!
 void fixDOFs(f); // (3)!
-void fixIthDOF(i); // (4)!
+void setSingleDOF(i, bool flag); // (4)!
 const auto& basis() const; // (5)!
 bool isConstrained(std::size_t i) const; // (6)!
 auto fixedDOFsize() const; // (7)!
 auto size() const ; // (8)!
+auto reset(); // (9)!
 ```
 
 1. Create class by inserting a global basis, [@sander2020dune] Chapter 10.
 2. Accepts a functor to fix boundary degrees of freedom. `f` is  a functor that will be called with the boolean vector of fixed boundary.
  degrees of freedom and the usual arguments of `Dune::Functions::forEachBoundaryDOF`, as defined on page 388 of the Dune
    [@sander2020dune] book.
 3. A more general version of `fixBoundaryDOFs`. Here, a functor is to be provided that accepts a basis and the corresponding boolean
-4. A function that helps to fix the $i$-th degree of freedom
+4. A function that helps to fix or unfix the $i$-th degree of freedom
    vector considering the Dirichlet degrees of freedom.
 5. Returns the underlying basis.
 6. Indicates whether the degree of freedom `i` is fixed.
 7. Returns the number of fixed degrees of freedom.
 8. Returns the number of all dirichlet degrees of freedom.
+9. Resets the whole container
 
 \bibliography

ikarus/python/dirichletvalues/dirichletvalues.hh

@@ -24,9 +24,8 @@
 // PYBIND11_MAKE_OPAQUE(std::vector<bool>);
 namespace Ikarus::Python {
 
-template <class DirichletValues, typename CppVisitor>
-void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::function& functor,
-                            CppVisitor&& cppFunction) {
+template <class DirichletValues>
+void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::function& functor, auto&& cppFunction) {
   using Basis        = typename DirichletValues::Basis;
   using Intersection = typename Basis::GridView::Intersection;
   using BackendType  = typename DirichletValues::BackendType;
@@ -46,24 +45,24 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu
   size_t numParams                = pybind11::len(result);
 
   if (numParams == 2) {
-    auto& function = functor.template cast<const FixBoundaryDOFsWithGlobalIndexFunction>();
-    auto lambda    = [&](BackendType& vec, const MultiIndex& indexGlobal) {
+    auto function = functor.template cast<const FixBoundaryDOFsWithGlobalIndexFunction>();
+    auto lambda   = [&](BackendType& vec, const MultiIndex& indexGlobal) {
       // we explicitly only allow flat indices
       function(vec.vector(), indexGlobal[0]);
     };
     cppFunction(lambda);
 
   } else if (numParams == 3) {
-    auto& function = functor.template cast<const FixBoundaryDOFsWithLocalViewFunction>();
-    auto lambda    = [&](BackendType& vec, int localIndex, auto& lv) {
+    auto function = functor.template cast<const FixBoundaryDOFsWithLocalViewFunction>();
+    auto lambda   = [&](BackendType& vec, int localIndex, auto& lv) {
       auto lvWrapper = LocalViewWrapper(lv.rootLocalView());
       function(vec.vector(), localIndex, lvWrapper);
     };
     cppFunction(lambda);
 
   } else if (numParams == 4) {
-    auto& function = functor.template cast<const FixBoundaryDOFsWithIntersectionFunction>();
-    auto lambda    = [&](BackendType& vec, int localIndex, auto& lv, const Intersection& intersection) {
+    auto function = functor.template cast<const FixBoundaryDOFsWithIntersectionFunction>();
+    auto lambda   = [&](BackendType& vec, int localIndex, auto& lv, const Intersection& intersection) {
       auto lvWrapper = LocalViewWrapper(lv.rootLocalView());
       function(vec.vector(), localIndex, lvWrapper, intersection);
     };
@@ -80,15 +79,21 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu
  * This function registers Python bindings for a DirichletValues class, allowing it to be used in Python scripts.
  * The registered class will have an initializer that takes a `Basis` object. It exposes several member functions to
  * Python:
- *   - `fixBoundaryDOFs(f)`: Fixes boundary degrees of freedom. This function can be overloaded with the following
- * arguments:
- *    - using a user-defined function `f`.
- *    -  using a user-defined function `f` with a `LocalView` argument.
- *    -  using a user-defined function `f` with `LocalView` and `Intersection` arguments.
- *    - `fixDOFs(f)`: Fixes boundary degrees of freedom using a user-defined function `f` with the boolean vector and
- * the basis as arguments.
- *  -  - `fixBoundaryDOFsOfSubSpaceBasis(f)`: Same implementation but you can pass a index of a child basis of a Power
- * Basis
+ *  - `fixBoundaryDOFs(f)`: Fixes boundary degrees of freedom using a user-defined function `f` than can be called
+ * with the following arguments:
+ *   -  with the boolean vector and the global index.
+ *   -  with the boolean vector, the local index and the `LocalView`.
+ *   -  with the boolean vector, the local index, the `LocalView` and the `Intersection`.
+ *  - `fixDOFs(f)`: Fixes boundary degrees of freedom using a user-defined function `f` with the basis and the boolean
+ * vector as arguments.
+ *  - `setSingleDOF(i, flag: bool): Fixes or unfixes DOF with index i
+ *  - `isConstrained(i)`: Checks whether index i is constrained
+ *  - `reset()`: Resets the whole container
+ *
+ * The following properties can be accessed:
+ *  - `container`: the underlying container of dirichlet flags (as a const reference)
+ *  - `size`: the size of the underlying basis
+ *  - `fixedDOFsize`: the amount of DOFs currently fixed
  *
  * \tparam DirichletValues The DirichletValues class to be registered.
  * \tparam options Variadic template parameters for additional options when defining the Python class.
@@ -120,12 +125,12 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_<DirichletV
 
   cls.def(pybind11::init([](const Basis& basis) { return new DirichletValues(basis); }), pybind11::keep_alive<1, 2>());
 
-  cls.def_property_readonly("container", &DirichletValues::container);
+  cls.def_property_readonly("container", &DirichletValues::container, pybind11::return_value_policy::reference);
   cls.def_property_readonly("size", &DirichletValues::size);
   cls.def_property_readonly("fixedDOFsize", &DirichletValues::fixedDOFsize);
   cls.def("isConstrained", [](DirichletValues& self, std::size_t i) -> bool { return self.isConstrained(i); });
-
-  cls.def("fixIthDOF", [](DirichletValues& self, std::size_t i) { self.fixIthDOF(MultiIndex(std::array{i})); });
+  cls.def("setSingleDOF", [](DirichletValues& self, std::size_t i, bool flag) { self.setSingleDOF(i, flag); });
+  cls.def("reset", &DirichletValues::reset);
 
   cls.def("fixDOFs",
           [](DirichletValues& self, const std::function<void(const Basis&, Eigen::Ref<Eigen::VectorX<bool>>)>& f) {

ikarus/python/test/testdirichletvalues.py

@@ -3,7 +3,7 @@
 
 import debug_info
 
-debug_info.unsetDebugFlags()
+debug_info.setDebugFlags()
 
 import ikarus as iks
 import dune.grid
@@ -42,7 +42,7 @@ def fixOneIndex(vec, globalIndex):
     # This is equivalent to values[0] == True, but this syntax is discouraged by PEP
     assert dirichletValues.container[0]
 
-    # Note that the result of localView.index(localIndex) is a multiIndex even for a flat basis, the localIndex appears to be a int
+    # Note that the result of localView.index(localIndex) is a multiIndex even for a flat basis, the localIndex is an int
     def fixAnotherIndexWithLocalView(vec, localIndex, localView):
         if localView.index(localIndex) == [4]:
             vec[localView.index(localIndex)] = True
@@ -52,7 +52,9 @@ def fixAnotherIndexWithLocalView(vec, localIndex, localView):
         assert isinstance(localView.index(localIndex)[0], int)
 
     dirichletValues.fixBoundaryDOFs(fixAnotherIndexWithLocalView)
-    assert sum(dirichletValues.container) == 2
+    container = dirichletValues.container
+
+    assert sum(container) == 2
     assert dirichletValues.fixedDOFsize == 2
 
     def fixTopSide(vec, localIndex, localView, intersection):
@@ -69,6 +71,9 @@ def fixTopSide(vec, localIndex, localView, intersection):
     indicesPerDirection: int = (math.sqrt(grid.size(0)) + 1) * 2
     assert dirichletValues.fixedDOFsize == 2 + indicesPerDirection
 
+    # This checks whether container is a reference
+    assert sum(container) == dirichletValues.fixedDOFsize
+
     # Test Subbasis
     dirichletValues2 = iks.dirichletValues(basis)
 
@@ -81,11 +86,18 @@ def fixTopSide(vec, localIndex, localView, intersection):
     dirichletValues2.fixBoundaryDOFs(fixTopSide, 1)
     assert dirichletValues2.fixedDOFsize == 2 + indicesPerDirection
 
-    dirichletValues2.fixIthDOF(1)
+    dirichletValues2.setSingleDOF(1, True)
     assert dirichletValues2.fixedDOFsize == 2 + indicesPerDirection + 1
     assert dirichletValues2.container[1]
     assert dirichletValues2.isConstrained(1)
-    
+
+    dirichletValues2.setSingleDOF(1, False)
+    assert dirichletValues2.fixedDOFsize == 2 + indicesPerDirection
+    assert not dirichletValues2.isConstrained(1)
+
+    dirichletValues2.reset()
+    assert dirichletValues2.fixedDOFsize == 0
+    assert sum(dirichletValues2.container) == 0
 
 
 if __name__ == "__main__":

ikarus/utils/dirichletvalues.hh

@@ -104,15 +104,21 @@ public:
    * \brief Fixes (set boolean value to true) a specific degree of freedom
    *
    * \param i An index indicating the DOF number to be fixed
+   * \param flag Boolean indicating whether the DOF should fixed or not
    */
-  void fixIthDOF(typename Basis::MultiIndex i) { dirichletFlagsBackend_[i] = true; }
+  template <typename MultiIndex>
+  requires(not std::integral<MultiIndex>)
+  void setSingleDOF(const MultiIndex i, bool flag) {
+    dirichletFlagsBackend_[i] = flag;
+  }
 
   /**
-   * \brief Unfixes (set boolean value to false) a specific degree of freedom
+   * \brief Fixes (set boolean value to true) a specific degree of freedom
    *
    * \param i An index indicating the DOF number to be fixed
+   * \param flag Boolean indicating whether the DOF should fixed or not
    */
-  void unfixIthDOF(typename Basis::MultiIndex i) { dirichletFlagsBackend_[i] = false; }
+  void setSingleDOF(std::size_t i, bool flag) { dirichletFlags_[i] = flag; }
 
   /**
    * \brief Resets all degrees of freedom

tests/src/testdirichletvalue.cpp

@@ -2,8 +2,6 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
 #include <config.h>
 
-#include "testhelpers.hh"
-
 #include <vector>
 
 #include <dune/common/float_cmp.hh>
@@ -44,22 +42,23 @@ static auto dirichletBCTest() {
   using namespace Dune::Functions::BasisFactory;
   auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>()));
 
-  auto basisP = std::make_shared<const decltype(basis)>(basis);
+  auto basisP    = std::make_shared<const decltype(basis)>(basis);
+  auto flatBasis = basisP->flat();
 
-  Ikarus::DirichletValues dirichletValues1(basisP->flat());
+  Ikarus::DirichletValues dirichletValues1(flatBasis);
   dirichletValues1.fixDOFs([](auto& basis_, auto& dirichFlags) {
     Dune::Functions::forEachBoundaryDOF(basis_, [&](auto&& indexGlobal) { dirichFlags[indexGlobal] = true; });
   });
 
-  Ikarus::DirichletValues dirichletValues2(basisP->flat());
+  Ikarus::DirichletValues dirichletValues2(flatBasis);
   dirichletValues2.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; });
 
-  for (std::size_t i = 0; i < basisP->flat().size(); ++i)
+  for (std::size_t i = 0; i < flatBasis.size(); ++i)
     t.check(dirichletValues1.isConstrained(i) == dirichletValues2.isConstrained(i))
         << "Different dirichlet value creations didn't provide the same result. Index: i=" << i;
 
   auto testSubSpaceBasis = [&](auto&& tree1, auto&& tree2) {
-    Ikarus::DirichletValues dirichletValues_SSB(basisP->flat());
+    Ikarus::DirichletValues dirichletValues_SSB(flatBasis);
     dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; },
                                         tree1);
 
@@ -94,10 +93,10 @@ static auto dirichletBCTest() {
     if (intersection.geometry().center()[0] > 4 - 1e-8)
       dirichletFlags[localView.index(localIndex)] = true;
   };
-  Ikarus::DirichletValues dirichletValues5(basisP->flat());
+  Ikarus::DirichletValues dirichletValues5(flatBasis);
   dirichletValues5.fixBoundaryDOFs(fixLambda);
 
-  Ikarus::DirichletValues dirichletValues_SSB2(basisP->flat());
+  Ikarus::DirichletValues dirichletValues_SSB2(flatBasis);
   dirichletValues_SSB2.fixBoundaryDOFs(fixLambda, treePath0);
 
   t.check(dirichletValues5.fixedDOFsize() == dirichletValues_SSB2.fixedDOFsize() * 2)
@@ -113,15 +112,13 @@ static auto dirichletBCTest() {
       << dirichletValues_SSB2.fixedDOFsize() << ", where as the full basis has " << dirichletValues5.fixedDOFsize()
       << " fixed DOFs";
 
-  for (std::size_t i = 0; i < basisP->flat().size(); ++i)
+  for (std::size_t i = 0; i < flatBasis.size(); ++i)
     t.check(dirichletValues_SSB2.isConstrained(i) == dirichletValues5.isConstrained(i))
         << "Different dirichlet value creations with subspace basis didn't provide the same result as with full basis. "
            "Index: i="
         << i;
 
-  auto sum = [](const auto& container_) {
-    return std::accumulate(container_.begin(), container_.end(), 0, std::plus());
-  };
+  auto sum        = [](const auto& container_) { return std::accumulate(container_.begin(), container_.end(), 0); };
   auto manual_sum = [](const auto& dv) {
     int sum = 0;
     for (auto i : Dune::range(dv.size()))
@@ -130,6 +127,18 @@ static auto dirichletBCTest() {
     return sum;
   };
 
+  // Fix non-boundaryDOFs
+  // Ikarus::DirichletValues dirichletValues6(flatBasis);
+  // dirichletValues6.fixDOFs([](const auto& basis, auto& dirichletFlags) {
+  //   Dune::Functions::forEachBoundaryDOF(basis, [&](auto&& localIndex, auto&& localView, auto&& intersection) {
+  //     dirichletFlags[localView.index(localIndex)] = true;
+  //   });
+  // });
+
+  // t.check(dirichletValues6.fixedDOFsize() == flatBasis.size())
+  //     << "All Dofs should be fixed but only " << dirichletValues6.fixedDOFsize() << "out of " << flatBasis.size()
+  //     << "are fixed";
+
   // Test container
   auto& container = dirichletValues2.container();
   static_assert(std::is_reference_v<decltype(container)>, "container should be a reference");
@@ -143,17 +152,21 @@ static auto dirichletBCTest() {
   t.check(sumManual == sumPre)
       << "Summing over container should yield the same amount of fixed DOFs then checking manually";
 
-  dirichletValues2.unfixIthDOF({1});
+  dirichletValues2.setSingleDOF(decltype(flatBasis)::MultiIndex{1}, false);
   t.check(sum(container) == sumPre - 1) << "The sum of fixed DOFs should be one less then after unfixing one";
 
   t.check(manual_sum(dirichletValues2) == sumPre - 1)
       << "Summing over container should yield the same amount of fixed DOFs then checking manually";
 
+  dirichletValues2.setSingleDOF(2, false);
+  t.check(sum(container) == sumPre - 2) << "The sum of fixed DOFs should be two less then after unfixing two";
+
   // Reset
   dirichletValues2.reset();
   t.check(sum(container) == 0) << "After resetting container should have only false entries";
-  t.check(manual_sum(dirichletValues2) == 0) << "After resetting all DOFs should be false";
+  t.check(manual_sum(dirichletValues2) == 0) << "After resetting, all DOFs should be false";
 
+  // Inhomogenious Boundary Conditions
   auto inhomogeneousDisplacement = []<typename T>(const auto& globalCoord, const T& lambda) {
     Eigen::Vector<T, 2> localInhomogeneous;
     if (globalCoord[0] > 4 - 1e-8) {
@@ -251,7 +264,7 @@ static auto dirichletBCTest() {
           (std::abs(nodalPos[i][1]) < tol) or (std::abs(nodalPos[i][1] - Ly) < tol))
         for (auto fixedDirection = 0; fixedDirection < 2; ++fixedDirection) {
           auto fixIndex = localView.index(localView.tree().child(fixedDirection).localIndex(i));
-          dirichletValues4.fixIthDOF(fixIndex);
+          dirichletValues4.setSingleDOF(fixIndex, true);
         }
   }