From 46b138c0404ce9a8b1240049b132da38d29d01f8 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 25 Jun 2024 18:33:21 +0200 Subject: [PATCH 01/18] add new test homogenize Python Dirichlet interface + subspace basis fixing Fix material fix3 docs --- .../python/dirichletvalues/dirichletvalues.hh | 87 +++++++++++++------ ikarus/python/finiteelements/fe.hh | 4 +- ikarus/python/test/CMakeLists.txt | 11 +++ ikarus/python/test/kltest.py | 2 +- ikarus/python/test/linearelastictest.py | 36 ++++---- ikarus/python/test/nonlinearelastictest.py | 28 +++--- ikarus/python/test/testdirichletvalues.py | 75 ++++++++++++++++ ikarus/python/test/testmaterials.py | 1 + ikarus/utils/dirichletvalues.hh | 27 ++++++ tests/src/testdirichletvalue.cpp | 22 +++++ 10 files changed, 236 insertions(+), 57 deletions(-) create mode 100644 ikarus/python/test/testdirichletvalues.py diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 2b6320457..5836aff84 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -20,6 +20,7 @@ #include #include + // PYBIND11_MAKE_OPAQUE(std::vector); namespace Ikarus::Python { @@ -57,6 +58,7 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_ LocalViewWrapper; auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; auto lv = Dune::Python::insertClass( @@ -66,32 +68,30 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_()); - // Eigen::Ref needed due to https://pybind11.readthedocs.io/en/stable/advanced/cast/eigen.html#pass-by-reference - cls.def("fixBoundaryDOFs", - [](DirichletValues& self, const std::function>, int)>& f) { - auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { - // we explicitly only allow flat indices - f(vec.vector(), indexGlobal[0]); - }; - self.fixBoundaryDOFs(lambda); - }); + auto fixBoundaryDOFs_ = [](DirichletValues& self, + const std::function>, int)>& f) { + auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { + // we explicitly only allow flat indices + f(vec.vector(), indexGlobal[0]); + }; + self.fixBoundaryDOFs(lambda); + }; - cls.def("fixBoundaryDOFsUsingLocalView", - [](DirichletValues& self, - const std::function>, int, LocalViewWrapper&)>& f) { - auto lambda = [&](BackendType& vec, int localIndex, LocalView& lv) { - auto lvWrapper = LocalViewWrapper(lv.globalBasis()); - // this can be simplified when - // https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/418 becomes available - pybind11::object obj = pybind11::cast(lv.element()); - lvWrapper.bind(obj); - f(vec.vector(), localIndex, lvWrapper); - }; - self.fixBoundaryDOFs(lambda); - }); + auto fixBoundaryDOFsUsingLocalView_ = + [](DirichletValues& self, + const std::function>, int, LocalViewWrapper&)>& f) { + auto lambda = [&](BackendType& vec, int localIndex, LocalView& lv) { + auto lvWrapper = LocalViewWrapper(lv.globalBasis()); + // this can be simplified when + // https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/418 becomes available + pybind11::object obj = pybind11::cast(lv.element()); + lvWrapper.bind(obj); + f(vec.vector(), localIndex, lvWrapper); + }; + self.fixBoundaryDOFs(lambda); + }; - cls.def( - "fixBoundaryDOFsUsingLocalViewAndIntersection", + auto fixBoundaryDOFsUsingLocalViewAndIntersection_ = [](DirichletValues& self, const std::function>, int, LocalViewWrapper&, const Intersection&)>& f) { auto lambda = [&](BackendType& vec, int localIndex, LocalView& lv, const Intersection& intersection) { @@ -103,7 +103,42 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_>, int)>; + using FixBoundaryDOFsWithLocalViewFunction = + std::function>, int, LocalViewWrapper&)>; + using FixBoundaryDOFsWithIntersectionFunction = + std::function>, int, LocalViewWrapper&, const Intersection&)>; + + // Disambiguate by number of arguments, as casting doesn't properly work with functions + pybind11::module inspect_module = pybind11::module::import("inspect"); + pybind11::object result = inspect_module.attr("signature")(functor).attr("parameters"); + size_t numParams = pybind11::len(result); + + if (numParams == 2) { + auto& function = functor.template cast(); + fixBoundaryDOFs_(self, function); + + } else if (numParams == 3) { + auto& function = functor.template cast(); + fixBoundaryDOFsUsingLocalView_(self, function); + + } else if (numParams == 4) { + auto& function = functor.template cast(); + fixBoundaryDOFsUsingLocalViewAndIntersection_(self, function); + + } else { + DUNE_THROW(Dune::NotImplemented, "fixBoundaryDOFs: A function with this signature is not supported"); + } + }, + pybind11::arg("functor")); + + cls.def("fixBoundaryDOFsOfSubSpaceBasis", [](DirichletValues& self, const pybind11::function& functor, + const pybind11::object& ssb) { pybind11::print(ssb); }); cls.def("fixDOFs", [](DirichletValues& self, const std::function>)>& f) { @@ -113,6 +148,8 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_ int { return self.size(); }); } } // namespace Ikarus::Python diff --git a/ikarus/python/finiteelements/fe.hh b/ikarus/python/finiteelements/fe.hh index ec29ad809..6f0e6a136 100644 --- a/ikarus/python/finiteelements/fe.hh +++ b/ikarus/python/finiteelements/fe.hh @@ -48,7 +48,7 @@ namespace Ikarus::Python { * \throws Dune::NotImplemented If the specified resultType is not supported by the finite element. */ template -void registerCalculateAt(pybind11::handle scope, pybind11::class_ cls, auto restultTypesTuple) { +void registerCalculateAt(pybind11::handle scope, pybind11::class_ cls, auto resultTypesTuple) { using Traits = typename FE::Traits; using FERequirements = typename FE::Requirement; cls.def( @@ -57,7 +57,7 @@ void registerCalculateAt(pybind11::handle scope, pybind11::class_(RT i) { + Dune::Hybrid::forEach(resultTypesTuple, [&](RT i) { if (resType == toString(i)) { success = true; result = self.template calculateAt(req, local).asVec(); diff --git a/ikarus/python/test/CMakeLists.txt b/ikarus/python/test/CMakeLists.txt index 1c869f9d1..e49554b13 100644 --- a/ikarus/python/test/CMakeLists.txt +++ b/ikarus/python/test/CMakeLists.txt @@ -34,6 +34,17 @@ dune_python_add_test( python ) +dune_python_add_test( + NAME + pydirichletvalues + SCRIPT + testdirichletvalues.py + WORKING_DIRECTORY + ${CMAKE_CURRENT_SOURCE_DIR} + LABELS + python +) + if(HAVE_DUNE_IGA) dune_python_add_test( NAME diff --git a/ikarus/python/test/kltest.py b/ikarus/python/test/kltest.py index 64abc2c73..e6957c480 100644 --- a/ikarus/python/test/kltest.py +++ b/ikarus/python/test/kltest.py @@ -82,7 +82,7 @@ def fixLeftAndRightEdge(vec, localIndex, localView, intersection): ): vec[localView.index(localIndex)] = True - dirichletValues.fixBoundaryDOFsUsingLocalViewAndIntersection(fixLeftAndRightEdge) + dirichletValues.fixBoundaryDOFs(fixLeftAndRightEdge) assembler = iks.assembler.sparseFlatAssembler(fes, dirichletValues) diff --git a/ikarus/python/test/linearelastictest.py b/ikarus/python/test/linearelastictest.py index 5afb65dd0..70e13e9c7 100644 --- a/ikarus/python/test/linearelastictest.py +++ b/ikarus/python/test/linearelastictest.py @@ -1,8 +1,8 @@ # SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de # SPDX-License-Identifier: LGPL-3.0-or-later - import debug_info + debug_info.setDebugFlags() import ikarus as iks @@ -57,20 +57,26 @@ def loadTopEdgePredicate(x): indexSet = grid.indexSet for v in grid.vertices: - neumannVertices[indexSet.index(v)]=loadTopEdgePredicate(v.geometry.center) + neumannVertices[indexSet.index(v)] = loadTopEdgePredicate(v.geometry.center) boundaryPatch = iks.utils.boundaryPatch(grid, neumannVertices) - nBLoad= iks.finite_elements.neumannBoundaryLoad(boundaryPatch,neumannLoad) + + nBLoad = iks.finite_elements.neumannBoundaryLoad(boundaryPatch, neumannLoad) linElastic = iks.finite_elements.linearElastic(youngs_modulus=1000, nu=0.2) - easF= iks.finite_elements.eas(4) + easF = iks.finite_elements.eas(4) for e in grid.elements: if easBool: - - fes.append(iks.finite_elements.makeFE(basisLagrange1,linElastic,easF,vLoad,nBLoad)) + fes.append( + iks.finite_elements.makeFE( + basisLagrange1, linElastic, easF, vLoad, nBLoad + ) + ) else: - fes.append(iks.finite_elements.makeFE(basisLagrange1,linElastic,vLoad,nBLoad)) + fes.append( + iks.finite_elements.makeFE(basisLagrange1, linElastic, vLoad, nBLoad) + ) fes[-1].bind(e) req = fes[0].createRequirement() @@ -90,8 +96,8 @@ def loadTopEdgePredicate(x): forces = np.zeros(8) stiffness = np.zeros((8, 8)) - fes[0].calculateVector(req,iks.VectorAffordance.forces, forces) - fes[0].calculateMatrix(req,iks.MatrixAffordance.stiffness, stiffness) + fes[0].calculateVector(req, iks.VectorAffordance.forces, forces) + fes[0].calculateMatrix(req, iks.MatrixAffordance.stiffness, stiffness) fes[0].localView() dirichletValues = iks.dirichletValues(flatBasis) @@ -108,8 +114,8 @@ def fixLeftHandEdge(vec, localIndex, localView, intersection): vec[localView.index(localIndex)] = True dirichletValues.fixBoundaryDOFs(fixFirstIndex) - dirichletValues.fixBoundaryDOFsUsingLocalView(fixAnotherVertex) - dirichletValues.fixBoundaryDOFsUsingLocalViewAndIntersection(fixLeftHandEdge) + dirichletValues.fixBoundaryDOFs(fixAnotherVertex) + dirichletValues.fixBoundaryDOFs(fixLeftHandEdge) assembler = iks.assembler.sparseFlatAssembler(fes, dirichletValues) assemblerDense = iks.assembler.denseFlatAssembler(fes, dirichletValues) @@ -134,17 +140,15 @@ def fixLeftHandEdge(vec, localIndex, localView, intersection): ) # Writing results into vtk file - fileName = "resultdisplacement"+ ("EAS" if easBool else "") + fileName = "resultdisplacement" + ("EAS" if easBool else "") - writer = vtkWriter( - grid, fileName, pointData={("displacement", (0, 1)): fx} - ) + vtkWriter(grid, fileName, pointData={("displacement", (0, 1)): fx}) writer2 = vtkUnstructuredGridWriter(grid) writer2.addCellData(stressFuncScalar, name="stress") writer2.addCellData(stressFuncVec, name="stress2") - writer2.write(name= "result"+ ("EAS" if easBool else "")) + writer2.write(name="result" + ("EAS" if easBool else "")) # Querying for a different ResultType should result in a runtime error try: diff --git a/ikarus/python/test/nonlinearelastictest.py b/ikarus/python/test/nonlinearelastictest.py index cc30e5075..8d8f44d7c 100644 --- a/ikarus/python/test/nonlinearelastictest.py +++ b/ikarus/python/test/nonlinearelastictest.py @@ -2,6 +2,7 @@ # SPDX-License-Identifier: LGPL-3.0-or-later import debug_info + debug_info.setDebugFlags() import ikarus as iks @@ -47,11 +48,11 @@ def loadTopEdgePredicate(x): indexSet = grid.indexSet for v in grid.vertices: - neumannVertices[indexSet.index(v)]=loadTopEdgePredicate(v.geometry.center) + neumannVertices[indexSet.index(v)] = loadTopEdgePredicate(v.geometry.center) boundaryPatch = iks.utils.boundaryPatch(grid, neumannVertices) - nBLoad= iks.finite_elements.neumannBoundaryLoad(boundaryPatch,neumannLoad) + nBLoad = iks.finite_elements.neumannBoundaryLoad(boundaryPatch, neumannLoad) svk = iks.materials.StVenantKirchhoff(E=1000, nu=0.3) @@ -61,7 +62,9 @@ def loadTopEdgePredicate(x): fes = [] for e in grid.elements: - fes.append(iks.finite_elements.makeFE(basisLagrange1,nonLinElastic,vLoad,nBLoad)) + fes.append( + iks.finite_elements.makeFE(basisLagrange1, nonLinElastic, vLoad, nBLoad) + ) fes[-1].bind(e) dirichletValues = iks.dirichletValues(flatBasis) @@ -78,8 +81,8 @@ def fixLeftHandEdge(vec, localIndex, localView, intersection): vec[localView.index(localIndex)] = True dirichletValues.fixBoundaryDOFs(fixFirstIndex) - dirichletValues.fixBoundaryDOFsUsingLocalView(fixAnotherVertex) - dirichletValues.fixBoundaryDOFsUsingLocalViewAndIntersection(fixLeftHandEdge) + dirichletValues.fixBoundaryDOFs(fixAnotherVertex) + dirichletValues.fixBoundaryDOFs(fixLeftHandEdge) assembler = iks.assembler.sparseFlatAssembler(fes, dirichletValues) @@ -91,13 +94,11 @@ def fixLeftHandEdge(vec, localIndex, localView, intersection): def energy(dRedInput): feReq = fes[0].createRequirement() - feReq.insertParameter( lambdaLoad) + feReq.insertParameter(lambdaLoad) dBig = assembler.createFullVector(dRedInput) - feReq.insertGlobalSolution( dBig) + feReq.insertGlobalSolution(dBig) feReq.globalSolution() - return assembler.scalar( - feReq, iks.ScalarAffordance.mechanicalPotentialEnergy - ) + return assembler.scalar(feReq, iks.ScalarAffordance.mechanicalPotentialEnergy) def gradient(dRedInput): feReq = fes[0].createRequirement() @@ -105,7 +106,8 @@ def gradient(dRedInput): dBig = assembler.createFullVector(dRedInput) feReq.insertGlobalSolution(dBig) return assembler.vector( - feReq, iks.VectorAffordance.forces, iks.DBCOption.Reduced) + feReq, iks.VectorAffordance.forces, iks.DBCOption.Reduced + ) def hess(dRedInput): feReq = fes[0].createRequirement() @@ -114,7 +116,7 @@ def hess(dRedInput): feReq.insertGlobalSolution(dBig) return assembler.matrix( feReq, iks.MatrixAffordance.stiffness, iks.DBCOption.Reduced - ).todense() # this is slow, but for this test we don't care + ).todense() # this is slow, but for this test we don't care resultd = minimize(energy, x0=dRed, options={"disp": True}, tol=1e-14) resultd2 = minimize( @@ -137,6 +139,6 @@ def hess(dRedInput): feReq = fes[0].createRequirement() fullD = assembler.createFullVector(resultd2.x) - feReq.insertGlobalSolution( fullD) + feReq.insertGlobalSolution(fullD) res1 = fes[0].calculateAt(feReq, np.array([0.5, 0.5]), "PK2Stress") diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py new file mode 100644 index 000000000..8d2168d5d --- /dev/null +++ b/ikarus/python/test/testdirichletvalues.py @@ -0,0 +1,75 @@ +# SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de +# SPDX-License-Identifier: LGPL-3.0-or-later + +import debug_info + +debug_info.setDebugFlags() + +import ikarus as iks +import dune.grid +import dune.functions + +import math + + +def makeGrid(): + lowerLeft = [] + upperRight = [] + elements = [] + for _ in range(2): + lowerLeft.append(-1) + upperRight.append(1) + elements.append(3) + + return dune.grid.structuredGrid(lowerLeft, upperRight, elements) + + +def testDirichletValues(): + grid = makeGrid() + basis = iks.basis(grid, dune.functions.Power(dune.functions.Lagrange(order=1), 2)) + basis = basis.flat() + + dirichletValues = iks.dirichletValues(basis) + + assert basis.size() == dirichletValues.size + + def fixOneIndex(vec, globalIndex): + if globalIndex == 0: + vec[globalIndex] = True + + dirichletValues.fixBoundaryDOFs(fixOneIndex) + + # 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 + def fixAnotherIndexWithLocalView(vec, localIndex, localView): + if localView.index(localIndex) == [4]: + vec[localView.index(localIndex)] = True + + assert isinstance(localIndex, int) + assert isinstance(localView.index(localIndex), list) + assert isinstance(localView.index(localIndex)[0], int) + + dirichletValues.fixBoundaryDOFs(fixAnotherIndexWithLocalView) + assert sum(dirichletValues.container) == 2 + + def fixTopSide(vec, localIndex, localView, intersection): + if intersection.geometry.center[0] == 1.0: + vec[localView.index(localIndex)] = True + + assert isinstance(localIndex, int) + assert isinstance(localView.index(localIndex), list) + assert isinstance(localView.index(localIndex)[0], int) + + dirichletValues.fixBoundaryDOFs(fixTopSide) + + # This assmues a structured grid + indicesPerDirection: int = (math.sqrt(grid.size(0)) + 1) * 2 + assert sum(dirichletValues.container) == 2 + indicesPerDirection + + ssb0 = dune.functions.subspaceBasis(basis, 0) + dirichletValues.fixBoundaryDOFsOfSubSpaceBasis(fixTopSide, ssb0) + +if __name__ == "__main__": + testDirichletValues() diff --git a/ikarus/python/test/testmaterials.py b/ikarus/python/test/testmaterials.py index 82ecfdae1..1b150ec30 100644 --- a/ikarus/python/test/testmaterials.py +++ b/ikarus/python/test/testmaterials.py @@ -2,6 +2,7 @@ # SPDX-License-Identifier: LGPL-3.0-or-later import debug_info + debug_info.setDebugFlags() import ikarus as iks diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index 18ac9da4e..e611c6050 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -102,6 +102,33 @@ public: } else DUNE_THROW(Dune::IOError, "Invalid callback function passed to fixBoundaryDOFs"); } + /** + * \brief Function to fix (set boolean values to true or false) degrees of freedom on the boundary using a + * SubSpaceBasis of the Basis used to create the DirichletValues + * + This function takes a callback function, which will be called with the boolean vector of fixed boundary + * degrees of freedom and the usual arguments of `Dune::Functions::forEachBoundaryDOF`. + * + * \param f A callback function + * \param ssb A subspace basis (compile-time-checked to be an actual subspace basis of basis) + */ + template + requires(std::is_same_v::RootBasis, Basis>) + void fixBoundaryDOFs(F&& f, SSB&& ssb) { + if constexpr (Concepts::IsFunctorWithArgs) { + auto lambda = [&](auto&& indexGlobal) { f(dirichletFlagsBackend_, indexGlobal); }; + Dune::Functions::forEachBoundaryDOF(std::forward(ssb), lambda); + } else if constexpr (Concepts::IsFunctorWithArgs) { + auto lambda = [&](auto&& localIndex, auto&& localView) { f(dirichletFlagsBackend_, localIndex, localView); }; + Dune::Functions::forEachBoundaryDOF(std::forward(ssb), lambda); + } else if constexpr (Concepts::IsFunctorWithArgs) { + auto lambda = [&](auto&& localIndex, auto&& localView, auto&& intersection) { + f(dirichletFlagsBackend_, localIndex, localView, intersection); + }; + Dune::Functions::forEachBoundaryDOF(std::forward(ssb), lambda); + } + } /** * \brief Function to fix (set boolean values to true or false) degrees of freedom. diff --git a/tests/src/testdirichletvalue.cpp b/tests/src/testdirichletvalue.cpp index 2559dc6ed..e9495c88a 100644 --- a/tests/src/testdirichletvalue.cpp +++ b/tests/src/testdirichletvalue.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -59,6 +60,27 @@ static auto dirichletBCTest() { "result. Index: i=" << i; + // Test subspacebasis + auto subBasis0 = Dune::Functions::subspaceBasis(basisP->flat(), Dune::Indices::_0); + auto subBasis1 = Dune::Functions::subspaceBasis(basisP->flat(), Dune::Indices::_1); + Ikarus::DirichletValues dirichletValues_SSB(basisP->flat()); + dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, + subBasis0); + + t.check(dirichletValues2.fixedDOFsize() == dirichletValues_SSB.fixedDOFsize() * 2) + << "DirichletValues with subspace basis should have half as many fixed DOFs as the full basis, but has " + << dirichletValues_SSB.fixedDOFsize() << ", where as the full basis has " << dirichletValues2.fixedDOFsize() + << " fixed DOFs"; + + dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, + subBasis1); + + for (std::size_t i = 0; i < basisP->flat().size(); ++i) + t.check(dirichletValues_SSB.isConstrained(i) == dirichletValues2.isConstrained(i)) + << "Different dirichlet value creations with subspace basis didn't provide the same result as with full basis. " + "Index: i=" + << i; + auto inhomogeneousDisplacement = [](const auto& globalCoord, const T& lambda) { Eigen::Vector localInhomogeneous; if (globalCoord[0] > 4 - 1e-8) { From 5477a8aa8025821addb358afacc33b0a637fb555 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Thu, 27 Jun 2024 14:03:31 +0200 Subject: [PATCH 02/18] support for power basis subspace basis dirichlet fixing fix changelog fixBoundaryDOFs now with treepath fix reworked python interface --- CHANGELOG.md | 1 + .../python/dirichletvalues/dirichletvalues.hh | 150 ++++++++---------- ikarus/python/test/testdirichletvalues.py | 25 ++- ikarus/utils/dirichletvalues.hh | 42 ++--- python/ikarus/dirichlet_values.py | 34 +++- tests/src/testdirichletvalue.cpp | 63 ++++++-- 6 files changed, 181 insertions(+), 134 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cfcf894fd..3c9462f6e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,6 +52,7 @@ SPDX-License-Identifier: LGPL-3.0-or-later ([#304](https://github.com/ikarus-project/ikarus/pull/304)) - This can be used, for instance, to apply concentrated forces or to add spring stiffness in a particular direction. - Furthermore, a helper function to get the global index of a Lagrange node at the given global position is added. +- Rework Python Interface for `DirichletValues` plus adding support to easily fix boundary DOFs of subspacebasis in C++ and Python ([#305](https://github.com/ikarus-project/ikarus/pull/305)) ## Release v0.4 (Ganymede) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 5836aff84..7b9dcec9d 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -24,19 +24,70 @@ // PYBIND11_MAKE_OPAQUE(std::vector); namespace Ikarus::Python { +template +void forwardCorrectFunction(DirichletValues& self, const pybind11::function& functor, CppVisitor&& cppFunction) { + using Basis = typename DirichletValues::Basis; + using Intersection = typename Basis::GridView::Intersection; + using BackendType = typename DirichletValues::BackendType; + using MultiIndex = typename Basis::MultiIndex; + + using LocalViewWrapper = Dune::Python::LocalViewWrapper; + + using FixBoundaryDOFsWithGlobalIndexFunction = std::function>, int)>; + using FixBoundaryDOFsWithLocalViewFunction = + std::function>, int, LocalViewWrapper&)>; + using FixBoundaryDOFsWithIntersectionFunction = + std::function>, int, LocalViewWrapper&, const Intersection&)>; + + // Disambiguate by number of arguments, as casting doesn't properly work with functions + pybind11::module inspect_module = pybind11::module::import("inspect"); + pybind11::object result = inspect_module.attr("signature")(functor).attr("parameters"); + size_t numParams = pybind11::len(result); + + if (numParams == 2) { + auto& function = functor.template cast(); + 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(); + 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(); + auto lambda = [&](BackendType& vec, int localIndex, auto& lv, const Intersection& intersection) { + auto lvWrapper = LocalViewWrapper(lv.rootLocalView()); + function(vec.vector(), localIndex, lvWrapper, intersection); + }; + cppFunction(lambda); + + } else { + DUNE_THROW(Dune::NotImplemented, "fixBoundaryDOFs: A function with this signature is not supported"); + } +} + /** * \brief Register Python bindings for a DirichletValues class. * * 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 using a user-defined function `f`. - * - `fixBoundaryDOFsUsingLocalView(f)`: Fixes boundary degrees of freedom using a user-defined function `f` with a - * `LocalView` argument. - * - `fixBoundaryDOFsUsingLocalViewAndIntersection(f)`: Fixes boundary degrees of freedom 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 + * - `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 * * \tparam DirichletValues The DirichletValues class to be registered. * \tparam options Variadic template parameters for additional options when defining the Python class. @@ -61,84 +112,19 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_ LocalViewWrapper; auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; - auto lv = Dune::Python::insertClass( - scopedf, "LocalView", - Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes) - .first; + auto lv_ = Dune::Python::insertClass( + scopedf, "LocalView", + Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes) + .first; cls.def(pybind11::init([](const Basis& basis) { return new DirichletValues(basis); }), pybind11::keep_alive<1, 2>()); - auto fixBoundaryDOFs_ = [](DirichletValues& self, - const std::function>, int)>& f) { - auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { - // we explicitly only allow flat indices - f(vec.vector(), indexGlobal[0]); - }; - self.fixBoundaryDOFs(lambda); - }; - - auto fixBoundaryDOFsUsingLocalView_ = - [](DirichletValues& self, - const std::function>, int, LocalViewWrapper&)>& f) { - auto lambda = [&](BackendType& vec, int localIndex, LocalView& lv) { - auto lvWrapper = LocalViewWrapper(lv.globalBasis()); - // this can be simplified when - // https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/418 becomes available - pybind11::object obj = pybind11::cast(lv.element()); - lvWrapper.bind(obj); - f(vec.vector(), localIndex, lvWrapper); - }; - self.fixBoundaryDOFs(lambda); - }; - - auto fixBoundaryDOFsUsingLocalViewAndIntersection_ = - [](DirichletValues& self, - const std::function>, int, LocalViewWrapper&, const Intersection&)>& f) { - auto lambda = [&](BackendType& vec, int localIndex, LocalView& lv, const Intersection& intersection) { - auto lvWrapper = LocalViewWrapper(lv.globalBasis()); - // this can be simplified when - // https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/418 becomes available - pybind11::object obj = pybind11::cast(lv.element()); - lvWrapper.bind(obj); - f(vec.vector(), localIndex, lvWrapper, intersection); - }; - self.fixBoundaryDOFs(lambda); - }; - - cls.def( - "fixBoundaryDOFs", - [&](DirichletValues& self, const pybind11::function& functor) { - using FixBoundaryDOFsWithGlobalIndexFunction = std::function>, int)>; - using FixBoundaryDOFsWithLocalViewFunction = - std::function>, int, LocalViewWrapper&)>; - using FixBoundaryDOFsWithIntersectionFunction = - std::function>, int, LocalViewWrapper&, const Intersection&)>; - - // Disambiguate by number of arguments, as casting doesn't properly work with functions - pybind11::module inspect_module = pybind11::module::import("inspect"); - pybind11::object result = inspect_module.attr("signature")(functor).attr("parameters"); - size_t numParams = pybind11::len(result); - - if (numParams == 2) { - auto& function = functor.template cast(); - fixBoundaryDOFs_(self, function); - - } else if (numParams == 3) { - auto& function = functor.template cast(); - fixBoundaryDOFsUsingLocalView_(self, function); - - } else if (numParams == 4) { - auto& function = functor.template cast(); - fixBoundaryDOFsUsingLocalViewAndIntersection_(self, function); - - } else { - DUNE_THROW(Dune::NotImplemented, "fixBoundaryDOFs: A function with this signature is not supported"); - } - }, - pybind11::arg("functor")); - - cls.def("fixBoundaryDOFsOfSubSpaceBasis", [](DirichletValues& self, const pybind11::function& functor, - const pybind11::object& ssb) { pybind11::print(ssb); }); + cls.def_property_readonly("container", &DirichletValues::container); + 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("fixDOFs", [](DirichletValues& self, const std::function>)>& f) { @@ -148,8 +134,6 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_ int { return self.size(); }); } } // namespace Ikarus::Python diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py index 8d2168d5d..5102b32ae 100644 --- a/ikarus/python/test/testdirichletvalues.py +++ b/ikarus/python/test/testdirichletvalues.py @@ -3,7 +3,7 @@ import debug_info -debug_info.setDebugFlags() +debug_info.unsetDebugFlags() import ikarus as iks import dune.grid @@ -53,6 +53,7 @@ def fixAnotherIndexWithLocalView(vec, localIndex, localView): dirichletValues.fixBoundaryDOFs(fixAnotherIndexWithLocalView) assert sum(dirichletValues.container) == 2 + assert dirichletValues.fixedDOFsize == 2 def fixTopSide(vec, localIndex, localView, intersection): if intersection.geometry.center[0] == 1.0: @@ -66,10 +67,26 @@ def fixTopSide(vec, localIndex, localView, intersection): # This assmues a structured grid indicesPerDirection: int = (math.sqrt(grid.size(0)) + 1) * 2 - assert sum(dirichletValues.container) == 2 + indicesPerDirection + assert dirichletValues.fixedDOFsize == 2 + indicesPerDirection + + # Test Subbasis + dirichletValues2 = iks.dirichletValues(basis) + + dirichletValues2.fixBoundaryDOFs(fixOneIndex, 0) + dirichletValues2.fixBoundaryDOFs(fixAnotherIndexWithLocalView, 0) + dirichletValues2.fixBoundaryDOFs(fixTopSide, 0) + + assert dirichletValues2.fixedDOFsize == int(2 + indicesPerDirection / 2) + + dirichletValues2.fixBoundaryDOFs(fixTopSide, 1) + assert dirichletValues2.fixedDOFsize == 2 + indicesPerDirection + + dirichletValues2.fixIthDOF(1) + assert dirichletValues2.fixedDOFsize == 2 + indicesPerDirection + 1 + assert dirichletValues2.container[1] + assert dirichletValues2.isConstrained(1) + - ssb0 = dune.functions.subspaceBasis(basis, 0) - dirichletValues.fixBoundaryDOFsOfSubSpaceBasis(fixTopSide, ssb0) if __name__ == "__main__": testDirichletValues() diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index e611c6050..4436cabbe 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -18,6 +18,7 @@ #include #include +#include #include @@ -84,49 +85,24 @@ public: * degrees of freedom and the usual arguments of `Dune::Functions::forEachBoundaryDOF`. * * \param f A callback function + * \param treePath An optional argument specifying a tree path to a subspacebasis, e.g. Dune::Indices::_0 */ - template - void fixBoundaryDOFs(F&& f) { - if constexpr (Concepts::IsFunctorWithArgs) { - auto lambda = [&](auto&& indexGlobal) { f(dirichletFlagsBackend_, indexGlobal); }; - Dune::Functions::forEachBoundaryDOF(basis_, lambda); - } else if constexpr (Concepts::IsFunctorWithArgs) { - auto lambda = [&](auto&& localIndex, auto&& localView) { f(dirichletFlagsBackend_, localIndex, localView); }; - Dune::Functions::forEachBoundaryDOF(basis_, lambda); - } else if constexpr (Concepts::IsFunctorWithArgs) { - auto lambda = [&](auto&& localIndex, auto&& localView, auto&& intersection) { - f(dirichletFlagsBackend_, localIndex, localView, intersection); - }; - Dune::Functions::forEachBoundaryDOF(basis_, lambda); - } else - DUNE_THROW(Dune::IOError, "Invalid callback function passed to fixBoundaryDOFs"); - } - /** - * \brief Function to fix (set boolean values to true or false) degrees of freedom on the boundary using a - * SubSpaceBasis of the Basis used to create the DirichletValues - * - This function takes a callback function, which will be called with the boolean vector of fixed boundary - * degrees of freedom and the usual arguments of `Dune::Functions::forEachBoundaryDOF`. - * - * \param f A callback function - * \param ssb A subspace basis (compile-time-checked to be an actual subspace basis of basis) - */ - template - requires(std::is_same_v::RootBasis, Basis>) - void fixBoundaryDOFs(F&& f, SSB&& ssb) { + template > + void fixBoundaryDOFs(F&& f, TreePath&& treePath = {}) { + using namespace Dune::Functions; + if constexpr (Concepts::IsFunctorWithArgs) { auto lambda = [&](auto&& indexGlobal) { f(dirichletFlagsBackend_, indexGlobal); }; - Dune::Functions::forEachBoundaryDOF(std::forward(ssb), lambda); + Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward(treePath)), lambda); } else if constexpr (Concepts::IsFunctorWithArgs) { auto lambda = [&](auto&& localIndex, auto&& localView) { f(dirichletFlagsBackend_, localIndex, localView); }; - Dune::Functions::forEachBoundaryDOF(std::forward(ssb), lambda); + Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward(treePath)), lambda); } else if constexpr (Concepts::IsFunctorWithArgs) { auto lambda = [&](auto&& localIndex, auto&& localView, auto&& intersection) { f(dirichletFlagsBackend_, localIndex, localView, intersection); }; - Dune::Functions::forEachBoundaryDOF(std::forward(ssb), lambda); + Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward(treePath)), lambda); } } diff --git a/python/ikarus/dirichlet_values.py b/python/ikarus/dirichlet_values.py index ee61fa5f1..b6ceb67a9 100644 --- a/python/ikarus/dirichlet_values.py +++ b/python/ikarus/dirichlet_values.py @@ -1,9 +1,35 @@ # SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de # SPDX-License-Identifier: LGPL-3.0-or-later +import types from dune.common.hashit import hashIt from .generator import MySimpleGenerator +from io import StringIO +from dune.generator.algorithm import run + + +def __fixBoundaryDOFs(dirichletValues): + def __fixBoundaryDOFsFunc(dirichletValues, f, *args: int): + prefixPathTypeName = "Dune::TypeTree::HybridTreePath<" + prefixPathTypeName += ",".join( + "Dune::index_constant<" + str(i) + ">" for i in args + ) + prefixPathTypeName += ">" + + runCode = """ + #include + template + void fixBoundaryDofs(DV& dirichletValues, const pybind11::function& functor) + {{ + Ikarus::Python::forwardCorrectFunction(dirichletValues, functor, [&](auto&& functor_) {{ return dirichletValues.fixBoundaryDOFs(functor_, {prefixPathType}{{}}); }}); + }} + """.format(prefixPathType=prefixPathTypeName) + return run("fixBoundaryDofs", StringIO(runCode), dirichletValues, f) + + return __fixBoundaryDOFsFunc + + def dirichletValues(basis): """ @brief Creates a Dirichlet values handler for the given basis. @@ -21,6 +47,10 @@ def dirichletValues(basis): includes += ["ikarus/python/dirichletvalues/dirichletvalues.hh"] moduleName = "dirichletValues_" + hashIt(element_type) module = generator.load( - includes=includes, typeName=element_type, moduleName=moduleName + includes=includes, typeName=element_type, moduleName=moduleName, dynamicAttr=True ) - return module.DirichletValues(basis) + + dv = module.DirichletValues(basis) + dv.fixBoundaryDOFs = types.MethodType(__fixBoundaryDOFs(dv), dv) + + return dv diff --git a/tests/src/testdirichletvalue.cpp b/tests/src/testdirichletvalue.cpp index e9495c88a..50801d71b 100644 --- a/tests/src/testdirichletvalue.cpp +++ b/tests/src/testdirichletvalue.cpp @@ -14,7 +14,6 @@ #include #include #include -#include #include #include @@ -60,23 +59,63 @@ static auto dirichletBCTest() { "result. Index: i=" << i; - // Test subspacebasis - auto subBasis0 = Dune::Functions::subspaceBasis(basisP->flat(), Dune::Indices::_0); - auto subBasis1 = Dune::Functions::subspaceBasis(basisP->flat(), Dune::Indices::_1); - Ikarus::DirichletValues dirichletValues_SSB(basisP->flat()); - dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, - subBasis0); + auto testSubSpaceBasis = [&](auto&& tree1, auto&& tree2) { + Ikarus::DirichletValues dirichletValues_SSB(basisP->flat()); + dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, + tree1); + + t.check(dirichletValues2.fixedDOFsize() == dirichletValues_SSB.fixedDOFsize() * 2) + << "DirichletValues with subspace basis should have half as many fixed DOFs as the full basis, but has " + << dirichletValues_SSB.fixedDOFsize() << ", where as the full basis has " << dirichletValues2.fixedDOFsize() + << " fixed DOFs"; + + dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, + tree2); + + for (std::size_t i = 0; i < basisP->flat().size(); ++i) + t.check(dirichletValues_SSB.isConstrained(i) == dirichletValues2.isConstrained(i)) + << "Different dirichlet value creations with subspace basis didn't provide the same result as with full " + "basis. " + "Index: i=" + << i; + }; + + // Static tree path + auto treePath0 = Dune::Indices::_0; + auto treePath1 = Dune::Indices::_1; + testSubSpaceBasis(treePath0, treePath1); + + // Dynamic tree path + auto treePath0d = 0; + auto subBasis1d = 1; + testSubSpaceBasis(treePath0d, subBasis1d); + + // Test with Intersection + auto fixLambda = [](auto& dirichletFlags, auto&& localIndex, auto&& localView, auto&& intersection) { + if (intersection.geometry().center()[0] > 4 - 1e-8) + dirichletFlags[localView.index(localIndex)] = true; + }; + Ikarus::DirichletValues dirichletValues5(basisP->flat()); + dirichletValues5.fixBoundaryDOFs(fixLambda); + + Ikarus::DirichletValues dirichletValues_SSB2(basisP->flat()); + dirichletValues_SSB2.fixBoundaryDOFs(fixLambda, treePath0); - t.check(dirichletValues2.fixedDOFsize() == dirichletValues_SSB.fixedDOFsize() * 2) + t.check(dirichletValues5.fixedDOFsize() == dirichletValues_SSB2.fixedDOFsize() * 2) << "DirichletValues with subspace basis should have half as many fixed DOFs as the full basis, but has " - << dirichletValues_SSB.fixedDOFsize() << ", where as the full basis has " << dirichletValues2.fixedDOFsize() + << dirichletValues_SSB2.fixedDOFsize() << ", where as the full basis has " << dirichletValues5.fixedDOFsize() << " fixed DOFs"; - dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, - subBasis1); + dirichletValues_SSB2.fixBoundaryDOFs(fixLambda, treePath1); + + t.check(dirichletValues_SSB2.fixedDOFsize() > 0); + t.check(dirichletValues5.fixedDOFsize() == dirichletValues_SSB2.fixedDOFsize()) + << "DirichletValues with subspace basis should have as many fixed DOFs as the full basis, but has " + << dirichletValues_SSB2.fixedDOFsize() << ", where as the full basis has " << dirichletValues5.fixedDOFsize() + << " fixed DOFs"; for (std::size_t i = 0; i < basisP->flat().size(); ++i) - t.check(dirichletValues_SSB.isConstrained(i) == dirichletValues2.isConstrained(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; From a4b78496024264ab9b377b4506f884f44aeec8d4 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Thu, 27 Jun 2024 15:52:34 +0200 Subject: [PATCH 03/18] add some more interface methods to dirichletvalues --- .../python/dirichletvalues/dirichletvalues.hh | 3 +- ikarus/utils/dirichletvalues.hh | 17 ++++++++- python/ikarus/dirichlet_values.py | 5 ++- tests/src/testdirichletvalue.cpp | 37 ++++++++++++++++++- 4 files changed, 58 insertions(+), 4 deletions(-) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 7b9dcec9d..450b5b9d6 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -25,7 +25,8 @@ namespace Ikarus::Python { template -void forwardCorrectFunction(DirichletValues& self, const pybind11::function& functor, CppVisitor&& cppFunction) { +void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::function& functor, + CppVisitor&& cppFunction) { using Basis = typename DirichletValues::Basis; using Intersection = typename Basis::GridView::Intersection; using BackendType = typename DirichletValues::BackendType; diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index 4436cabbe..2c30ab999 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -120,12 +120,27 @@ public: } /** - * \brief Function to fix (set boolean values to true or false) of degrees 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 */ void fixIthDOF(typename Basis::MultiIndex i) { dirichletFlagsBackend_[i] = true; } + /** + * \brief Unfixes (set boolean value to false) a specific degree of freedom + * + * \param i An index indicating the DOF number to be fixed + */ + void unfixIthDOF(typename Basis::MultiIndex i) { dirichletFlagsBackend_[i] = false; } + + /** + * \brief Resets all degrees of freedom + */ + void reset() { + std::fill(dirichletFlags_.begin(), dirichletFlags_.end(), false); + inhomogeneousBoundaryVectorDummy_.setZero(static_cast(basis_.size())); + } + /* \brief Returns the local basis object */ const auto& basis() const { return basis_; } diff --git a/python/ikarus/dirichlet_values.py b/python/ikarus/dirichlet_values.py index b6ceb67a9..3a2e8c063 100644 --- a/python/ikarus/dirichlet_values.py +++ b/python/ikarus/dirichlet_values.py @@ -22,7 +22,10 @@ def __fixBoundaryDOFsFunc(dirichletValues, f, *args: int): template void fixBoundaryDofs(DV& dirichletValues, const pybind11::function& functor) {{ - Ikarus::Python::forwardCorrectFunction(dirichletValues, functor, [&](auto&& functor_) {{ return dirichletValues.fixBoundaryDOFs(functor_, {prefixPathType}{{}}); }}); + Ikarus::Python::forwardCorrectFunction(dirichletValues, functor, + [&](auto&& functor_) {{ + return dirichletValues.fixBoundaryDOFs(functor_, {prefixPathType}{{}}); + }}); }} """.format(prefixPathType=prefixPathTypeName) return run("fixBoundaryDofs", StringIO(runCode), dirichletValues, f) diff --git a/tests/src/testdirichletvalue.cpp b/tests/src/testdirichletvalue.cpp index 50801d71b..3ae6e0d09 100644 --- a/tests/src/testdirichletvalue.cpp +++ b/tests/src/testdirichletvalue.cpp @@ -29,7 +29,7 @@ using Dune::TestSuite; static auto dirichletBCTest() { - TestSuite t("dirichletBCTest"); + TestSuite t("DirichletValueTest"); using Grid = Dune::YaspGrid<2>; const double Lx = 4.0; @@ -120,6 +120,41 @@ static auto dirichletBCTest() { "Index: i=" << i; + auto sum = [](const auto& container_) { + return std::accumulate(container_.begin(), container_.end(), 0, std::plus()); + }; + auto manual_sum = [](const auto& dv) { + int sum = 0; + for (auto i : Dune::range(dv.size())) + if (dv.isConstrained(i)) + ++sum; + return sum; + }; + + // Test container + auto& container = dirichletValues2.container(); + static_assert(std::is_reference_v, "container should be a reference"); + static_assert(std::is_const_v>, "container should be const"); + + t.check(std::ranges::all_of(container, [](auto&& flag) { return flag; })) << "All values should be fixed"; + + auto sumPre = sum(container); + auto sumManual = manual_sum(dirichletValues2); + + t.check(sumManual == sumPre) + << "Summing over container should yield the same amount of fixed DOFs then checking manually"; + + dirichletValues2.unfixIthDOF({1}); + 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"; + + // 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"; + auto inhomogeneousDisplacement = [](const auto& globalCoord, const T& lambda) { Eigen::Vector localInhomogeneous; if (globalCoord[0] > 4 - 1e-8) { From 659e1bee6b29c935af1834b20f447046f92d5382 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Wed, 3 Jul 2024 15:42:20 +0200 Subject: [PATCH 04/18] docs refs wording additions from review wording some work --- docs/website/01_framework/dirichletBCs.md | 6 ++- .../python/dirichletvalues/dirichletvalues.hh | 47 ++++++++++--------- ikarus/python/test/testdirichletvalues.py | 22 +++++++-- ikarus/utils/dirichletvalues.hh | 12 +++-- tests/src/testdirichletvalue.cpp | 45 +++++++++++------- 5 files changed, 85 insertions(+), 47 deletions(-) diff --git a/docs/website/01_framework/dirichletBCs.md b/docs/website/01_framework/dirichletBCs.md index 0db59b14f..45083e31b 100644 --- a/docs/website/01_framework/dirichletBCs.md +++ b/docs/website/01_framework/dirichletBCs.md @@ -24,11 +24,12 @@ 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. @@ -36,11 +37,12 @@ auto size() const ; // (8)! 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 diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 450b5b9d6..0759e8dc8 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -24,9 +24,8 @@ // PYBIND11_MAKE_OPAQUE(std::vector); namespace Ikarus::Python { -template -void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::function& functor, - CppVisitor&& cppFunction) { +template +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(); - auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { + auto function = functor.template cast(); + 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(); - auto lambda = [&](BackendType& vec, int localIndex, auto& lv) { + auto function = functor.template cast(); + 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(); - auto lambda = [&](BackendType& vec, int localIndex, auto& lv, const Intersection& intersection) { + auto function = functor.template cast(); + 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_()); - 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>)>& f) { diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py index 5102b32ae..4d04020a8 100644 --- a/ikarus/python/test/testdirichletvalues.py +++ b/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__": diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index 2c30ab999..04378c1fc 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -123,15 +123,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 + requires(not std::integral) + 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 diff --git a/tests/src/testdirichletvalue.cpp b/tests/src/testdirichletvalue.cpp index 3ae6e0d09..31433e297 100644 --- a/tests/src/testdirichletvalue.cpp +++ b/tests/src/testdirichletvalue.cpp @@ -2,8 +2,6 @@ // SPDX-License-Identifier: LGPL-3.0-or-later #include -#include "testhelpers.hh" - #include #include @@ -43,24 +41,25 @@ static auto dirichletBCTest() { using namespace Dune::Functions::BasisFactory; auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>())); - auto basisP = std::make_shared(basis); + auto basisP = std::make_shared(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 (dirichletValues1 and dirichletValues2) 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); @@ -95,10 +94,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) @@ -114,15 +113,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())) @@ -131,6 +128,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, "container should be a reference"); @@ -144,17 +153,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 = [](const auto& globalCoord, const T& lambda) { Eigen::Vector localInhomogeneous; if (globalCoord[0] > 4 - 1e-8) { @@ -252,7 +265,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); } } From 7570209056d011be78dcfc0077cbac5ff6ac2b19 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Wed, 3 Jul 2024 15:49:59 +0200 Subject: [PATCH 05/18] fix --- tests/src/testdirichletvalue.cpp | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/tests/src/testdirichletvalue.cpp b/tests/src/testdirichletvalue.cpp index 31433e297..3b7739d1a 100644 --- a/tests/src/testdirichletvalue.cpp +++ b/tests/src/testdirichletvalue.cpp @@ -128,18 +128,6 @@ 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, "container should be a reference"); From fbe112ec61517f9a7613fc082bdf7266c9288f99 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Mon, 15 Jul 2024 10:07:19 +0200 Subject: [PATCH 06/18] code review --- ikarus/python/dirichletvalues/dirichletvalues.hh | 14 +++++--------- ikarus/python/test/testdirichletvalues.py | 4 ++-- ikarus/utils/dirichletvalues.hh | 14 ++++++++++++-- tests/src/testdirichletvalue.cpp | 13 ++++++------- 4 files changed, 25 insertions(+), 20 deletions(-) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 0759e8dc8..48c99e774 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -46,10 +46,7 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu if (numParams == 2) { auto function = functor.template cast(); - auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { - // we explicitly only allow flat indices - function(vec.vector(), indexGlobal[0]); - }; + auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { function(vec.vector(), indexGlobal); }; cppFunction(lambda); } else if (numParams == 3) { @@ -125,19 +122,18 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_()); - cls.def_property_readonly("container", &DirichletValues::container, pybind11::return_value_policy::reference); + cls.def_property_readonly("container", &DirichletValues::container); 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("setSingleDOF", [](DirichletValues& self, std::size_t i, bool flag) { self.setSingleDOF(i, flag); }); + cls.def("isConstrained", [](DirichletValues& self, MultiIndex i) -> bool { return self.isConstrained(i); }); + cls.def("setSingleDOF", [](DirichletValues& self, MultiIndex i, bool flag) { self.setSingleDOF(i, flag); }); cls.def("reset", &DirichletValues::reset); cls.def("fixDOFs", [](DirichletValues& self, const std::function>)>& f) { - auto lambda = [&](const Basis& basis, BackendType& vec) { - // we explicitly only allow flat indices - f(basis, vec.vector()); - }; + auto lambda = [&](const Basis& basis, BackendType& vec) { f(basis, vec.vector()); }; self.fixDOFs(lambda); }); } diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py index 4d04020a8..20dd6a86d 100644 --- a/ikarus/python/test/testdirichletvalues.py +++ b/ikarus/python/test/testdirichletvalues.py @@ -91,9 +91,9 @@ def fixTopSide(vec, localIndex, localView, intersection): assert dirichletValues2.container[1] assert dirichletValues2.isConstrained(1) - dirichletValues2.setSingleDOF(1, False) + dirichletValues2.setSingleDOF((1), False) # via MultiIndex assert dirichletValues2.fixedDOFsize == 2 + indicesPerDirection - assert not dirichletValues2.isConstrained(1) + assert not dirichletValues2.isConstrained((1)) # via MultiIndex dirichletValues2.reset() assert dirichletValues2.fixedDOFsize == 0 diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index 04378c1fc..032213828 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -18,6 +18,7 @@ #include #include +#include #include #include @@ -137,7 +138,12 @@ public: * \param i An index indicating the DOF number to be fixed * \param flag Boolean indicating whether the DOF should fixed or not */ - void setSingleDOF(std::size_t i, bool flag) { dirichletFlags_[i] = flag; } + + void setSingleDOF(std::size_t i, bool flag) + requires(std::same_as>) + { + dirichletFlags_[i] = flag; + } /** * \brief Resets all degrees of freedom @@ -158,7 +164,11 @@ public: } /* \brief Returns a boolean values, if the i-th degree of freedom is constrained */ - [[nodiscard]] bool isConstrained(std::size_t i) const { return dirichletFlags_[i]; } + [[nodiscard]] bool isConstrained(std::size_t i) const + requires(std::same_as>) + { + return dirichletFlags_[i]; + } /* \brief Returns how many degrees of freedoms are fixed */ auto fixedDOFsize() const { return std::ranges::count(dirichletFlags_, true); } diff --git a/tests/src/testdirichletvalue.cpp b/tests/src/testdirichletvalue.cpp index 3b7739d1a..1a1ad5984 100644 --- a/tests/src/testdirichletvalue.cpp +++ b/tests/src/testdirichletvalue.cpp @@ -41,8 +41,7 @@ static auto dirichletBCTest() { using namespace Dune::Functions::BasisFactory; auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>())); - auto basisP = std::make_shared(basis); - auto flatBasis = basisP->flat(); + auto flatBasis = basis.flat(); Ikarus::DirichletValues dirichletValues1(flatBasis); dirichletValues1.fixDOFs([](auto& basis_, auto& dirichFlags) { @@ -71,7 +70,7 @@ static auto dirichletBCTest() { dirichletValues_SSB.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }, tree2); - for (std::size_t i = 0; i < basisP->flat().size(); ++i) + for (std::size_t i = 0; i < flatBasis.size(); ++i) t.check(dirichletValues_SSB.isConstrained(i) == dirichletValues2.isConstrained(i)) << "Different dirichlet value creations with subspace basis didn't provide the same result as with full " "basis. " @@ -194,7 +193,7 @@ static auto dirichletBCTest() { << "Values differ dispDerivs[i]: " << dispDerivs[globalIndex0]; } }; - Dune::Functions::forEachBoundaryDOF(basisP->flat(), lambdaCheck); + Dune::Functions::forEachBoundaryDOF(flatBasis, lambdaCheck); // Check that we can store lambda from python std::string inhomogeneousDisplacementFunction = @@ -229,7 +228,7 @@ static auto dirichletBCTest() { return Dune::toEigen(pythonFunc(globalCoord, lambda_)); }; // - Ikarus::DirichletValues dirichletValues3(basisP->flat()); + Ikarus::DirichletValues dirichletValues3(flatBasis); dirichletValues3.fixBoundaryDOFs([](auto& dirichFlags, auto&& indexGlobal) { dirichFlags[indexGlobal] = true; }); Eigen::VectorXd disps2, dispDerivs2; @@ -240,7 +239,7 @@ static auto dirichletBCTest() { t.check(disps2.isApprox(dispDerivs2 * lambda3)); // Check if all boundary DOFs found manually using obtainLagrangeNodePositions is the same as using forEachBoundaryDOF - Ikarus::DirichletValues dirichletValues4(basisP->flat()); + Ikarus::DirichletValues dirichletValues4(flatBasis); constexpr double tol = 1e-8; auto localView = basis.flat().localView(); for (auto& ele : elements(gridView)) { @@ -257,7 +256,7 @@ static auto dirichletBCTest() { } } - 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) == dirichletValues4.isConstrained(i)) << "Different dirichlet value creations (dirichletValues1 and dirichletValues4) didn't provide the same " "result. Index: i=" From 065c65b9b25c98429d8bd2dfa6debefc9baf81fd Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Mon, 15 Jul 2024 11:01:09 +0200 Subject: [PATCH 07/18] added size to assembler --- ikarus/python/assembler/flatassembler.hh | 3 +++ ikarus/python/dirichletvalues/dirichletvalues.hh | 3 ++- ikarus/python/test/testdirichletvalues.py | 1 + 3 files changed, 6 insertions(+), 1 deletion(-) diff --git a/ikarus/python/assembler/flatassembler.hh b/ikarus/python/assembler/flatassembler.hh index 00da4fbc1..977077652 100644 --- a/ikarus/python/assembler/flatassembler.hh +++ b/ikarus/python/assembler/flatassembler.hh @@ -114,6 +114,8 @@ void registerFlatAssembler(pybind11::handle scope, pybind11::class_ int { return self.size(); }); } #define MAKE_ASSEMBLER_REGISTERY_FUNCTION(name) \ @@ -122,6 +124,7 @@ void registerFlatAssembler(pybind11::handle scope, pybind11::class_>, int, LocalViewWrapper&, const Intersection&)>; - // Disambiguate by number of arguments, as casting doesn't properly work with functions + // Disambiguate by number of arguments pybind11::module inspect_module = pybind11::module::import("inspect"); pybind11::object result = inspect_module.attr("signature")(functor).attr("parameters"); size_t numParams = pybind11::len(result); @@ -124,6 +124,7 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_ int { return self.size(); }); cls.def_property_readonly("fixedDOFsize", &DirichletValues::fixedDOFsize); cls.def("isConstrained", [](DirichletValues& self, std::size_t i) -> bool { return self.isConstrained(i); }); cls.def("setSingleDOF", [](DirichletValues& self, std::size_t i, bool flag) { self.setSingleDOF(i, flag); }); diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py index 20dd6a86d..22fb46649 100644 --- a/ikarus/python/test/testdirichletvalues.py +++ b/ikarus/python/test/testdirichletvalues.py @@ -32,6 +32,7 @@ def testDirichletValues(): dirichletValues = iks.dirichletValues(basis) assert basis.size() == dirichletValues.size + assert basis.size() == len(dirichletValues) def fixOneIndex(vec, globalIndex): if globalIndex == 0: From 03da2ea545243ebbbc6e63b14c9899098ef72995 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 16 Jul 2024 12:16:17 +0200 Subject: [PATCH 08/18] not working python --- .../python/dirichletvalues/dirichletvalues.hh | 68 ++++++++++++------- ikarus/python/test/testdirichletvalues.py | 1 - ikarus/utils/dirichletvalues.hh | 6 +- 3 files changed, 49 insertions(+), 26 deletions(-) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 17be6e6f5..6c47df835 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -24,6 +24,29 @@ // PYBIND11_MAKE_OPAQUE(std::vector); namespace Ikarus::Python { +namespace Impl { + using FixBoundaryDOFsWithGlobalIndexFunction = std::function>, int)>; + + template + using FixBoundaryDOFsWithLocalViewFunction = std::function>, int, LV&)>; + + template + using FixBoundaryDOFsWithIntersectionFunction = + std::function>, int, LV&, const IS&)>; + + template + auto registerLocalView() { + pybind11::module scopedf = pybind11::module::import("dune.functions"); + using LocalViewWrapper = Dune::Python::LocalViewWrapper; + + auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; + return Dune::Python::insertClass( + scopedf, "LocalView", + Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes) + .first; + } +} // namespace Impl + template void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::function& functor, auto&& cppFunction) { using Basis = typename DirichletValues::Basis; @@ -31,13 +54,7 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu using BackendType = typename DirichletValues::BackendType; using MultiIndex = typename Basis::MultiIndex; - using LocalViewWrapper = Dune::Python::LocalViewWrapper; - - using FixBoundaryDOFsWithGlobalIndexFunction = std::function>, int)>; - using FixBoundaryDOFsWithLocalViewFunction = - std::function>, int, LocalViewWrapper&)>; - using FixBoundaryDOFsWithIntersectionFunction = - std::function>, int, LocalViewWrapper&, const Intersection&)>; + // using LocalViewWrapper = Dune::Python::LocalViewWrapper; // Disambiguate by number of arguments pybind11::module inspect_module = pybind11::module::import("inspect"); @@ -45,22 +62,34 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu size_t numParams = pybind11::len(result); if (numParams == 2) { - auto function = functor.template cast(); + auto function = functor.template cast(); auto lambda = [&](BackendType& vec, const MultiIndex& indexGlobal) { function(vec.vector(), indexGlobal); }; cppFunction(lambda); } else if (numParams == 3) { - auto function = functor.template cast(); - auto lambda = [&](BackendType& vec, int localIndex, auto& lv) { - auto lvWrapper = LocalViewWrapper(lv.rootLocalView()); + auto lambda = [&](BackendType& vec, int localIndex, auto&& lv) { + using SubSpaceBasis = typename std::remove_cvref_t::GlobalBasis; + Impl::registerLocalView(); + + using SubSpaceLocalViewWrapper = Dune::Python::LocalViewWrapper; + auto lvWrapper = SubSpaceLocalViewWrapper(lv); + + auto function = + functor.template cast>(); function(vec.vector(), localIndex, lvWrapper); }; cppFunction(lambda); } else if (numParams == 4) { - auto function = functor.template cast(); - auto lambda = [&](BackendType& vec, int localIndex, auto& lv, const Intersection& intersection) { - auto lvWrapper = LocalViewWrapper(lv.rootLocalView()); + auto lambda = [&](BackendType& vec, int localIndex, auto&& lv, const Intersection& intersection) { + using SubSpaceBasis = typename std::remove_cvref_t::GlobalBasis; + Impl::registerLocalView(); + + using SubSpaceLocalViewWrapper = Dune::Python::LocalViewWrapper; + auto lvWrapper = SubSpaceLocalViewWrapper(lv); + + auto function = functor.template cast< + const Impl::FixBoundaryDOFsWithIntersectionFunction>(); function(vec.vector(), localIndex, lvWrapper, intersection); }; cppFunction(lambda); @@ -111,20 +140,13 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_ LocalViewWrapper; - auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; - auto lv_ = Dune::Python::insertClass( - scopedf, "LocalView", - Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes) - .first; + auto lv_ = Impl::registerLocalView(); 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("size", &DirichletValues::size); - cls.def("__len__", [](DirichletValues& self) -> int { return self.size(); }); + cls.def("__len__", &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("setSingleDOF", [](DirichletValues& self, std::size_t i, bool flag) { self.setSingleDOF(i, flag); }); diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py index 22fb46649..63f3ccb00 100644 --- a/ikarus/python/test/testdirichletvalues.py +++ b/ikarus/python/test/testdirichletvalues.py @@ -100,6 +100,5 @@ def fixTopSide(vec, localIndex, localView, intersection): assert dirichletValues2.fixedDOFsize == 0 assert sum(dirichletValues2.container) == 0 - if __name__ == "__main__": testDirichletValues() diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index 032213828..4fda645be 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -91,14 +91,16 @@ public: template > void fixBoundaryDOFs(F&& f, TreePath&& treePath = {}) { using namespace Dune::Functions; + using SubSpaceLocalView = + typename std::remove_cvref_t(treePath)))>::LocalView; if constexpr (Concepts::IsFunctorWithArgs) { auto lambda = [&](auto&& indexGlobal) { f(dirichletFlagsBackend_, indexGlobal); }; Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward(treePath)), lambda); - } else if constexpr (Concepts::IsFunctorWithArgs) { + } else if constexpr (Concepts::IsFunctorWithArgs) { auto lambda = [&](auto&& localIndex, auto&& localView) { f(dirichletFlagsBackend_, localIndex, localView); }; Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward(treePath)), lambda); - } else if constexpr (Concepts::IsFunctorWithArgs) { auto lambda = [&](auto&& localIndex, auto&& localView, auto&& intersection) { f(dirichletFlagsBackend_, localIndex, localView, intersection); From f279a25e767e49c61a17af66c8af5860c9213562 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Tue, 16 Jul 2024 15:40:44 +0200 Subject: [PATCH 09/18] Fix? --- .../python/dirichletvalues/dirichletvalues.hh | 40 ++++++++++++++----- 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 6c47df835..aa1dcea59 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -34,16 +35,37 @@ namespace Impl { using FixBoundaryDOFsWithIntersectionFunction = std::function>, int, LV&, const IS&)>; - template + template auto registerLocalView() { pybind11::module scopedf = pybind11::module::import("dune.functions"); - using LocalViewWrapper = Dune::Python::LocalViewWrapper; + using LocalView = Dune::Python::LocalViewWrapper; auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; - return Dune::Python::insertClass( - scopedf, "LocalView", - Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes) - .first; + + // also register subspace basis + if constexpr (registerBasis) { + auto construct = [](const Basis& basis) { return new Basis(basis); }; + + auto [basisCls, isNotRegistered] = Dune::Python::insertClass( + scopedf, "SubspaceBasis", Dune::Python::GenerateTypeName(Dune::className()), includes); + if (isNotRegistered) + Dune::Python::registerSubspaceBasis(scopedf, basisCls); + // Dune::Python::registerBasisType(scopedf, basisCls, construct, std::false_type{}); + } + + auto [lv, isNotRegistered] = Dune::Python::insertClass( + scopedf, "LocalView", Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), + includes); + + if (isNotRegistered) { + lv.def("bind", &LocalView::bind); + lv.def("unbind", &LocalView::unbind); + lv.def("index", [](const LocalView& localView, int index) { return localView.index(index); }); + lv.def("__len__", [](LocalView& self) -> int { return self.size(); }); + + Dune::Python::Functions::registerTree(lv); + lv.def("tree", [](const LocalView& view) { return view.tree(); }); + } } } // namespace Impl @@ -69,7 +91,7 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu } else if (numParams == 3) { auto lambda = [&](BackendType& vec, int localIndex, auto&& lv) { using SubSpaceBasis = typename std::remove_cvref_t::GlobalBasis; - Impl::registerLocalView(); + Impl::registerLocalView(); using SubSpaceLocalViewWrapper = Dune::Python::LocalViewWrapper; auto lvWrapper = SubSpaceLocalViewWrapper(lv); @@ -83,7 +105,7 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu } else if (numParams == 4) { auto lambda = [&](BackendType& vec, int localIndex, auto&& lv, const Intersection& intersection) { using SubSpaceBasis = typename std::remove_cvref_t::GlobalBasis; - Impl::registerLocalView(); + Impl::registerLocalView(); using SubSpaceLocalViewWrapper = Dune::Python::LocalViewWrapper; auto lvWrapper = SubSpaceLocalViewWrapper(lv); @@ -140,7 +162,7 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_(); + Impl::registerLocalView(); cls.def(pybind11::init([](const Basis& basis) { return new DirichletValues(basis); }), pybind11::keep_alive<1, 2>()); From fb2de1d9ffdcc0040897bf8327d0a7af0cb79f40 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Tue, 16 Jul 2024 19:13:10 +0200 Subject: [PATCH 10/18] not a fix --- .../python/dirichletvalues/dirichletvalues.hh | 43 +++++++++++-------- 1 file changed, 25 insertions(+), 18 deletions(-) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index aa1dcea59..3082f7d2a 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -8,6 +8,10 @@ #pragma once +#include +#include + +#include "dune/common/classname.hh" #include #include #include @@ -46,25 +50,28 @@ namespace Impl { if constexpr (registerBasis) { auto construct = [](const Basis& basis) { return new Basis(basis); }; - auto [basisCls, isNotRegistered] = Dune::Python::insertClass( - scopedf, "SubspaceBasis", Dune::Python::GenerateTypeName(Dune::className()), includes); - if (isNotRegistered) - Dune::Python::registerSubspaceBasis(scopedf, basisCls); + // This if statement does absolutly nothing + if (Dune::Python::findInTypeRegistry().second) { + auto [basisCls, isNotRegistered] = Dune::Python::insertClass( + scopedf, "SubspaceBasis", Dune::Python::GenerateTypeName(Dune::className()), includes); + if (isNotRegistered) + Dune::Python::registerSubspaceBasis(scopedf, basisCls); + } // Dune::Python::registerBasisType(scopedf, basisCls, construct, std::false_type{}); - } - - auto [lv, isNotRegistered] = Dune::Python::insertClass( - scopedf, "LocalView", Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), - includes); - - if (isNotRegistered) { - lv.def("bind", &LocalView::bind); - lv.def("unbind", &LocalView::unbind); - lv.def("index", [](const LocalView& localView, int index) { return localView.index(index); }); - lv.def("__len__", [](LocalView& self) -> int { return self.size(); }); - - Dune::Python::Functions::registerTree(lv); - lv.def("tree", [](const LocalView& view) { return view.tree(); }); + } else { + // auto [lv, isNotRegistered] = Dune::Python::insertClass( + // scopedf, "LocalView", + // Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes); + + // if (isNotRegistered) { + // lv.def("bind", &LocalView::bind); + // lv.def("unbind", &LocalView::unbind); + // lv.def("index", [](const LocalView& localView, int index) { return localView.index(index); }); + // lv.def("__len__", [](LocalView& self) -> int { return self.size(); }); + + // Dune::Python::Functions::registerTree(lv); + // lv.def("tree", [](const LocalView& view) { return view.tree(); }); + // } } } } // namespace Impl From c93b98f4b31287b242a5cdf9e59d98f0c65bec62 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Wed, 17 Jul 2024 15:21:45 +0200 Subject: [PATCH 11/18] Fix subspacelocalview --- .../python/dirichletvalues/dirichletvalues.hh | 70 ++++++++++--------- ikarus/python/finiteelements/fe.hh | 29 ++++---- .../finiteelements/registerferequirements.hh | 18 ++--- ikarus/python/test/testdirichletvalues.py | 1 + 4 files changed, 62 insertions(+), 56 deletions(-) diff --git a/ikarus/python/dirichletvalues/dirichletvalues.hh b/ikarus/python/dirichletvalues/dirichletvalues.hh index 3082f7d2a..28ff9cf7b 100644 --- a/ikarus/python/dirichletvalues/dirichletvalues.hh +++ b/ikarus/python/dirichletvalues/dirichletvalues.hh @@ -39,39 +39,27 @@ namespace Impl { using FixBoundaryDOFsWithIntersectionFunction = std::function>, int, LV&, const IS&)>; - template - auto registerLocalView() { + template + auto registerSubSpaceLocalView() { pybind11::module scopedf = pybind11::module::import("dune.functions"); - using LocalView = Dune::Python::LocalViewWrapper; + using LocalViewWrapper = Dune::Python::LocalViewWrapper; auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; - // also register subspace basis - if constexpr (registerBasis) { - auto construct = [](const Basis& basis) { return new Basis(basis); }; - - // This if statement does absolutly nothing - if (Dune::Python::findInTypeRegistry().second) { - auto [basisCls, isNotRegistered] = Dune::Python::insertClass( - scopedf, "SubspaceBasis", Dune::Python::GenerateTypeName(Dune::className()), includes); - if (isNotRegistered) - Dune::Python::registerSubspaceBasis(scopedf, basisCls); - } - // Dune::Python::registerBasisType(scopedf, basisCls, construct, std::false_type{}); - } else { - // auto [lv, isNotRegistered] = Dune::Python::insertClass( - // scopedf, "LocalView", - // Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes); - - // if (isNotRegistered) { - // lv.def("bind", &LocalView::bind); - // lv.def("unbind", &LocalView::unbind); - // lv.def("index", [](const LocalView& localView, int index) { return localView.index(index); }); - // lv.def("__len__", [](LocalView& self) -> int { return self.size(); }); - - // Dune::Python::Functions::registerTree(lv); - // lv.def("tree", [](const LocalView& view) { return view.tree(); }); - // } + Dune::Python::insertClass(scopedf, "SubspaceBasis_" + Dune::className(), + Dune::Python::GenerateTypeName(Dune::className()), includes); + + auto [lv, isNew] = Dune::Python::insertClass( + scopedf, "LocalView_" + Dune::className(), + Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), includes); + if (isNew) { + lv.def("bind", &LocalViewWrapper::bind); + lv.def("unbind", &LocalViewWrapper::unbind); + lv.def("index", [](const LocalViewWrapper& localView, int index) { return localView.index(index); }); + lv.def("__len__", [](LocalViewWrapper& self) -> int { return self.size(); }); + + Dune::Python::Functions::registerTree(lv); + lv.def("tree", [](const LocalViewWrapper& view) { return view.tree(); }); } } } // namespace Impl @@ -83,8 +71,6 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu using BackendType = typename DirichletValues::BackendType; using MultiIndex = typename Basis::MultiIndex; - // using LocalViewWrapper = Dune::Python::LocalViewWrapper; - // Disambiguate by number of arguments pybind11::module inspect_module = pybind11::module::import("inspect"); pybind11::object result = inspect_module.attr("signature")(functor).attr("parameters"); @@ -98,7 +84,7 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu } else if (numParams == 3) { auto lambda = [&](BackendType& vec, int localIndex, auto&& lv) { using SubSpaceBasis = typename std::remove_cvref_t::GlobalBasis; - Impl::registerLocalView(); + Impl::registerSubSpaceLocalView(); using SubSpaceLocalViewWrapper = Dune::Python::LocalViewWrapper; auto lvWrapper = SubSpaceLocalViewWrapper(lv); @@ -112,7 +98,7 @@ void forwardCorrectFunction(DirichletValues& dirichletValues, const pybind11::fu } else if (numParams == 4) { auto lambda = [&](BackendType& vec, int localIndex, auto&& lv, const Intersection& intersection) { using SubSpaceBasis = typename std::remove_cvref_t::GlobalBasis; - Impl::registerLocalView(); + Impl::registerSubSpaceLocalView(); using SubSpaceLocalViewWrapper = Dune::Python::LocalViewWrapper; auto lvWrapper = SubSpaceLocalViewWrapper(lv); @@ -169,7 +155,23 @@ void registerDirichletValues(pybind11::handle scope, pybind11::class_(); + pybind11::module scopedf = pybind11::module::import("dune.functions"); + using LocalViewWrapper = Dune::Python::LocalViewWrapper; + + auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; + auto [lv, isNew] = Dune::Python::insertClass( + scopedf, "LocalView", Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapper", Dune::MetaType()), + includes); + + if (isNew) { + lv.def("bind", &LocalViewWrapper::bind); + lv.def("unbind", &LocalViewWrapper::unbind); + lv.def("index", [](const LocalViewWrapper& localView, int index) { return localView.index(index); }); + lv.def("__len__", [](LocalViewWrapper& self) -> int { return self.size(); }); + + Dune::Python::Functions::registerTree(lv); + lv.def("tree", [](const LocalViewWrapper& view) { return view.tree(); }); + } cls.def(pybind11::init([](const Basis& basis) { return new DirichletValues(basis); }), pybind11::keep_alive<1, 2>()); diff --git a/ikarus/python/finiteelements/fe.hh b/ikarus/python/finiteelements/fe.hh index 6f0e6a136..34eb02e49 100644 --- a/ikarus/python/finiteelements/fe.hh +++ b/ikarus/python/finiteelements/fe.hh @@ -113,21 +113,22 @@ void registerFE(pybind11::handle scope, pybind11::class_ cls) { pybind11::arg("Requirement"), pybind11::arg("MatrixAffordance"), pybind11::arg("elementMatrix").noconvert()); pybind11::module scopedf = pybind11::module::import("dune.functions"); + using LocalViewWrapper = Dune::Python::LocalViewWrapper; - typedef Dune::Python::LocalViewWrapper LocalViewWrapper; - auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; - auto lv = Dune::Python::insertClass( - scopedf, "LocalViewWrapper", - Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapperWrapper", Dune::MetaType()), - includes) - .first; - lv.def("bind", &LocalViewWrapper::bind); - lv.def("unbind", &LocalViewWrapper::unbind); - lv.def("index", [](const LocalViewWrapper& localView, int index) { return localView.index(index); }); - lv.def("__len__", [](LocalViewWrapper& self) -> int { return self.size(); }); - - Dune::Python::Functions::registerTree(lv); - lv.def("tree", [](const LocalViewWrapper& view) { return view.tree(); }); + auto includes = Dune::Python::IncludeFiles{"dune/python/functions/globalbasis.hh"}; + auto [lv, isNew] = Dune::Python::insertClass( + scopedf, "LocalViewWrapper", + Dune::Python::GenerateTypeName("Dune::Python::LocalViewWrapperWrapper", Dune::MetaType()), includes); + + if (isNew) { + lv.def("bind", &LocalViewWrapper::bind); + lv.def("unbind", &LocalViewWrapper::unbind); + lv.def("index", [](const LocalViewWrapper& localView, int index) { return localView.index(index); }); + lv.def("__len__", [](LocalViewWrapper& self) -> int { return self.size(); }); + + Dune::Python::Functions::registerTree(lv); + lv.def("tree", [](const LocalViewWrapper& view) { return view.tree(); }); + } cls.def( "localView", diff --git a/ikarus/python/finiteelements/registerferequirements.hh b/ikarus/python/finiteelements/registerferequirements.hh index 77db64de7..83a2b531b 100644 --- a/ikarus/python/finiteelements/registerferequirements.hh +++ b/ikarus/python/finiteelements/registerferequirements.hh @@ -3,6 +3,7 @@ #pragma once +#include #include #include @@ -22,18 +23,19 @@ void registerFERequirement(pybind11::handle scope, pybind11::class_( + auto includes = Dune::Python::IncludeFiles{"ikarus/finiteelements/ferequirements.hh"}; + auto [req, isNew] = Dune::Python::insertClass( scope, "FERequirements", Dune::Python::GenerateTypeName(Dune::className()), includes); - if (isNotRegistered) { - lv.def(pybind11::init()); - lv.def(pybind11::init()); - lv.def( + if (isNew) { + req.def(pybind11::init()); + req.def(pybind11::init()); + + req.def( "insertGlobalSolution", [](FERequirements& self, SolutionVectorType solVec) { self.insertGlobalSolution(solVec); }, "solutionVector"_a.noconvert()); - lv.def( + req.def( "globalSolution", [](FERequirements& self) { return self.globalSolution(); }, pybind11::return_value_policy::reference_internal); lv.def( @@ -41,7 +43,7 @@ void registerFERequirement(pybind11::handle scope, pybind11::class_& parVal) { self.insertParameter(parVal.value()); }, pybind11::keep_alive<1, 2>(), "parameterValue"_a.noconvert()); - lv.def("parameter", [](const FERequirements& self) { return self.parameter(); }); + req.def("parameter", [](const FERequirements& self) { return self.parameter(); }); } } } // namespace Ikarus::Python \ No newline at end of file diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/testdirichletvalues.py index 63f3ccb00..22fb46649 100644 --- a/ikarus/python/test/testdirichletvalues.py +++ b/ikarus/python/test/testdirichletvalues.py @@ -100,5 +100,6 @@ def fixTopSide(vec, localIndex, localView, intersection): assert dirichletValues2.fixedDOFsize == 0 assert sum(dirichletValues2.container) == 0 + if __name__ == "__main__": testDirichletValues() From 94329a72c0981a6233f041a3fb74a20674b72211 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 18 Jul 2024 10:14:10 +0200 Subject: [PATCH 12/18] fixes from Taruns Review --- CHANGELOG.md | 1 + docs/website/01_framework/dirichletBCs.md | 12 ++++++------ ikarus/python/test/CMakeLists.txt | 2 +- ...{testdirichletvalues.py => dirichletvaluetest.py} | 0 ikarus/utils/dirichletvalues.hh | 6 ++++-- 5 files changed, 12 insertions(+), 9 deletions(-) rename ikarus/python/test/{testdirichletvalues.py => dirichletvaluetest.py} (100%) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c9462f6e..8e9969bed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -53,6 +53,7 @@ SPDX-License-Identifier: LGPL-3.0-or-later - This can be used, for instance, to apply concentrated forces or to add spring stiffness in a particular direction. - Furthermore, a helper function to get the global index of a Lagrange node at the given global position is added. - Rework Python Interface for `DirichletValues` plus adding support to easily fix boundary DOFs of subspacebasis in C++ and Python ([#305](https://github.com/ikarus-project/ikarus/pull/305)) +- Rework the Python Interface for `DirichletValues` plus add support to easily fix boundary DOFs of `Subspacebasis` in C++ and Python ([#305](https://github.com/ikarus-project/ikarus/pull/305)) ## Release v0.4 (Ganymede) diff --git a/docs/website/01_framework/dirichletBCs.md b/docs/website/01_framework/dirichletBCs.md index 45083e31b..22cd30457 100644 --- a/docs/website/01_framework/dirichletBCs.md +++ b/docs/website/01_framework/dirichletBCs.md @@ -24,25 +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 setSingleDOF(i, bool flag); // (4)! +void setSingleDOF(i, flag); // (4)! const auto& basis() const; // (5)! -bool isConstrained(std::size_t i) const; // (6)! +bool isConstrained(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. +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 +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 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. +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. +8. Returns the number of all Dirichlet degrees of freedom. 9. Resets the whole container \bibliography diff --git a/ikarus/python/test/CMakeLists.txt b/ikarus/python/test/CMakeLists.txt index e49554b13..ddcf6e290 100644 --- a/ikarus/python/test/CMakeLists.txt +++ b/ikarus/python/test/CMakeLists.txt @@ -38,7 +38,7 @@ dune_python_add_test( NAME pydirichletvalues SCRIPT - testdirichletvalues.py + dirichletvaluetest.py WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} LABELS diff --git a/ikarus/python/test/testdirichletvalues.py b/ikarus/python/test/dirichletvaluetest.py similarity index 100% rename from ikarus/python/test/testdirichletvalues.py rename to ikarus/python/test/dirichletvaluetest.py diff --git a/ikarus/utils/dirichletvalues.hh b/ikarus/utils/dirichletvalues.hh index 4fda645be..4591214b8 100644 --- a/ikarus/utils/dirichletvalues.hh +++ b/ikarus/utils/dirichletvalues.hh @@ -106,6 +106,8 @@ public: f(dirichletFlagsBackend_, localIndex, localView, intersection); }; Dune::Functions::forEachBoundaryDOF(subspaceBasis(basis_, std::forward(treePath)), lambda); + } else { + static_assert(Dune::AlwaysFalse(), "fixBoundaryDOFs: A function with this signature is not supported"); } } @@ -123,7 +125,7 @@ public: } /** - * \brief Fixes (set boolean value to true) a specific degree of freedom + * \brief Fixes and unfixes (set boolean value to true or false) 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 @@ -135,7 +137,7 @@ public: } /** - * \brief Fixes (set boolean value to true) a specific degree of freedom + * \brief Fixes or unfixes (set boolean value to true or flase) 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 From c12a365d41d2587d08fab226e5ef8149a9da27c3 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 18 Jul 2024 13:47:41 +0200 Subject: [PATCH 13/18] WIP on feature/dirichletValues --- ikarus/python/assembler/flatassembler.hh | 5 ++--- ikarus/utils/dirichletvalues.hh | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/ikarus/python/assembler/flatassembler.hh b/ikarus/python/assembler/flatassembler.hh index 977077652..444611da2 100644 --- a/ikarus/python/assembler/flatassembler.hh +++ b/ikarus/python/assembler/flatassembler.hh @@ -114,8 +114,8 @@ void registerFlatAssembler(pybind11::handle scope, pybind11::class_ int { return self.size(); }); + cls.def_property_readonly("size", &Assembler::size); + cls.def("__len__", [](Assembler& self) -> int { return self.size(); }); } #define MAKE_ASSEMBLER_REGISTERY_FUNCTION(name) \ @@ -124,7 +124,6 @@ void registerFlatAssembler(pybind11::handle scope, pybind11::class_ Date: Thu, 18 Jul 2024 14:00:18 +0200 Subject: [PATCH 14/18] fix --- tests/src/testassembler.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/src/testassembler.cpp b/tests/src/testassembler.cpp index 798a96d7f..c46c8dfb9 100644 --- a/tests/src/testassembler.cpp +++ b/tests/src/testassembler.cpp @@ -89,7 +89,7 @@ auto SimpleAssemblersTest(const PreBasis& preBasis) { } else { const auto fixIndices = utils::globalIndexFromGlobalPosition(basis.flat(), pos); for (const auto idx : fixIndices) - dirichletValues.fixIthDOF(idx); + dirichletValues.setSingleDOF(idx, true); } using SparseAssembler = SparseFlatAssembler; From 4ba604c9dc4474b1ecda4450b8e8863e17c4bee7 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 18 Jul 2024 14:17:03 +0200 Subject: [PATCH 15/18] fix --- ikarus/python/finiteelements/registerferequirements.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ikarus/python/finiteelements/registerferequirements.hh b/ikarus/python/finiteelements/registerferequirements.hh index 83a2b531b..83f2158ca 100644 --- a/ikarus/python/finiteelements/registerferequirements.hh +++ b/ikarus/python/finiteelements/registerferequirements.hh @@ -38,7 +38,7 @@ void registerFERequirement(pybind11::handle scope, pybind11::class_& parVal) { self.insertParameter(parVal.value()); }, pybind11::keep_alive<1, 2>(), "parameterValue"_a.noconvert()); From deb479b06b7be38888579a6102b5f80525514423 Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 18 Jul 2024 15:02:35 +0200 Subject: [PATCH 16/18] fix? --- ikarus/python/assembler/flatassembler.hh | 2 -- 1 file changed, 2 deletions(-) diff --git a/ikarus/python/assembler/flatassembler.hh b/ikarus/python/assembler/flatassembler.hh index 444611da2..00da4fbc1 100644 --- a/ikarus/python/assembler/flatassembler.hh +++ b/ikarus/python/assembler/flatassembler.hh @@ -114,8 +114,6 @@ void registerFlatAssembler(pybind11::handle scope, pybind11::class_ int { return self.size(); }); } #define MAKE_ASSEMBLER_REGISTERY_FUNCTION(name) \ From 83951825334d81251ed351c2b8ae1fc1eac3c10a Mon Sep 17 00:00:00 2001 From: Henrik Jakob Date: Thu, 18 Jul 2024 17:16:34 +0200 Subject: [PATCH 17/18] add test for fixDOFs --- ikarus/python/test/dirichletvaluetest.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ikarus/python/test/dirichletvaluetest.py b/ikarus/python/test/dirichletvaluetest.py index 22fb46649..128382f03 100644 --- a/ikarus/python/test/dirichletvaluetest.py +++ b/ikarus/python/test/dirichletvaluetest.py @@ -100,6 +100,12 @@ def fixTopSide(vec, localIndex, localView, intersection): assert dirichletValues2.fixedDOFsize == 0 assert sum(dirichletValues2.container) == 0 + def fixDOFFunction(basis, vec): + vec[0] = True + + dirichletValues2.fixDOFs(fixDOFFunction) + assert dirichletValues2.fixedDOFsize == 1 + if __name__ == "__main__": testDirichletValues() From 605ae02613c289f6577c1d218da2cfa64b8b6bf5 Mon Sep 17 00:00:00 2001 From: henrij22 <96132706+henrij22@users.noreply.github.com> Date: Thu, 18 Jul 2024 18:54:27 +0200 Subject: [PATCH 18/18] trigger CI --- ikarus/python/test/dirichletvaluetest.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ikarus/python/test/dirichletvaluetest.py b/ikarus/python/test/dirichletvaluetest.py index 128382f03..3ab49a59a 100644 --- a/ikarus/python/test/dirichletvaluetest.py +++ b/ikarus/python/test/dirichletvaluetest.py @@ -10,6 +10,7 @@ import dune.functions import math +import numpy as np def makeGrid(): @@ -101,10 +102,10 @@ def fixTopSide(vec, localIndex, localView, intersection): assert sum(dirichletValues2.container) == 0 def fixDOFFunction(basis, vec): - vec[0] = True + vec[:] = np.ones(dirichletValues2.size) dirichletValues2.fixDOFs(fixDOFFunction) - assert dirichletValues2.fixedDOFsize == 1 + assert dirichletValues2.fixedDOFsize == dirichletValues2.size if __name__ == "__main__":