diff -Nru dune-localfunctions-2.7.1/CHANGELOG.md dune-localfunctions-2.8.0/CHANGELOG.md --- dune-localfunctions-2.7.1/CHANGELOG.md 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/CHANGELOG.md 2021-08-31 14:03:38.000000000 +0000 @@ -1,11 +1,43 @@ -# Release 2.7 +# Release 2.8 -* The header `lagrange.hh` now includes all headers of all Lagrange implementations, - not just the ones with run-time order. +* Passing functions that support `f.evaluate(x,y)` to `interpolate()` + is deprecated. Instead the functions should now provide `operator()`. + Passing functions providing the old interface is still supported in 2.8. + * `LocalFiniteElementFunctionBase` is deprecated. You can rely + on duck-typing when passing functions with the new interface. + * The virtual interface for interpolating functions in `LocalFiniteElementVirtualInterface` + now uses `std::function` instead of the deprecated `VirtualFunction` + for the passed function. + * The virtual interface wrapper `LocalFiniteElementVirtualImp` now + requires that the wrapped `LocalFiniteElement` implementation + supports the new `operator()` based interpolation interface. + +* Add an implementation of the Nédélec element of the first kind, + as introduced in "Nédélec, Mixed finite elements in R^3, 1980, + DOI: http://dx.doi.org/10.1007/BF01396415". + Only the first-order case for triangles, tetrahedra, squares and cubes is implemented. * Fix a bug in a shape function of the second-order Lagrange element on the three-dimensional pyramid. +* Add an implementation of the Raviart-Thomas element for tetrehedra with order 0. + +* Remove deprecated `GenericLocalFiniteElement::topologyId()`, use + `type().id()` instead. + +* Imported the Python bindings from the 2.7 branch of dune-python. + +* Replaced the combination of function arguments `topologyId` and `dim` with a single `GeometryType` argument. + Tagged the old versions of: `numLagrangePoints`, `equidistantLagrangePoints`, `RTL2InterpolationBuilder::topologyId()`, + `VirtualMonomialBasis(topologyId)`, `VirtualMonomialBasis::topologyId()` as deprecated. + +* Add a construction algorithm for high order Nédélec elements on triangles and tetrahedra. + +# Release 2.7 + +* The header `lagrange.hh` now includes all headers of all Lagrange implementations, + not just the ones with run-time order. + * Introduce a run-time polymorphic container `LocalFiniteElementVariant`. Much like `std::variant`, it implements a type-safe union of different `LocalFiniteElement` implementations. Elements of type @@ -50,8 +82,6 @@ * Introduce a convenience header `dune/localfunctions/brezzidouglasmarini.hh` that includes all BDM implementations. -* Fix bugs in the implementation of `interpolate` for Raviart-Thomas elements. - # Release 2.6 * The `diffOrder` value has disappeared from the `LocalBasisTraits` class. diff -Nru dune-localfunctions-2.7.1/CMakeLists.txt dune-localfunctions-2.8.0/CMakeLists.txt --- dune-localfunctions-2.7.1/CMakeLists.txt 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -2,7 +2,7 @@ project("dune-localfunctions" C CXX) # general stuff -cmake_minimum_required(VERSION 3.1) +cmake_minimum_required(VERSION 3.13) # guess dune-common build dir if(NOT (dune-common_DIR OR dune-common_ROOT OR @@ -25,4 +25,10 @@ add_subdirectory(doc) add_subdirectory(dune) +if( DUNE_ENABLE_PYTHONBINDINGS ) + add_subdirectory(python) + dune_python_install_package(PATH python) +endif() + + finalize_dune_project(GENERATE_CONFIG_H_CMAKE) diff -Nru dune-localfunctions-2.7.1/COPYING dune-localfunctions-2.8.0/COPYING --- dune-localfunctions-2.7.1/COPYING 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/COPYING 2021-08-31 14:03:38.000000000 +0000 @@ -7,33 +7,38 @@ 2016 Lukas Böger 2014 Andreas Buhr 2014--2017 Ansgar Burchardt -2009--2019 Andreas Dedner +2009--2021 Andreas Dedner +2020 Nils-Arne Dreier 2010 Martin Drohmann -2008--2018 Christian Engwer -2008--2017 Jorrit Fahlke +2008--2020 Christian Engwer +2008--2020 Jorrit Fahlke 2011--2013 Bernd Flemisch 2015 Elisa Friebel 2012 Christoph Gersbacher 2019 Janick Gerstenberger -2009--2019 Carsten Gräser +2009--2021 Carsten Gräser 2015--2018 Felix Gruber -2011--2019 Christoph Grüninger +2011--2021 Christoph Grüninger +2020 René Heß 2017 Patrick Jaap 2013 Guillaume Jouvet -2015--2019 Dominic Kempf +2015--2020 Dominic Kempf 2015 Angela Klewinghaus +2020 Robert Klöfkorn 2014 Arne Morten Kvarving 2015 Jizhou Li 2013--2017 Tobias Malkmus 2009 Sven Marnach 2013--2018 Steffen Müthing +2020 Lisa Julia Nebel 2012 Rebecca Neumann 2008--2018 Martin Nolte 2011--2016 Elias Pipping -2016--2019 Simon Praetorius +2016--2021 Simon Praetorius 2012--2013 Human Rezaijafari 2011--2012 Uli Sack -2008--2019 Oliver Sander +2008--2021 Oliver Sander +2020--2021 Henrik Stolzmann 2011--2012 Matthias Wohlmuth 2010--2015 Jonathan Youett diff -Nru dune-localfunctions-2.7.1/debian/changelog dune-localfunctions-2.8.0/debian/changelog --- dune-localfunctions-2.7.1/debian/changelog 2021-01-11 21:31:42.000000000 +0000 +++ dune-localfunctions-2.8.0/debian/changelog 2021-11-04 09:39:18.000000000 +0000 @@ -1,3 +1,48 @@ +dune-localfunctions (2.8.0-5ubuntu2) jammy; urgency=medium + + * Use -O1 CXXFLAG instead of gcc-10 + + -- Gianfranco Costamagna Thu, 04 Nov 2021 10:39:18 +0100 + +dune-localfunctions (2.8.0-5) unstable; urgency=medium + + * Upload to unstable. + + -- Ansgar Thu, 21 Oct 2021 18:40:46 +0200 + +dune-localfunctions (2.8.0-4) experimental; urgency=medium + + * Increase tolerance for test-orthonormal even more. This should + also fix the build failure on armel. + + -- Ansgar Tue, 19 Oct 2021 22:19:02 +0200 + +dune-localfunctions (2.8.0-3) experimental; urgency=medium + + * d/patches: Increase tolerance for test-orthonormal even more. + + -- Patrick Jaap Mon, 18 Oct 2021 13:25:43 +0200 + +dune-localfunctions (2.8.0-2) experimental; urgency=medium + + * d/patches: Add patch to increase tolerance for a unit test + + -- Patrick Jaap Fri, 15 Oct 2021 16:07:35 +0200 + +dune-localfunctions (2.8.0-1) experimental; urgency=medium + + * New upstream release. + * d/control: Set Architecture to "any" + + -- Patrick Jaap Fri, 08 Oct 2021 13:39:34 +0200 + +dune-localfunctions (2.8.0~rc1-1) experimental; urgency=medium + + * New upstream release. + * debian/rules: Explicitly set buildsystem to CMake + + -- Patrick Jaap Fri, 20 Aug 2021 07:09:59 -0400 + dune-localfunctions (2.7.1-2) unstable; urgency=medium * Upload to unstable. diff -Nru dune-localfunctions-2.7.1/debian/control dune-localfunctions-2.8.0/debian/control --- dune-localfunctions-2.7.1/debian/control 2021-01-07 13:18:47.000000000 +0000 +++ dune-localfunctions-2.8.0/debian/control 2021-11-04 09:39:17.000000000 +0000 @@ -8,19 +8,19 @@ Vcs-Git: https://salsa.debian.org/science-team/dune-localfunctions.git Homepage: https://www.dune-project.org/ Build-Depends: debhelper-compat (= 13), - cmake (>= 3.1), gfortran, mpi-default-bin, mpi-default-dev, pkg-config, python3, - libdune-common-dev (>= 2.7.1), - libdune-geometry-dev (>= 2.7.1) + cmake (>= 3.13), gfortran, mpi-default-bin, mpi-default-dev, pkg-config, python3, + libdune-common-dev (>= 2.8.0), + libdune-geometry-dev (>= 2.8.0) Build-Depends-Indep: doxygen, ghostscript, graphviz, imagemagick, texlive-latex-extra, texlive-latex-recommended, texlive-pictures Rules-Requires-Root: no Package: libdune-localfunctions-dev Section: libdevel -Architecture: all +Architecture: any Multi-Arch: foreign Depends: ${misc:Depends}, - libdune-common-dev (>= 2.7.1), - libdune-geometry-dev (>= 2.7.1) + libdune-common-dev (>= 2.8.0), + libdune-geometry-dev (>= 2.8.0) Suggests: libdune-localfunctions-doc (= ${source:Version}) Description: toolbox for solving PDEs -- local basis (development files) DUNE, the Distributed and Unified Numerics Environment is a modular toolbox diff -Nru dune-localfunctions-2.7.1/debian/copyright dune-localfunctions-2.8.0/debian/copyright --- dune-localfunctions-2.7.1/debian/copyright 2020-05-25 14:22:01.000000000 +0000 +++ dune-localfunctions-2.8.0/debian/copyright 2021-10-13 18:20:01.000000000 +0000 @@ -11,33 +11,38 @@ 2016 Lukas Böger 2014 Andreas Buhr 2014--2017 Ansgar Burchardt - 2009--2019 Andreas Dedner + 2009--2021 Andreas Dedner + 2020 Nils-Arne Dreier 2010 Martin Drohmann - 2008--2018 Christian Engwer - 2008--2017 Jorrit Fahlke + 2008--2020 Christian Engwer + 2008--2020 Jorrit Fahlke 2011--2013 Bernd Flemisch 2015 Elisa Friebel 2012 Christoph Gersbacher 2019 Janick Gerstenberger - 2009--2019 Carsten Gräser + 2009--2021 Carsten Gräser 2015--2018 Felix Gruber - 2011--2019 Christoph Grüninger + 2011--2021 Christoph Grüninger + 2020 René Heß 2017 Patrick Jaap 2013 Guillaume Jouvet - 2015--2019 Dominic Kempf + 2015--2020 Dominic Kempf 2015 Angela Klewinghaus + 2020 Robert Klöfkorn 2014 Arne Morten Kvarving 2015 Jizhou Li 2013--2017 Tobias Malkmus 2009 Sven Marnach 2013--2018 Steffen Müthing + 2020 Lisa Julia Nebel 2012 Rebecca Neumann 2008--2018 Martin Nolte 2011--2016 Elias Pipping - 2016--2019 Simon Praetorius + 2016--2021 Simon Praetorius 2012--2013 Human Rezaijafari 2011--2012 Uli Sack - 2008--2019 Oliver Sander + 2008--2021 Oliver Sander + 2020--2021 Henrik Stolzmann 2011--2012 Matthias Wohlmuth 2010--2015 Jonathan Youett License: GPL-2 with DUNE exception diff -Nru dune-localfunctions-2.7.1/debian/patches/increase-test-orthonormal-tolerance.patch dune-localfunctions-2.8.0/debian/patches/increase-test-orthonormal-tolerance.patch --- dune-localfunctions-2.7.1/debian/patches/increase-test-orthonormal-tolerance.patch 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/debian/patches/increase-test-orthonormal-tolerance.patch 2021-10-19 20:17:30.000000000 +0000 @@ -0,0 +1,13 @@ +diff --git a/dune/localfunctions/test/test-orthonormal.cc b/dune/localfunctions/test/test-orthonormal.cc +index be3b73a..ddadcdd 100644 +--- a/dune/localfunctions/test/test-orthonormal.cc ++++ b/dune/localfunctions/test/test-orthonormal.cc +@@ -73,7 +73,7 @@ bool test(unsigned int order) + for( unsigned int j = 0; j < size; ++j ) + { + const double value = m[ i*size + j ]; +- if( std::abs( value - double( i == j ) ) > 1200.*Dune::Zero::epsilon() ) { ++ if( std::abs( value - double( i == j ) ) > 5000.*Dune::Zero::epsilon() ) { + std::cout << "i = " << i << ", j = " << j << ": " << std::abs( value - double( i == j ) ) << std::endl; + ret = false; + } diff -Nru dune-localfunctions-2.7.1/debian/patches/series dune-localfunctions-2.8.0/debian/patches/series --- dune-localfunctions-2.7.1/debian/patches/series 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/debian/patches/series 2021-10-15 16:18:24.000000000 +0000 @@ -0,0 +1 @@ +increase-test-orthonormal-tolerance.patch diff -Nru dune-localfunctions-2.7.1/debian/rules dune-localfunctions-2.8.0/debian/rules --- dune-localfunctions-2.7.1/debian/rules 2020-05-25 14:22:01.000000000 +0000 +++ dune-localfunctions-2.8.0/debian/rules 2021-11-04 09:39:18.000000000 +0000 @@ -2,5 +2,9 @@ include /usr/share/dune/dune-debian.mk +ifneq (,$(filter ppc64el, $(DEB_HOST_ARCH))) + export CXXFLAGS=-O1 +endif + %: - dh $@ --builddirectory=build + dh $@ --builddirectory=build --buildsystem=cmake diff -Nru dune-localfunctions-2.7.1/dune/CMakeLists.txt dune-localfunctions-2.8.0/dune/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/CMakeLists.txt 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -1 +1,5 @@ add_subdirectory(localfunctions) + +if( DUNE_ENABLE_PYTHONBINDINGS ) + add_subdirectory("python") +endif() diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d/brezzidouglasmarini1simplex2dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d/brezzidouglasmarini1simplex2dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d/brezzidouglasmarini1simplex2dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d/brezzidouglasmarini1simplex2dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -36,6 +36,7 @@ */ BDM1Simplex2DLocalInterpolation (unsigned int s) { + using std::sqrt; sign0 = sign1 = sign2 = 1.0; if (s & 1) { diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/CMakeLists.txt dune-localfunctions-2.8.0/dune/localfunctions/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/localfunctions/CMakeLists.txt 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -7,6 +7,7 @@ add_subdirectory(meta) add_subdirectory(mimetic) add_subdirectory(monomial) +add_subdirectory(nedelec) add_subdirectory(orthonormal) add_subdirectory(rannacherturek) add_subdirectory(raviartthomas) @@ -23,6 +24,7 @@ lagrange.hh mimetic.hh monomial.hh + nedelec.hh orthonormal.hh rannacherturek.hh raviartthomas.hh diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/common/interfaceswitch.hh dune-localfunctions-2.8.0/dune/localfunctions/common/interfaceswitch.hh --- dune-localfunctions-2.7.1/dune/localfunctions/common/interfaceswitch.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/common/interfaceswitch.hh 2021-08-31 14:03:38.000000000 +0000 @@ -9,7 +9,8 @@ #include #include -#include +#include +#include namespace Dune { @@ -44,25 +45,32 @@ //! Type for storing finite elements /** - * Some algorithms use one variable to store (a pointer) a finite element - * and update that pointer while iterating through the grid. This works - * well for local finite elements, since they exists in a finite number of - * variants, which can be stored somewhere and don't need to change for - * the duration of the algorithm, so we can always store a simple pointer. - * For global finite elements we have to store the object itself however, - * and we must make sure that we destroy the object when we are done with - * it. Since global finite elements are not assignable in general, we - * needs to copy-construct them for each grid element we visit. + * Some algorithms use one variable to store (as a shared pointer) + * a finite element and update that pointer while iterating + * through the grid. This works well as long as there is only a + * moderate number of different finite elements, which can be + * stored somewhere and don't change for the duration of the + * algorithm. This is the case for most local finite elements, + * since they exists in a finite number of variants. * - * To accommodate both interfaces, we define a store: for local finite - * elements it is a simple pointer, and if we want to store a finite - * element in it we simply store its address. For global finite elements - * we use a shared_ptr, and if we want to store a finite element in it we - * allocate a new object and initialise it with the copy-constructor. For - * local finite elements we don't need to do anything when we are done - * with it, global finite elements are automatically destructed by the - * shared_ptr when we store a new one or when the shared_ptr itself is - * destroyed. Access to the finite element is done by simply + * If the number of possible finite element realizations grows to + * big, e.g. for global finite elements or also for p-adaptive + * local finite elements, these are only created on the + * fly. Therefore we need to store a copy. Since finite elements + * in general are not assignable, we either copy-construct or + * move-construt them for each grid element we visit. + * + * To accommodate both interfaces, we store in a `shared_ptr`. + * Different ways to initialize a possible, from an l-value + * reference, an r-value reference or from a shared_ptr. + * + * For backwards compatibility we assume that an l-value reference + * to a local finite element is persistent and that we can simply + * store the pointer using `stackobject_to_shared_ptr`, while in + * the case of global finite elements we always need to copy + * construct using `make_shared`. If a local finite element is not + * persistent, it should be passed in as an r-value + * reference. Access to the finite element is done by simply * dereferencing the store in both cases. */ typedef std::shared_ptr Store; @@ -73,7 +81,13 @@ * with allocation and copy-construction and storing that. */ static void setStore(Store& store, const FiniteElement& fe) - { store.reset(new FiniteElement(fe)); } + { store = std::make_shared(fe); } + //! Store a finite element in the store. + static void setStore(Store& store, FiniteElement&& fe) + { store = std::make_shared(std::move(fe)); } + //! Store a finite element in the store. + static void setStore(Store& store, const Store& fe) + { store = fe; } }; #ifndef DOXYGEN @@ -82,7 +96,7 @@ template struct FiniteElementInterfaceSwitch< FiniteElement, - typename std::enable_if::value>::type > { @@ -105,10 +119,16 @@ { return fe.localCoefficients(); } //! Type for storing finite elements - typedef const FiniteElement *Store; + typedef std::shared_ptr Store; //! Store a finite element in the store. static void setStore(Store& store, const FiniteElement& fe) - { store = &fe; } + { store = stackobject_to_shared_ptr(fe); } + //! Store a finite element in the store. + static void setStore(Store& store, FiniteElement&& fe) + { store = std::make_shared(std::move(fe)); } + //! Store a finite element in the store. + static void setStore(Store& store, const Store& fe) + { store = fe; } }; #endif // !DOXYGEN @@ -171,7 +191,7 @@ template struct BasisInterfaceSwitch #include +#include #include #include -#include #include #include @@ -30,8 +30,8 @@ template void visitIf(Visitor&& visitor, Variant&& variant) { - auto visitorWithFallback = overload([&](Std::monostate& impl) {}, [&](const Std::monostate& impl) {}, visitor); - Std::visit(visitorWithFallback, variant); + auto visitorWithFallback = overload([&](std::monostate& impl) {}, [&](const std::monostate& impl) {}, visitor); + std::visit(visitorWithFallback, variant); } template @@ -123,7 +123,7 @@ } private: - Std::variant impl_; + std::variant impl_; std::size_t size_; std::size_t order_; }; @@ -161,13 +161,13 @@ // an l-value reference, we use a default constructed // dummy LocalKey value. static const Dune::LocalKey dummyLocalKey; - return Std::visit(overload( - [&](const Std::monostate& impl) -> decltype(auto) { return (dummyLocalKey);}, + return std::visit(overload( + [&](const std::monostate& impl) -> decltype(auto) { return (dummyLocalKey);}, [&](const auto* impl) -> decltype(auto) { return impl->localKey(i); }), impl_); } private: - Std::variant impl_; + std::variant impl_; std::size_t size_; }; @@ -195,7 +195,7 @@ } private: - Std::variant impl_; + std::variant impl_; }; } // namespace Impl @@ -209,16 +209,14 @@ * implementations that this class can hold have to be provided as * template parameter. * - * The implementation is based on Std::variant - * which is either std::variant or a drop-in replacement if the former is - * not available. - * Notice that this prepends Std::monostate to the Implementations - * list for the internally stored Std::variant such that + * The implementation is based on std::variant. + * Notice that this prepends std::monostate to the Implementations + * list for the internally stored std::variant such that * LocalFiniteElementVariant can be empty and is default-constructible. - * As a consequence providing Std::monostate manually to + * As a consequence providing std::monostate manually to * LocalFiniteElementVariant is neither necessary nor allowed. * Access to the stored implementation is internally implemented - * using Std::visit(). To avoid multiple trivial Std::visit() + * using std::visit(). To avoid multiple trivial std::visit() * calls, the results of size(), order(), and type() are cached * on creation and assignment. * @@ -235,7 +233,7 @@ class LocalFiniteElementVariant { - // In each LocalFooVariant we store a Std::variant, i.e. a Std::variant + // In each LocalFooVariant we store a std::variant, i.e. a std::variant // with the pointer to the Foo implementation unless LocalFiniteElementVariant stores a monostate. In this // case each LocalFooVariant also stores a monostate (and not a monostate*). using LocalBasis = Impl::LocalBasisVariant; @@ -245,8 +243,8 @@ // Update members after changing impl_ void updateMembers() { - Std::visit(overload( - [&](Std::monostate&) { + std::visit(overload( + [&](std::monostate&) { localBasis_ = LocalBasis(); localCoefficients_ = LocalCoefficients(); localInterpolation_ = LocalInterpolation(); @@ -276,7 +274,7 @@ /** * \brief Construct empty LocalFiniteElementVariant */ - LocalFiniteElementVariant(const Std::monostate& monostate) + LocalFiniteElementVariant(const std::monostate& monostate) {} /** @@ -286,7 +284,7 @@ * copy of the provided implementation. */ template, Implementations>...>::value, int> = 0> + std::enable_if_t, Implementations>...>::value, int> = 0> LocalFiniteElementVariant(Implementation&& impl) : impl_(std::forward(impl)) { @@ -335,7 +333,7 @@ * \brief Assignment from implementation */ template, Implementations>...>::value, int> = 0> + std::enable_if_t, Implementations>...>::value, int> = 0> LocalFiniteElementVariant& operator=(Implementation&& impl) { impl_ = std::forward(impl); @@ -385,15 +383,15 @@ } /** - * \brief Provide access to underlying Std::variant + * \brief Provide access to underlying std::variant * - * This allows to use Std::visit on a higher level + * This allows to use std::visit on a higher level * which allows to avoid the indirection of the - * Std::variant - polymorphism inside the visitor code. - * Notice that the provided Std::variant contains - * Std::monostate in its type list. Hence any + * std::variant - polymorphism inside the visitor code. + * Notice that the provided std::variant contains + * std::monostate in its type list. Hence any * visitor used to access the variant has to be - * Std::monostate-aware. + * std::monostate-aware. */ const auto& variant() const { @@ -407,11 +405,11 @@ */ operator bool () const { - return not(Std::holds_alternative(variant())); + return not(std::holds_alternative(variant())); } private: - Std::variant impl_; + std::variant impl_; std::size_t size_; GeometryType geometryType_; LocalBasis localBasis_; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/common/localinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/common/localinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/common/localinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/common/localinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -6,7 +6,6 @@ #include #include -#include @@ -34,24 +33,6 @@ ); }; - // Create function supporting f.evaluate(Domain, Range&) - // If the argument already does this, just forward it. - template, F>(), int> = 0> - decltype(auto) makeFunctionWithEvaluate(const F& f) - { - return f; - } - - // Create function supporting f.evaluate(Domain, Range&) - // If the argument does not support this, wrap it as VirtualFunction - template, F>(), int> = 0> - decltype(auto) makeFunctionWithEvaluate(const F& f) - { - return makeVirtualFunction(std::cref(f)); - } - // Create function supporting Range = f(Domain) // If the argument already does this, just forward it. template >, F>(), int> = 0> +#ifndef DUNE_DEPRECATED_INTERPOLATE_CHECK + [[deprecated( "Passing functions only supporting 'f.evaluate(x,y)' to interpolate() is deprecated." + "Use functions supporting operator(), i.e. f(x) instead!")]] +#endif decltype(auto) makeFunctionWithCallOperator(const F& f) { return [&](auto&& x) { diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/common/virtualinterface.hh dune-localfunctions-2.8.0/dune/localfunctions/common/virtualinterface.hh --- dune-localfunctions-2.7.1/dune/localfunctions/common/virtualinterface.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/common/virtualinterface.hh 2021-08-31 14:03:38.000000000 +0000 @@ -6,8 +6,7 @@ #include #include #include - -#include +#include #include @@ -31,27 +30,52 @@ * @brief Return a proper base class for functions to use with LocalInterpolation. * * @tparam FE A FiniteElement type + * + * \deprecated + * This class is deprecated. + * To keep this traits class working it exports a simple + * look-a-like of the old Dune::Function base class. + * However, you should stop using this and pass functions with + * plain operator() interface to interpolate() from now on. */ template - class LocalFiniteElementFunctionBase + class + [[deprecated("Dune::LocalFiniteElementFunctionBase is deprecated after Dune 2.7. You can now pass functions providing operator() to interpolate.")]] + LocalFiniteElementFunctionBase { - typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType; - typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType; + typedef typename FE::Traits::LocalBasisType::Traits::DomainType Domain; + typedef typename FE::Traits::LocalBasisType::Traits::RangeType Range; + + // Hack: Keep a copy of Dune::Function here. This allows to avoid depending + // on the deprecated dune-common header while still keeping the LocalFiniteElementFunctionBase + // mechanism working during its deprecation period. + class FunctionBaseDummy + { + public: + + using RangeType = Range; + using DomainType = Domain; + + struct Traits + { + using RangeType = Range; + using DomainType = Domain; + }; - typedef LocalInterpolationVirtualInterface Interface; - typedef typename FE::Traits::LocalInterpolationType Implementation; + void evaluate(const DomainType& x, RangeType& y) const; + }; public: - typedef VirtualFunction VirtualFunctionBase; - typedef Function FunctionBase; + using VirtualFunctionBase = FunctionBaseDummy; + using FunctionBase = FunctionBaseDummy; /** \brief Base class type for functions to use with LocalInterpolation * - * This is the VirtualFunction interface class if FE implements the virtual + * This is just a dummy providing the old typedefs. * interface and Function base class otherwise. */ - typedef typename std::conditional::value, VirtualFunctionBase, FunctionBase>::type type; + using type = FunctionBaseDummy; }; @@ -133,8 +157,8 @@ { public: - //! type of virtual function to interpolate - typedef Dune::VirtualFunction FunctionType; + //! type of function to interpolate + using FunctionType = std::function; //! type of the coefficient vector in the interpolate method typedef typename RangeType::field_type CoefficientType; @@ -164,8 +188,8 @@ { public: - //! type of virtual function to interpolate - typedef Dune::VirtualFunction FunctionType; + //! type of function to interpolate + using FunctionType = std::function; //! type of the coefficient vector in the interpolate method typedef typename RangeType::field_type CoefficientType; @@ -196,7 +220,7 @@ const auto& f = Impl::makeFunctionWithCallOperator(ff); const LocalInterpolationVirtualInterfaceBase& asBase = *this; - asBase.interpolate(makeVirtualFunction(std::cref(f)),out); + asBase.interpolate(FunctionType(std::cref(f)),out); } /** \brief determine coefficients interpolating a given function @@ -211,7 +235,7 @@ std::vector outDummy; const LocalInterpolationVirtualInterfaceBase& asBase = *this; - asBase.interpolate(makeVirtualFunction(std::cref(f)),outDummy); + asBase.interpolate(FunctionType(std::cref(f)),outDummy); out.resize(outDummy.size()); for(typename std::vector::size_type i=0; i -#include - #include #include #include @@ -139,8 +137,8 @@ * @brief class for wrapping a local interpolation * using the virtual interface * - * @tparam DomainType domain type of the Dune::VirtualFunction to interpolate - * @tparam RangeType range type of the Dune::VirtualFunction to interpolate + * @tparam DomainType domain type of the function to interpolate + * @tparam RangeType range type of the function to interpolate * \tparam Imp LocalInterpolationVirtualInterface implementation */ template diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/emptypoints.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/emptypoints.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/emptypoints.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/emptypoints.hh 2021-08-31 14:03:38.000000000 +0000 @@ -36,8 +36,14 @@ return localKey_; } + const Field weight () const + { + return weight_; + } + Vector point_; LocalKey localKey_; + Field weight_; }; // EmptyPointSet diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/equidistantpoints.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/equidistantpoints.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/equidistantpoints.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/equidistantpoints.hh 2021-08-31 14:03:38.000000000 +0000 @@ -6,7 +6,7 @@ #include #include -#include +#include #include #include @@ -18,52 +18,57 @@ // numLagrangePoints // ----------------- - inline std::size_t numLagrangePoints ( unsigned int topologyId, int dim, std::size_t order ) + inline std::size_t numLagrangePoints ( const GeometryType& gt, std::size_t order ) { - assert( topologyId < Impl::numTopologies( dim ) ); - + const int dim = gt.dim(); if( dim > 0 ) { - const unsigned int baseId = Impl::baseTopologyId( topologyId, dim ); - if( Impl::isPyramid( topologyId, dim ) ) + const GeometryType baseGeometryType = Impl::getBase( gt ); + if( gt.isConical() ) { std::size_t size = 0; for( unsigned int o = 0; o <= order; ++o ) - size += numLagrangePoints( baseId, dim-1, o ); + size += numLagrangePoints( baseGeometryType, o ); return size; } else - return numLagrangePoints( baseId, dim-1, order ) * (order+1); + return numLagrangePoints( baseGeometryType, order ) * (order+1); } else return 1; } + [[deprecated("Use numLagrangePoints(const GeometryType& gt, std::size_t order ) instead.")]] + inline std::size_t numLagrangePoints ( unsigned int topologyId, unsigned int dim, std::size_t order ) + { + return numLagrangePoints ( GeometryType(topologyId, dim), order); + } + // equidistantLagrangePoints // ------------------------- template< class ct, unsigned int cdim > - inline static unsigned int equidistantLagrangePoints ( unsigned int topologyId, unsigned int dim, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points ) + inline static unsigned int equidistantLagrangePoints ( const GeometryType& gt, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points ) { + const unsigned int dim = gt.dim(); assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) ); - assert( topologyId < Impl::numTopologies( dim ) ); if( dim > 0 ) { - const unsigned int baseId = Impl::baseTopologyId( topologyId, dim ); - const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseId, dim-1, codim ) : 0); - const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseId, dim-1, codim-1 ) : 0); + const GeometryType baseGeometryType = Impl::getBase( gt ); + const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim ) : 0); + const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim-1 ) : 0); - if( Impl::isPrism( topologyId, dim ) ) + if( gt.isPrismatic() ) { unsigned int size = 0; if( codim < dim ) { for( unsigned int i = 1; i < order; ++i ) { - const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim, order, count, points ); + const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, order, count, points ); for( unsigned int j = 0; j < n; ++j ) { LocalKey &key = points->localKey_; @@ -77,7 +82,7 @@ if( codim > 0 ) { - const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim-1, order, count+numBaseN, points ); + const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim-1, order, count+numBaseN, points ); for( unsigned int j = 0; j < n; ++j ) { LocalKey &key = points[ j ].localKey_; @@ -95,7 +100,7 @@ } else { - unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseId, dim-1, codim-1, order, count, points ) : 0); + unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseGeometryType, codim-1, order, count, points ) : 0); LagrangePoint< ct, cdim > *const end = points + size; for( ; points != end; ++points ) points->localKey_ = LocalKey( points->localKey_.subEntity(), codim, points->localKey_.index() ); @@ -104,7 +109,7 @@ { for( unsigned int i = order-1; i > 0; --i ) { - const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim, i, count+numBaseM, points ); + const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, i, count+numBaseM, points ); LagrangePoint< ct, cdim > *const end = points + n; for( ; points != end; ++points ) { @@ -135,6 +140,13 @@ } } + template< class ct, unsigned int cdim > + [[deprecated("Use equidistantLagrangePoints ( GeometryType gt, ... ) instead.")]] + inline static unsigned int equidistantLagrangePoints ( unsigned int topologyId, unsigned int dim, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points ) + { + return equidistantLagrangePoints ( GeometryType(topologyId, dim), codim, order, *count, *points ); + } + // EquidistantPointSet @@ -156,7 +168,7 @@ void build ( GeometryType gt ) { assert( gt.dim() == dimension ); - points_.resize( numLagrangePoints( gt.id(), dimension, order() ) ); + points_.resize( numLagrangePoints( gt, order() ) ); typename Base::LagrangePoint *p = points_.data(); std::vector< unsigned int > count; @@ -164,19 +176,31 @@ { count.resize( Geo::Impl::size( gt.id(), dimension, dimension-mydim ) ); std::fill( count.begin(), count.end(), 0u ); - p += equidistantLagrangePoints( gt.id(), dimension, dimension-mydim, order(), count.data(), p ); + p += equidistantLagrangePoints( gt, dimension-mydim, order(), count.data(), p ); } + const auto &refElement = referenceElement(gt); + F weight = refElement.volume()/F(double(points_.size())); + for (auto &p : points_) + p.weight_ = weight; } - template< class T > + template< GeometryType::Id geometryId > bool build () { - build( GeometryType( T() ) ); + build( GeometryType( geometryId ) ); return true; } - template< class T > - static bool supports ( std::size_t order ) { return true; } + bool buildCube () + { + return build< GeometryTypes::cube(dim) > (); + } + + static bool supports ( GeometryType gt, std::size_t order ) { return true; } + template< GeometryType::Id geometryId> + static bool supports ( std::size_t order ) { + return supports( GeometryType( geometryId ), order ); + } private: using Base::points_; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/interpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/interpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/interpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/interpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -7,7 +7,6 @@ #include #include -#include #include #include @@ -46,7 +45,7 @@ template< class Fn, class Vector > auto interpolate ( const Fn &fn, Vector &coefficients, PriorityTag< 1 > ) const - -> std::enable_if_t< Std::is_invocable< const Fn &, decltype( this->lagrangePoints_.begin()->point() ) >::value > + -> std::enable_if_t< std::is_invocable_v< const Fn &, decltype( this->lagrangePoints_.begin()->point() ) > > { unsigned int index = 0; for( const auto &lp : lagrangePoints_ ) @@ -103,17 +102,17 @@ typedef typename LagrangePointSetFactory::Key Key; typedef const LocalLagrangeInterpolation< LP,dim,F > Object; - template< class Topology > + template< GeometryType::Id geometryId > static Object *create ( const Key &key ) { const LagrangePointSet *lagrangeCoeff - = LagrangePointSetFactory::template create< Topology >( key ); + = LagrangePointSetFactory::template create< geometryId >( key ); if ( lagrangeCoeff == 0 ) return 0; else return new Object( *lagrangeCoeff ); } - template< class Topology > + template< GeometryType::Id geometryId > static bool supports ( const Key &key ) { return true; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/lagrangecoefficients.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/lagrangecoefficients.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/lagrangecoefficients.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/lagrangecoefficients.hh 2021-08-31 14:03:38.000000000 +0000 @@ -23,14 +23,14 @@ const typedef LP Object; typedef std::size_t Key; - template< class T > + template< GeometryType::Id geometryId > static Object *create ( const Key &order ) { - if (order == 0 || !Object::template supports(order)) + if (order == 0 || !Object::template supports(order)) return 0; typedef typename std::remove_const::type LagrangeCoefficients; LagrangeCoefficients *object = new LagrangeCoefficients(order); - if ( !object->template build() ) + if ( !object->template build() ) { delete object; object = nullptr; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/lagrangesimplex.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/lagrangesimplex.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/lagrangesimplex.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/lagrangesimplex.hh 2021-08-31 14:03:38.000000000 +0000 @@ -400,7 +400,7 @@ // Helper method: Return the derivative of a single Lagrangian factor of l_ij evaluated at x // direction: Derive in x-direction if this is 0, otherwise derive in y direction auto lagrangianFactorDerivative = [&lagrangeNode] - (const int direction, const int no, const int i, const int j, const typename Traits::DomainType& x) + (const int direction, const int no, const int i, const int j, const typename Traits::DomainType&) -> typename Traits::RangeType { using T = typename Traits::RangeType; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/p0/p0localbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/p0/p0localbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/p0/p0localbasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/p0/p0localbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -38,7 +38,7 @@ } //! \brief Evaluate all shape functions - inline void evaluateFunction (const typename Traits::DomainType& in, + inline void evaluateFunction (const typename Traits::DomainType&, std::vector& out) const { out.resize(1); @@ -47,7 +47,7 @@ //! \brief Evaluate Jacobian of all shape functions inline void - evaluateJacobian (const typename Traits::DomainType& in, // position + evaluateJacobian (const typename Traits::DomainType&, // position std::vector& out) const // return value { out.resize(1); diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/p0/p0localcoefficients.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/p0/p0localcoefficients.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/p0/p0localcoefficients.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/p0/p0localcoefficients.hh 2021-08-31 14:03:38.000000000 +0000 @@ -32,7 +32,7 @@ } //! get i'th index - const LocalKey& localKey (std::size_t i) const + const LocalKey& localKey ([[maybe_unused]] std::size_t i) const { return index; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/lagrange/pqkfactory.hh dune-localfunctions-2.8.0/dune/localfunctions/lagrange/pqkfactory.hh --- dune-localfunctions-2.7.1/dune/localfunctions/lagrange/pqkfactory.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/lagrange/pqkfactory.hh 2021-08-31 14:03:38.000000000 +0000 @@ -32,7 +32,7 @@ typedef typename P0LocalFiniteElement::Traits::LocalBasisType::Traits T; //! create finite element for given GeometryType - static LocalFiniteElementVirtualInterface* create(const GeometryType& gt) + static LocalFiniteElementVirtualInterface* create(const GeometryType&) { return nullptr; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/meta/power/interpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/meta/power/interpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/meta/power/interpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/meta/power/interpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -8,7 +8,6 @@ #include #include #include -#include #include namespace Dune { @@ -43,8 +42,7 @@ private: template - class ComponentEvaluator : - public Function + class ComponentEvaluator { const F &f; std::size_t comp; @@ -54,11 +52,12 @@ f(f_), comp(comp_) { } - void evaluate(const typename Backend::Traits::DomainLocal &x, - typename Backend::Traits::Range &y) const + typename Backend::Traits::Range operator()(const typename Backend::Traits::DomainLocal &x) const { typename Traits::Range fy = f(x); + typename Backend::Traits::Range y; y[0] = fy[comp]; + return y; } }; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/CMakeLists.txt dune-localfunctions-2.8.0/dune/localfunctions/nedelec/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,4 @@ +install(FILES + nedelec1stkindsimplex.hh + nedelec1stkindcube.hh + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/localfunctions/nedelec) diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelec1stkindcube.hh dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelec1stkindcube.hh --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelec1stkindcube.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelec1stkindcube.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,550 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH +#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include // For deprecated makeFunctionWithCallOperator +#include + +namespace Dune +{ +namespace Impl +{ + /** \brief Nedelec shape functions of the first kind on reference cubes + * + * \note These shape functions are implemented for the reference cubes only! + * The transformation to other cubes has to be done by the user. + * + * \tparam D Number type used for domain coordinates + * \tparam R Number type used for function values + * \tparam dim Dimension of the reference element + * \tparam k Order of the Nedelec element + */ + template + class Nedelec1stKindCubeLocalBasis + { + // Number of edges of the reference cube + constexpr static std::size_t numberOfEdges = power(2,dim-1)*dim; + + public: + using Traits = LocalBasisTraits, + R,dim,FieldVector, + FieldMatrix >; + + /** \brief Constructor with default edge orientation + * + * Default orientation means that all edges point from the vertex with lower index + * to the vertex with higher index. The tangential parts of all shape functions + * point in the direction of the edge. + */ + Nedelec1stKindCubeLocalBasis() + { + std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0); + } + + /** \brief Constructor for a given edge orientation + */ + Nedelec1stKindCubeLocalBasis(std::bitset edgeOrientation) + : Nedelec1stKindCubeLocalBasis() + { + for (std::size_t i=0; i& out) const + { + static_assert(k==1, "Evaluating Nédélec shape functions is implemented only for first order."); + out.resize(size()); + + if (dim==2) + { + // First-order Nédélec shape functions on a square are of the form + // + // (a, b)^T + (c y, d x)^T, a, b, c, d \in R + // + // The following coefficients create the four basis vectors + // that are dual to the edge degrees of freedom: + // + // a[0] = 0 b[0] = 1 c[0] = 0 d[0] = -1 + // a[1] = 0 b[1] = 0 c[1] = 0 d[1] = 1 + // a[2] = 1 b[2] = 0 c[2] = 0 d[2] = -1 + // a[3] = 0 b[3] = 0 c[3] = 0 d[3] = 1 + + out[0] = { 0, D(1) - in[0]}; + out[1] = { 0, in[0]}; + out[2] = { D(1) - in[1], 0}; + out[3] = { in[1], 0}; + } + + if constexpr (dim==3) + { + // First-order Nédélec shape functions on a cube are of the form + // + // (e1 yz) + // a + b x + c y + d z + (e2 xz) , a, b, c, d \in R^3 and b[0]=c[1]=d[2]=0 + // (e3 xy) + // + // The following coefficients create the twelve basis vectors + // that are dual to the edge degrees of freedom: + // + // a[0] = { 0, 0, 1} b[0] = { 0, 0, -1} c[0] = { 0, 0, -1} d[0] = { 0, 0, 0} e[0] = { 0, 0, 1} + // a[1] = { 0, 0, 0} b[1] = { 0, 0, 1} c[1] = { 0, 0, 0} d[1] = { 0, 0, 0} e[1] = { 0, 0, -1} + // a[2] = { 0, 0, 0} b[2] = { 0, 0, 0} c[2] = { 0, 0, 1} d[2] = { 0, 0, 0} e[2] = { 0, 0, -1} + // a[3] = { 0, 0, 0} b[3] = { 0, 0, 0} c[3] = { 0, 0, 0} d[3] = { 0, 0, 0} e[3] = { 0, 0, 1} + // + // The following implementation uses these values, and simply + // skips all the zeros. + + for (std::size_t i=0; i& out) const + { + out.resize(size()); + if (dim==2) + { + for (std::size_t i=0; i& order, + const typename Traits::DomainType& in, + std::vector& out) const + { + auto totalOrder = std::accumulate(order.begin(), order.end(), 0); + if (totalOrder == 0) { + evaluateFunction(in, out); + } else if (totalOrder == 1) { + auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1)); + out.resize(size()); + + if (dim==2) + { + if (direction==0) + { + out[0] = { 0, -1}; + out[1] = { 0, 1}; + out[2] = { 0, 0}; + out[3] = { 0, 0}; + } + else + { + out[0] = { 0, 0}; + out[1] = { 0, 0}; + out[2] = { -1, 0}; + out[3] = { 1, 0}; + } + } + + if (dim==3) + { + switch (direction) + { + case 0: + out[0] = { 0, 0, -1 +in[1]}; + out[1] = { 0, 0, 1 -in[1]}; + out[2] = { 0, 0, -in[1]}; + out[3] = { 0, 0, in[1]}; + + out[4] = { 0, -1 +in[2], 0}; + out[5] = { 0, 1 -in[2], 0}; + out[8] = { 0, -in[2], 0}; + out[9] = { 0, in[2], 0}; + + out[6] = {0,0,0}; + out[7] = {0,0,0}; + out[10] = {0,0,0}; + out[11] = {0,0,0}; + break; + + case 1: + out[0] = { 0, 0, -1 + in[0]}; + out[1] = { 0, 0, - in[0]}; + out[2] = { 0, 0, 1 - in[0]}; + out[3] = { 0, 0, in[0]}; + + out[4] = {0,0,0}; + out[5] = {0,0,0}; + out[8] = {0,0,0}; + out[9] = {0,0,0}; + + out[6] = { -1 + in[2], 0, 0}; + out[7] = { 1 - in[2], 0, 0}; + out[10] = { - in[2], 0, 0}; + out[11] = { in[2], 0, 0}; + break; + + case 2: + out[0] = {0,0,0}; + out[1] = {0,0,0}; + out[2] = {0,0,0}; + out[3] = {0,0,0}; + + out[4] = { 0, -1 + in[0], 0}; + out[5] = { 0, - in[0], 0}; + out[8] = { 0, 1 - in[0], 0}; + out[9] = { 0, in[0], 0}; + + out[6] = { -1 + in[1], 0, 0}; + out[7] = { - in[1], 0, 0}; + out[10] = { 1 - in[1], 0, 0}; + out[11] = { in[1], 0, 0}; + break; + } + } + + for (std::size_t i=0; i edgeOrientation_; + }; + + + /** \brief Assignment of Nedelec degrees of freedom to subentities of the reference cube + * \tparam dim Dimension of the domain + * \tparam k Order of the Nedelec finite element + */ + template + class Nedelec1stKindCubeLocalCoefficients + { + public: + //! \brief Default constructor + Nedelec1stKindCubeLocalCoefficients () + : localKey_(size()) + { + static_assert(k==1, "Only first-order Nédélec local coefficients are implemented."); + // Assign all degrees of freedom to edges + // TODO: This is correct only for first-order Nédélec elements + for (std::size_t i=0; i localKey_; + }; + + /** \brief Project a given function into the Nedelec space + * + * \tparam LB The LocalBasis that spans the space + */ + template + class Nedelec1stKindCubeLocalInterpolation + { + static constexpr auto dim = LB::Traits::dimDomain; + static constexpr auto size = LB::size(); + + // Number of edges of the reference cube + constexpr static std::size_t numberOfEdges = power(2,dim-1)*dim; + + public: + + //! \brief Constructor with given set of edge orientations + Nedelec1stKindCubeLocalInterpolation (std::bitset s = 0) + { + auto refElement = Dune::referenceElement(GeometryTypes::cube(dim)); + + for (std::size_t i=0; iv1) + std::swap(v0,v1); + edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim); + edge_[i] *= (s[i]) ? -1.0 : 1.0; + } + } + + /** \brief Project a given function into the Nedelec space + * + * \param f The function to interpolate, must implement RangeField operator(DomainType) + * \param[out] out The coefficients of the projection + */ + template + void interpolate (const F& ff, std::vector& out) const + { + out.resize(size); + auto&& f = Impl::makeFunctionWithCallOperator(ff); + + for (std::size_t i=0; i m_; + // Edges of the reference cube + std::array edge_; + }; + +} + + + /** + * \brief Nédélec elements of the first kind for cube elements + * + * These elements have been described in : + * + * J.C. Nédélec, "Mixed finite elements in R^3", Numer.Math., 35(3):315-341,1980. + * DOI: http://dx.doi.org/10.1007/BF01396415 + * + * The order count starts at '1'. This is the counting used, e.g., by Nédélec himself, + * and by Kirby, Logg, Rognes, Terrel, "Common and unusual finite elements", + * https://doi.org/10.1007/978-3-642-23099-8_3 + * + * \note These shape functions are implemented for the reference cube only! + * The transformation to other cubes has to be done by the user. + * + * \ingroup Nedelec + * + * \tparam D Number type used for domain coordinates + * \tparam R Number type use for shape function values + * \tparam dim Dimension of the domain + * \tparam k Order of the Nedelec finite element (lowest is 1) + * + */ + template + class Nedelec1stKindCubeLocalFiniteElement + { + public: + using Traits = LocalFiniteElementTraits, + Impl::Nedelec1stKindCubeLocalCoefficients, + Impl::Nedelec1stKindCubeLocalInterpolation > >; + + static_assert(dim==2 || dim==3, "Nedelec elements are only implemented for 2d and 3d elements."); + static_assert(k==1, "Nedelec elements of the first kind are currently only implemented for order k==1."); + + /** \brief Default constructor + */ + Nedelec1stKindCubeLocalFiniteElement() = default; + + /** + * \brief Constructor with explicitly given edge orientations + * + * \param s Edge orientation indicator + */ + Nedelec1stKindCubeLocalFiniteElement (std::bitset s) : + basis_(s), + interpolation_(s) + {} + + const typename Traits::LocalBasisType& localBasis () const + { + return basis_; + } + + const typename Traits::LocalCoefficientsType& localCoefficients () const + { + return coefficients_; + } + + const typename Traits::LocalInterpolationType& localInterpolation () const + { + return interpolation_; + } + + static constexpr unsigned int size () + { + return Traits::LocalBasisType::size(); + } + + static constexpr GeometryType type () + { + return GeometryTypes::cube(dim); + } + + private: + typename Traits::LocalBasisType basis_; + typename Traits::LocalCoefficientsType coefficients_; + typename Traits::LocalInterpolationType interpolation_; + }; + +} + +#endif diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelec1stkindsimplex.hh dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelec1stkindsimplex.hh --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelec1stkindsimplex.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelec1stkindsimplex.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,457 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDSIMPLEX_HH +#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDSIMPLEX_HH + +#include + +#include +#include + +#include +#include + +#include +#include +#include // For deprecated makeFunctionWithCallOperator +#include + +namespace Dune +{ +namespace Impl +{ + /** \brief Nedelec shape functions of the first kind on reference simplices + * + * \note These shape functions are implemented for the reference simplices only! + * The transformation to other simplices has to be done by the user. + * + * \tparam D Number type used for domain coordinates + * \tparam R Number type used for function values + * \tparam dim Dimension of the reference element + * \tparam k Order of the Nedelec element + */ + template + class Nedelec1stKindSimplexLocalBasis + { + // Number of edges of the reference simplex + constexpr static std::size_t numberOfEdges = dim*(dim+1)/2; + + public: + using Traits = LocalBasisTraits, + R,dim,FieldVector, + FieldMatrix >; + + /** \brief Constructor with default edge orientation + * + * Default orientation means that all edges point from the vertex with lower index + * to the vertex with higher index. The tangential parts of all shape functions + * point in the direction of the edge. + */ + Nedelec1stKindSimplexLocalBasis() + { + std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0); + } + + /** \brief Constructor for a given edge orientation + */ + Nedelec1stKindSimplexLocalBasis(std::bitset edgeOrientation) + : Nedelec1stKindSimplexLocalBasis() + { + for (std::size_t i=0; i& out) const + { + static_assert(k==1, "Evaluating Nédélec shape functions is implemented only for first order."); + out.resize(size()); + + if (dim==2) + { + // First-order Nédélec shape functions on a triangle are of the form + // + // (a1, a2) + b(-x2, x1)^T, a_1, a_2, b \in R + out[0] = {D(1) - in[1], in[0]}; + out[1] = {in[1], -in[0]+D(1)}; + out[2] = {-in[1], in[0]}; + } + + if constexpr (dim==3) + { + // First-order Nédélec shape functions on a tetrahedron are of the form + // + // a + b \times x, a, b \in R^3 + // + // The following coefficients create the six basis vectors + // that are dual to the edge degrees of freedom: + // + // a[0] = { 1, 0, 0} b[0] = { 0, -1, 1} + // a[1] = { 0, 1, 0} b[1] = { 1, 0, -1} + // a[2] = { 0, 0, 0} b[2] = { 0, 0, 1} + // a[3] = { 0, 0, 1} b[3] = {-1, 1, 0} + // a[4] = { 0, 0, 0} b[4] = { 0, -1, 0} + // a[5] = { 0, 0, 0} b[5] = { 1, 0, 0} + // + // The following implementation uses these values, and simply + // skips all the zeros. + + out[0] = { 1 - in[1] - in[2], in[0] , in[0] }; + out[1] = { in[1] , 1 - in[0] - in[2], in[1]}; + out[2] = { - in[1] , in[0] , 0 }; + out[3] = { in[2], in[2], 1 - in[0] - in[1]}; + out[4] = { -in[2], 0 , in[0] }; + out[5] = { 0 , -in[2], in[1]}; + } + + for (std::size_t i=0; i& out) const + { + out.resize(size()); + if (dim==2) + { + out[0][0] = { 0, -1}; + out[0][1] = { 1, 0}; + + out[1][0] = { 0, 1}; + out[1][1] = {-1, 0}; + + out[2][0] = { 0, -1}; + out[2][1] = { 1, 0}; + } + if (dim==3) + { + out[0][0] = { 0,-1,-1}; + out[0][1] = { 1, 0, 0}; + out[0][2] = { 1, 0, 0}; + + out[1][0] = { 0, 1, 0}; + out[1][1] = {-1, 0, -1}; + out[1][2] = { 0, 1, 0}; + + out[2][0] = { 0, -1, 0}; + out[2][1] = { 1, 0, 0}; + out[2][2] = { 0, 0, 0}; + + out[3][0] = { 0, 0, 1}; + out[3][1] = { 0, 0, 1}; + out[3][2] = {-1, -1, 0}; + + out[4][0] = { 0, 0, -1}; + out[4][1] = { 0, 0, 0}; + out[4][2] = { 1, 0, 0}; + + out[5][0] = { 0, 0, 0}; + out[5][1] = { 0, 0, -1}; + out[5][2] = { 0, 1, 0}; + } + + for (std::size_t i=0; i& order, + const typename Traits::DomainType& in, + std::vector& out) const + { + auto totalOrder = std::accumulate(order.begin(), order.end(), 0); + if (totalOrder == 0) { + evaluateFunction(in, out); + } else if (totalOrder == 1) { + auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1)); + out.resize(size()); + + if (dim==2) + { + if (direction==0) + { + out[0] = {0, 1}; + out[1] = {0, -1}; + out[2] = {0, 1}; + } + else + { + out[0] = {-1, 0}; + out[1] = { 1, 0}; + out[2] = {-1, 0}; + } + } + + if (dim==3) + { + switch (direction) + { + case 0: + out[0] = { 0, 1, 1}; + out[1] = { 0,-1, 0}; + out[2] = { 0, 1, 0}; + out[3] = { 0, 0,-1}; + out[4] = { 0, 0, 1}; + out[5] = { 0, 0, 0}; + break; + + case 1: + out[0] = {-1, 0, 0}; + out[1] = { 1, 0, 1}; + out[2] = {-1, 0, 0}; + out[3] = { 0, 0,-1}; + out[4] = { 0, 0, 0}; + out[5] = { 0, 0, 1}; + break; + + case 2: + out[0] = {-1, 0, 0}; + out[1] = { 0,-1, 0}; + out[2] = { 0, 0, 0}; + out[3] = { 1, 1, 0}; + out[4] = {-1, 0, 0}; + out[5] = { 0,-1, 0}; + break; + } + } + + for (std::size_t i=0; i edgeOrientation_; + }; + + + /** \brief Assignment of Nedelec degrees of freedom to subentities of the reference simplex + * \tparam dim Dimension of the domain + * \tparam k Order of the Nedelec finite element + */ + template + class Nedelec1stKindSimplexLocalCoefficients + { + public: + //! \brief Default constructor + Nedelec1stKindSimplexLocalCoefficients () + : localKey_(size()) + { + static_assert(k==1, "Only first-order Nédélec local coefficients are implemented."); + // Assign all degrees of freedom to edges + // TODO: This is correct only for first-order Nédélec elements + for (std::size_t i=0; i localKey_; + }; + + /** \brief Project a given function into the Nedelec space + * + * \tparam LB The LocalBasis that spans the space + */ + template + class Nedelec1stKindSimplexLocalInterpolation + { + static constexpr auto dim = LB::Traits::dimDomain; + static constexpr auto size = LB::size(); + + // Number of edges of the reference simplex + constexpr static std::size_t numberOfEdges = dim*(dim+1)/2; + + public: + + //! \brief Constructor with given set of edge orientations + Nedelec1stKindSimplexLocalInterpolation (std::bitset s = 0) + { + auto refElement = Dune::referenceElement(GeometryTypes::simplex(dim)); + + for (std::size_t i=0; iv1) + std::swap(v0,v1); + edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim); + edge_[i] *= (s[i]) ? -1.0 : 1.0; + } + } + + /** \brief Project a given function into the Nedelec space + * + * \param f The function to interpolate, must implement RangeField operator(DomainType) + * \param[out] out The coefficients of the projection + */ + template + void interpolate (const F& ff, std::vector& out) const + { + out.resize(size); + auto&& f = Impl::makeFunctionWithCallOperator(ff); + + for (std::size_t i=0; i m_; + // Edges of the reference simplex + std::array edge_; + }; + +} + + + /** + * \brief Nédélec elements of the first kind for simplex elements + * + * These elements have been described in : + * + * J.C. Nédélec, "Mixed finite elements in R^3", Numer.Math., 35(3):315-341,1980. + * DOI: http://dx.doi.org/10.1007/BF01396415 + * + * The order count starts at '1'. This is the counting used, e.g., by Nédélec himself, + * and by Kirby, Logg, Rognes, Terrel, "Common and unusual finite elements", + * https://doi.org/10.1007/978-3-642-23099-8_3 + * + * \note These shape functions are implemented for the reference simplices only! + * The transformation to other simplices has to be done by the user. + * One way of doing that is using the implementation of the covariant + * Piola transform in globalvaluedlocalfinitelement.hh in dune-functions. + * + * \ingroup Nedelec + * + * \tparam D Number type used for domain coordinates + * \tparam R Number type use for shape function values + * \tparam dim Dimension of the domain + * \tparam k Order of the Nedelec finite element (lowest is 1) + * + */ + template + class Nedelec1stKindSimplexLocalFiniteElement + { + public: + using Traits = LocalFiniteElementTraits, + Impl::Nedelec1stKindSimplexLocalCoefficients, + Impl::Nedelec1stKindSimplexLocalInterpolation > >; + + static_assert(dim==2 || dim==3, "Nedelec elements are only implemented for 2d and 3d elements."); + static_assert(k==1, "Nedelec elements of the first kind are currently only implemented for order k==1."); + + /** \brief Default constructor + */ + Nedelec1stKindSimplexLocalFiniteElement() = default; + + /** + * \brief Constructor with explicitly given edge orientations + * + * \param s Edge orientation indicator + */ + Nedelec1stKindSimplexLocalFiniteElement (std::bitset s) : + basis_(s), + interpolation_(s) + {} + + const typename Traits::LocalBasisType& localBasis () const + { + return basis_; + } + + const typename Traits::LocalCoefficientsType& localCoefficients () const + { + return coefficients_; + } + + const typename Traits::LocalInterpolationType& localInterpolation () const + { + return interpolation_; + } + + static constexpr unsigned int size () + { + return Traits::LocalBasisType::size(); + } + + static constexpr GeometryType type () + { + return GeometryTypes::simplex(dim); + } + + private: + typename Traits::LocalBasisType basis_; + typename Traits::LocalCoefficientsType coefficients_; + typename Traits::LocalInterpolationType interpolation_; + }; + +} + +#endif diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/CMakeLists.txt dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,5 @@ +install(FILES + nedelecsimplexbasis.hh + nedelecsimplexinterpolation.hh + nedelecsimplexprebasis.hh + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/localfunctions/nedelec/nedelecsimplex) diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexbasis.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,42 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXBASIS_HH + +#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXBASIS_HH + +#include +#include + +#include +#include "nedelecsimplexinterpolation.hh" +#include "nedelecsimplexprebasis.hh" + +namespace Dune +{ + /* + * `NedelecPreBasisFactory` provides a basis for the Nedelec function space. + * `NedelecL2InterpolationFactory` provides the linear functionals. + * + * `Defaultbasisfactory::create` first builds the function space and the linear functionals. + * Then the constructor of `BasisMatrix` gets called. There the matrix + * + * \begin{equation} + * A_{i,j} := N_j(\phi_i) + * \end{equation} + * + * with linear functionals $N_j$ and basisfunctions $\phi_i$ gets assembled. + * Then the matrix gets inverted and is then used as a coefficent matrix for the standard monomial basis. + * + * For more details on the theory see the first chapter "Construction of Local Finite Element Spaces Using the Generic Reference Elements" + * of the book "Advances in Dune" by Dedner, Flemisch and Klöfkorn published in 2012. + */ + + template< unsigned int dim, class SF, class CF > + struct NedelecBasisFactory + : public DefaultBasisFactory< NedelecPreBasisFactory, + NedelecL2InterpolationFactory, + dim,dim,SF,CF > + {}; +} + +#endif // #ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXBASIS_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexinterpolation.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,687 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH +#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace Dune +{ + + // Internal Forward Declarations + // ----------------------------- + + template < unsigned int dim, class Field > + struct NedelecL2InterpolationFactory; + + + + // LocalCoefficientsContainer + // -------------------------- + + class LocalCoefficientsContainer + { + typedef LocalCoefficientsContainer This; + + public: + template + LocalCoefficientsContainer ( const Setter &setter ) + { + setter.setLocalKeys(localKey_); + } + + const LocalKey &localKey ( const unsigned int i ) const + { + assert( i < size() ); + return localKey_[ i ]; + } + + std::size_t size () const + { + return localKey_.size(); + } + + private: + std::vector< LocalKey > localKey_; + }; + + + + // NedelecCoefficientsFactory + // -------------------------------- + + template < unsigned int dim > + struct NedelecCoefficientsFactory + { + typedef std::size_t Key; + typedef const LocalCoefficientsContainer Object; + + template< GeometryType::Id geometryId > + static Object *create( const Key &key ) + { + typedef NedelecL2InterpolationFactory< dim, double > InterpolationFactory; + if( !supports< geometryId >( key ) ) + return nullptr; + typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key ); + Object *localKeys = new Object( *interpolation ); + InterpolationFactory::release( interpolation ); + return localKeys; + } + + template< GeometryType::Id geometryId > + static bool supports ( const Key &key ) + { + GeometryType gt = geometryId; + return gt.isTriangle() || gt.isTetrahedron() ; + } + static void release( Object *object ) { delete object; } + }; + + + + // NedelecL2InterpolationBuilder + // ------------------------ + + // L2 Interpolation requires: + // - for element + // - test basis + // - for each face (dynamic) + // - test basis + // - tangents + // - for each edge (dynamic) + // - test basis + // - tangent + template< unsigned int dim, class Field > + struct NedelecL2InterpolationBuilder + { + static const unsigned int dimension = dim; + + // for the dofs associated to the element + typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory; + typedef typename TestBasisFactory::Object TestBasis; + + // for the dofs associated to the faces + typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory; + typedef typename TestFaceBasisFactory::Object TestFaceBasis; + + // for the dofs associated to the edges + typedef OrthonormalBasisFactory< 1, Field > TestEdgeBasisFactory; + typedef typename TestEdgeBasisFactory::Object TestEdgeBasis; + + // the tangent of the edges + typedef FieldVector< Field, dimension > Tangent; + + // the normal and the tangents of the faces + typedef FieldVector< Field, dimension > Normal; + typedef std::array,dim-1> FaceTangents; + + NedelecL2InterpolationBuilder () = default; + + NedelecL2InterpolationBuilder ( const NedelecL2InterpolationBuilder & ) = delete; + NedelecL2InterpolationBuilder ( NedelecL2InterpolationBuilder && ) = delete; + + ~NedelecL2InterpolationBuilder () + { + TestBasisFactory::release( testBasis_ ); + for( FaceStructure &f : faceStructure_ ) + TestFaceBasisFactory::release( f.basis_ ); + for( EdgeStructure& e : edgeStructure_ ) + TestEdgeBasisFactory::release( e.basis_ ); + } + + unsigned int topologyId () const + { + return geometry_.id(); + } + + GeometryType type () const + { + return geometry_; + } + + std::size_t order () const + { + return order_; + } + + // number of faces + unsigned int faceSize () const + { + return numberOfFaces_; + } + + // number of edges + unsigned int edgeSize () const + { + return numberOfEdges_; + } + + // basis associated to the element + TestBasis *testBasis () const + { + return testBasis_; + } + + // basis associated to face f + TestFaceBasis *testFaceBasis ( unsigned int f ) const + { + assert( f < faceSize() ); + return faceStructure_[ f ].basis_; + } + + // basis associated to edge e + TestEdgeBasis *testEdgeBasis ( unsigned int e ) const + { + assert( e < edgeSize() ); + return edgeStructure_[ e ].basis_; + } + + const Tangent& edgeTangent ( unsigned int e ) const + { + assert( e < edgeSize() ); + return edgeStructure_[ e ].tangent_; + } + + const FaceTangents& faceTangents ( unsigned int f ) const + { + assert( f < faceSize() ); + return faceStructure_[ f ].faceTangents_; + } + + const Normal &normal ( unsigned int f ) const + { + assert( f < faceSize() ); + return faceStructure_[ f ].normal_; + } + + template< GeometryType::Id geometryId > + void build ( std::size_t order ) + { + constexpr GeometryType geometry = geometryId; + order_ = order; + geometry_ = geometry; + + /* + * The Nedelec parameter begins at 1. + * This is the numbering used by J.C. Nedelec himself. + * See "Mixed Finite Elements in \R^3" published in 1980. + * + * This construction is based on the construction of Raviart-Thomas elements. + * There the numbering starts at 0. + * Because of this we reduce the order internally by 1. + */ + order--; + + // if dimension == 2: order-1 on element + // if dimension == 3: order-2 on element + int requiredOrder = static_cast(dimension==3); + testBasis_ = (order > requiredOrder ? TestBasisFactory::template create< geometry >( order-1-requiredOrder ) : nullptr); + + const auto &refElement = ReferenceElements< Field, dimension >::general( type() ); + + numberOfFaces_ = refElement.size( 1 ); + faceStructure_.reserve( numberOfFaces_ ); + + // compute the basis, tangents and normals of each face + for (std::size_t i=0; i zero(0); + FaceTangents faceTangents; + faceTangents.fill(zero); + + // use the first dim-1 vertices of a face to compute the tangents + auto vertices = refElement.subEntities(i,1,dim).begin(); + auto vertex1 = *vertices; + for(int j=1; jvertex2) + faceTangents[j-1] *=-1; + + vertex1 = vertex2; + } + + /* For simplices or cubes of arbitrary dimension you could just use + * + * ``` + * GeometryType faceGeometry = Impl::getBase(geometry_); + * TestFaceBasis *faceBasis = ( dim == 3 && order > 0 ? TestFaceBasisFactory::template create< faceGeometry >( order-1 ) : nullptr); + * ``` + * + * For i.e. Prisms and Pyramids in 3d this does not work because they contain squares and triangles as faces. + * And depending on the dynamic face index a different face geometry is needed. + * + */ + TestFaceBasis *faceBasis = ( dim == 3 && order > 0 ? Impl::IfGeometryType< CreateFaceBasis, dimension-1 >::apply( refElement.type( i, 1 ), order-1 ) : nullptr); + faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal(i), faceTangents ); + } + assert( faceStructure_.size() == numberOfFaces_ ); + + numberOfEdges_ = refElement.size( dim-1 ); + edgeStructure_.reserve( numberOfEdges_ ); + + // compute the basis and tangent of each edge + for (std::size_t i=0; iv1) + std::swap(v0,v1); + auto tangent = std::move(refElement.position(v1,dim) - refElement.position(v0,dim)); + + TestEdgeBasis *edgeBasis = Impl::IfGeometryType< CreateEdgeBasis, 1 >::apply( refElement.type( i, dim-1 ), order ); + edgeStructure_.emplace_back( edgeBasis, tangent ); + } + assert( edgeStructure_.size() == numberOfEdges_ ); + } + + private: + + // helper struct for edges + struct EdgeStructure + { + EdgeStructure( TestEdgeBasis *teb, const Tangent &t ) + : basis_( teb ), tangent_( t ) + {} + + TestEdgeBasis *basis_; + const Dune::FieldVector< Field, dimension > tangent_; + }; + + template< GeometryType::Id edgeGeometryId > + struct CreateEdgeBasis + { + static TestEdgeBasis *apply ( std::size_t order ) { return TestEdgeBasisFactory::template create< edgeGeometryId >( order ); } + }; + + // helper struct for faces + struct FaceStructure + { + FaceStructure( TestFaceBasis *tfb, const Normal& normal, const FaceTangents& faceTangents ) + : basis_( tfb ), normal_(normal), faceTangents_( faceTangents ) + {} + + TestFaceBasis *basis_; + const Dune::FieldVector< Field, dimension > normal_; + const FaceTangents faceTangents_; + }; + + template< GeometryType::Id faceGeometryId > + struct CreateFaceBasis + { + static TestFaceBasis *apply ( std::size_t order ) { return TestFaceBasisFactory::template create< faceGeometryId >( order ); } + }; + + TestBasis *testBasis_ = nullptr; + std::vector< FaceStructure > faceStructure_; + unsigned int numberOfFaces_; + std::vector< EdgeStructure > edgeStructure_; + unsigned int numberOfEdges_; + GeometryType geometry_; + std::size_t order_; + }; + + + + // NedelecL2Interpolation + // ---------------------------- + + /** + * \class NedelecL2Interpolation + * \brief An L2-based interpolation for Nedelec + * + **/ + template< unsigned int dimension, class F> + class NedelecL2Interpolation + : public InterpolationHelper< F ,dimension > + { + typedef NedelecL2Interpolation< dimension, F > This; + typedef InterpolationHelper Base; + + public: + typedef F Field; + typedef NedelecL2InterpolationBuilder Builder; + typedef typename Builder::FaceTangents FaceTangents; + + NedelecL2Interpolation() + : order_(0), + size_(0) + {} + + template< class Function, class Vector > + auto interpolate ( const Function &function, Vector &coefficients ) const + -> std::enable_if_t< std::is_same< decltype(std::declval().resize(1) ),void >::value,void> + { + coefficients.resize(size()); + typename Base::template Helper func( function,coefficients ); + interpolate(func); + } + + template< class Basis, class Matrix > + auto interpolate ( const Basis &basis, Matrix &matrix ) const + -> std::enable_if_t< std::is_same< + decltype(std::declval().rowPtr(0)), typename Matrix::Field* >::value,void> + { + matrix.resize( size(), basis.size() ); + typename Base::template Helper func( basis,matrix ); + interpolate(func); + } + + std::size_t order() const + { + return order_; + } + std::size_t size() const + { + return size_; + } + + template + void build( std::size_t order ) + { + size_ = 0; + order_ = order; + builder_.template build(order_); + if (builder_.testBasis()) + size_ += dimension*builder_.testBasis()->size(); + + for ( unsigned int f=0; fsize(); + + for ( unsigned int e=0; esize(); + } + + void setLocalKeys(std::vector< LocalKey > &keys) const + { + keys.resize(size()); + unsigned int row = 0; + for (unsigned int e=0; esize(); ++i,++row) + keys[row] = LocalKey(e,dimension-1,i); + } + for (unsigned int f=0; fsize()*(dimension-1); ++i,++row) + keys[row] = LocalKey(f,1,i); + } + + if (builder_.testBasis()) + for (unsigned int i=0; isize()*dimension; ++i,++row) + keys[row] = LocalKey(0,0,i); + assert( row == size() ); + } + + protected: + template< class Func, class Container, bool type > + void interpolate ( typename Base::template Helper &func ) const + { + const Dune::GeometryType geoType( builder_.topologyId(), dimension ); + + std::vector testBasisVal; + + for (unsigned int i=0; i EdgeQuadrature; + typedef Dune::QuadratureRules EdgeQuadratureRules; + + const auto &refElement = Dune::ReferenceElements< Field, dimension >::general( geoType ); + + for (unsigned int e=0; esize()); + + const auto &geometry = refElement.template geometry< dimension-1 >( e ); + const Dune::GeometryType subGeoType( geometry.type().id(), 1 ); + const EdgeQuadrature &edgeQuad = EdgeQuadratureRules::rule( subGeoType, 2*order_+2 ); + + const unsigned int quadratureSize = edgeQuad.size(); + for( unsigned int qi = 0; qi < quadratureSize; ++qi ) + { + if (dimension>1) + builder_.testEdgeBasis(e)->template evaluate<0>(edgeQuad[qi].position(),testBasisVal); + else + testBasisVal[0] = 1.; + computeEdgeDofs(row, + testBasisVal, + func.evaluate( geometry.global( edgeQuad[qi].position() ) ), + builder_.edgeTangent(e), + edgeQuad[qi].weight(), + func); + } + + row += builder_.testEdgeBasis(e)->size(); + } + + // face dofs: + typedef Dune::QuadratureRule FaceQuadrature; + typedef Dune::QuadratureRules FaceQuadratureRules; + + for (unsigned int f=0; fsize()); + + const auto &geometry = refElement.template geometry< 1 >( f ); + const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 ); + const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 ); + + const unsigned int quadratureSize = faceQuad.size(); + for( unsigned int qi = 0; qi < quadratureSize; ++qi ) + { + if (dimension>1) + builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal); + else + testBasisVal[0] = 1.; + + computeFaceDofs( row, + testBasisVal, + func.evaluate( geometry.global( faceQuad[qi].position() ) ), + builder_.faceTangents(f), + builder_.normal(f), + faceQuad[qi].weight(), + func); + } + + row += builder_.testFaceBasis(f)->size()*(dimension-1); + } + } + + // element dofs + if (builder_.testBasis()) + { + testBasisVal.resize(builder_.testBasis()->size()); + + typedef Dune::QuadratureRule Quadrature; + typedef Dune::QuadratureRules QuadratureRules; + const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 ); + + const unsigned int quadratureSize = elemQuad.size(); + for( unsigned int qi = 0; qi < quadratureSize; ++qi ) + { + builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal); + computeInteriorDofs(row, + testBasisVal, + func.evaluate(elemQuad[qi].position()), + elemQuad[qi].weight(), + func ); + } + + row += builder_.testBasis()->size()*dimension; + } + assert(row==size()); + } + + private: + /** /brief evaluate functionals associated to the edges + * + * \param startRow row of matrix to start + * \param mVal value of the testBasis at a quadrature point on an edge + * \param nedVal value of the nedelecBasis at a quadrature point on an edge + * \param tangent the tangent of the edge + * \param weight quadrature weight + * \param matrix result gets written into matrix starting with row: row + */ + template + void computeEdgeDofs (unsigned int startRow, + const MVal &mVal, + const NedVal &nedVal, + const FieldVector &tangent, + const Field &weight, + Matrix &matrix) const + { + const unsigned int endRow = startRow+mVal.size(); + typename NedVal::const_iterator nedIter = nedVal.begin(); + for ( unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col) + { + Field cFactor = (*nedIter)*tangent; + typename MVal::const_iterator mIter = mVal.begin(); + for (unsigned int row = startRow; row!=endRow; ++mIter, ++row ) + matrix.add(row,col, (weight*cFactor)*(*mIter) ); + + assert( mIter == mVal.end() ); + } + } + + /** /brief evaluate functionals associated to the faces + * + * \param startRow row of matrix to start + * \param mVal value of the testBasis at a quadrature point on a face + * \param nedVal value of the nedelecBasis at a quadrature point on a face + * \param faceTangents the tangents of the face + * \param normal the normal of the face + * \param weight quadrature weight + * \param matrix result gets written into matrix starting with row: row + */ + template + void computeFaceDofs (unsigned int startRow, + const MVal &mVal, + const NedVal &nedVal, + const FaceTangents& faceTangents, + const FieldVector &normal, + const Field &weight, + Matrix &matrix) const + { + const unsigned int endRow = startRow+mVal.size()*(dimension-1); + typename NedVal::const_iterator nedIter = nedVal.begin(); + for ( unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col) + { + auto const& u=*nedIter; + auto const& n=normal; + FieldVector nedTimesNormal = { u[1]*n[2]-u[2]*n[1], + u[2]*n[0]-u[0]*n[2], + u[0]*n[1]-u[1]*n[0]}; + typename MVal::const_iterator mIter = mVal.begin(); + for (unsigned int row = startRow; row!=endRow; ++mIter) + { + for(int i=0; i + void computeInteriorDofs (unsigned int startRow, + const MVal &mVal, + const NedVal &nedVal, + Field weight, + Matrix &matrix) const + { + const unsigned int endRow = startRow+mVal.size()*dimension; + typename NedVal::const_iterator nedIter = nedVal.begin(); + for ( unsigned int col = 0; col < nedVal.size() ; ++nedIter,++col) + { + typename MVal::const_iterator mIter = mVal.begin(); + for (unsigned int row = startRow; row!=endRow; ++mIter,row+=dimension ) + for (unsigned int i=0; i + struct NedelecL2InterpolationFactory + { + typedef NedelecL2InterpolationBuilder Builder; + typedef const NedelecL2Interpolation Object; + typedef std::size_t Key; + typedef typename std::remove_const::type NonConstObject; + + template + static Object *create( const Key &key ) + { + if ( !supports(key) ) + return 0; + NonConstObject *interpol = new NonConstObject(); + interpol->template build(key); + return interpol; + } + + template + static bool supports( const Key &key ) + { + GeometryType gt = geometryId; + return gt.isTriangle() || gt.isTetrahedron() ; + } + static void release( Object *object ) { delete object; } + }; + +} // namespace Dune + +#endif // #ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXINTERPOLATION_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexprebasis.hh dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexprebasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexprebasis.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/nedelec/nedelecsimplex/nedelecsimplexprebasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,276 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_HH +#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_HH + +#include +#include + +#include + +#include + +namespace Dune +{ + template < GeometryType::Id geometryId, class Field > + struct NedelecVecMatrix; + + template + struct NedelecPreBasisFactory + { + typedef MonomialBasisProvider MBasisFactory; + typedef typename MBasisFactory::Object MBasis; + typedef StandardEvaluator EvalMBasis; + typedef PolynomialBasisWithMatrix > Basis; + + typedef const Basis Object; + typedef std::size_t Key; + + template + struct EvaluationBasisFactory + { + typedef MonomialBasisProvider Type; + }; + + template< GeometryType::Id geometryId > + static Object *create ( Key order ) + { + /* + * The nedelec parameter begins at 1. + * This is the numbering used by J.C. Nedelec himself. + * See "Mixed Finite Elements in \R^3" published in 1980. + * + * This construction is based on the construction of Raviart-Thomas elements. + * There the numbering starts at 0. + * Because of this we reduce the order internally by 1. + */ + order--; + NedelecVecMatrix vecMatrix(order); + MBasis *mbasis = MBasisFactory::template create(order+1); + std::remove_const_t* tmBasis = new std::remove_const_t(*mbasis); + tmBasis->fill(vecMatrix); + return tmBasis; + } + static void release( Object *object ) { delete object; } + }; + + template + struct NedelecVecMatrix + { + static constexpr GeometryType geometry = geometryId; + static const unsigned int dim = geometry.dim(); + typedef MultiIndex MI; + typedef MonomialBasis MIBasis; + NedelecVecMatrix(std::size_t order) + { + /* + * Construction of Nedelec elements see "Mixed Finite Elements in \R^3" by Nedelec, 1980. + * + * Let $\P_{n,k}$ be the space of polynomials in $n$ variables with degree $\leq k$. + * The space of Nedelec functions in $n$ dimensions with index $k$ is defined as + * + * \begin{equation*} + * Ned_k := (\P_{n,k-1})^n \oplus \{p \in (\P_{n,k})^n: =0 \} + * \end{equation*} + * with $x=(x,y)$ in two dimensions and $x=(x,y,z)$ in three dimensions. + * + * For $Ned_k$ holds + * \begin{equation*} + * (\P_{n,k-1})^n \subset Ned_k \subset (\P_{n,k})^n. + * \end{equation*} + * + * We construct $(\P_{n,k})^n$ and and only use the monomials contained in $Ned$. + * + */ + if( dim!=2 && dim!=3 || !geometry.isSimplex()) + DUNE_THROW(Dune::NotImplemented,"High order nedelec elements are only supported by simplices in 2d and 3d"); + + MIBasis basis(order+1); + FieldVector< MI, dim > x; + /* + * Init MultiIndices + * x[0]=(1,0,0) x + * x[1]=(0,1,0) y + * x[2]=(0,0,1) z + */ + for( unsigned int i = 0; i < dim; ++i ) + x[i].set(i,1); + std::vector< MI > val( basis.size() ); + + // val now contains all monomials in $n$ dimensions with degree $\leq order+1$ + basis.evaluate( x, val ); + + col_ = basis.size(); + + // get $\dim (\P_{n,order-1})$ + unsigned int notHomogen = 0; + if (order>0) + notHomogen = basis.sizes()[order-1]; + + // the number of basis functions for the set of homogeneous polynomials with degree $order$ + unsigned int homogen = basis.sizes()[order]-notHomogen; + + /* + * 2D: + * \begin{equation*} + * Ned_{order} = (\P_{order-1})^2 \oplus (-y,x)^T \widetilde \P_{order-1} + * \end{equation*} + * + * It gets more complicated in higher dimensions. + * + * 3D: + * \begin{equation*} + * Ned_{order} = (\P_{n,order-1})^3 \oplus (z,0,-x)^T \widetilde \P_{n,order-1} \oplus (-y,x,0)^T \widetilde \P_{n,order-1} \oplus (0,-z,y)^T \widetilde \P_{n-1,order-1} + * \end{equation*} + * + * Note the last term. The index $n-1$ is on purpose. + * Else i.e. k=2 + * + * (0,z,-y)^T x = (z,0,-x)^T y - (y,-x,0)^T z. + * + */ + + /* + * compute the number of rows for the coefficient matrix + * + * row_ = dim* \dim Ned_{order} + */ + if (dim == 2) + row_ = (notHomogen*dim+homogen*(dim+1))*dim; + else if (dim==3) + { + // get dim \P_{n-1,order-1} + int homogenTwoVariables = 0; + for( int w = notHomogen; w0; + /* + * Init 9 rows to zero. + * If the homogeneous monomial has a positive x-exponent (0,-z,y) gets skipped ( see example for the Nedelec space in 3D ) + * In this case only 6 rows get initialised. + */ + for (unsigned int r=0; r + void row( const unsigned int row, Vector &vec ) const + { + const unsigned int N = cols(); + assert( vec.size() == N ); + for (unsigned int i=0; i +#include + +#endif // #ifndef DUNE_LOCALFUNCTIONS_NEDELEC_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/orthonormal/orthonormalbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/orthonormal/orthonormalbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/orthonormal/orthonormalbasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/orthonormal/orthonormalbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -36,17 +36,17 @@ typedef unsigned int Key; typedef const Basis Object; - typedef typename Impl::SimplexTopology< dim >::type SimplexTopology; + static constexpr GeometryType SimplexGeometry = GeometryTypes::simplex(dim); - template< class Topology > + template< GeometryType::Id geometryId > static Object *create ( const unsigned int order ) { - const MonomialBasisType &monomialBasis = *MonomialBasisProviderType::template create< SimplexTopology >( order ); + const MonomialBasisType &monomialBasis = *MonomialBasisProviderType::template create< SimplexGeometry >( order ); static CoefficientMatrix _coeffs; if( _coeffs.size() <= monomialBasis.size() ) { - ONBCompute::ONBMatrix< Topology, ComputeField > matrix( order ); + ONBCompute::ONBMatrix< geometryId, ComputeField > matrix( order ); _coeffs.fill( matrix ); } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/orthonormal/orthonormalcompute.hh dune-localfunctions-2.8.0/dune/localfunctions/orthonormal/orthonormalcompute.hh --- dune-localfunctions-2.7.1/dune/localfunctions/orthonormal/orthonormalcompute.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/orthonormal/orthonormalcompute.hh 2021-08-31 14:03:38.000000000 +0000 @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -35,64 +36,62 @@ // Integral // -------- - template< class Topology > - struct Integral; - - template< class Base > - struct Integral< Dune::Impl::Pyramid< Base > > + template< Dune::GeometryType::Id geometryId > + struct Integral { - template< int dim, class scalar_t > - static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha, - scalar_t &p, scalar_t &q ) - { - const int dimension = Base::dimension+1; - int i = alpha.z( Base::dimension ); - int ord = Integral< Base >::compute( alpha, p, q ); - p *= factorial< scalar_t >( 1, i ); - q *= factorial< scalar_t >( dimension + ord, dimension + ord + i ); - return ord + i; - } - }; + static constexpr Dune::GeometryType geometry = geometryId; + static constexpr int dimension = geometry.dim(); - template< class Base > - struct Integral< Dune::Impl::Prism< Base > > - { template< int dim, class scalar_t > static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q ) { - int i = alpha.z( Base::dimension ); - int ord = Integral< Base >::compute( alpha, p, q ); - //Integral< Base >::compute( alpha, p, q ); - //p *= scalar_t( 1 ); - q *= scalar_t( i+1 ); - return ord + i; + return compute(alpha, p, q, std::make_integer_sequence{}); } - }; - template<> - struct Integral< Dune::Impl::Point > - { - template< int dim, class scalar_t > + template< int dim, class scalar_t , int ...ints> static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha, - scalar_t &p, scalar_t &q ) + scalar_t &p, scalar_t &q, std::integer_sequence intS) { p = scalar_t( 1 ); q = scalar_t( 1 ); - return 0; + + int ord = 0; + ((computeIntegral(alpha,p,q,ord)),...); + + return ord; + } + + template< int step, int dim, class scalar_t > + static void computeIntegral ( const Dune::MultiIndex< dim, scalar_t > &alpha, + scalar_t &p, scalar_t &q, int& ord) + { + int i = alpha.z( step ); + + if constexpr ( geometry.isPrismatic(step)) + { + //p *= scalar_t( 1 ); + q *= scalar_t( i+1 ); + } + else + { + p *= factorial< scalar_t >( 1, i ); + q *= factorial< scalar_t >( step+1 + ord, step+1 + ord + i ); + } + ord +=i; } - }; + }; // ONBMatrix // --------- - template< class Topology, class scalar_t > + template< Dune::GeometryType::Id geometryId, class scalar_t > class ONBMatrix : public Dune::LFEMatrix< scalar_t > { - typedef ONBMatrix< Topology, scalar_t > This; + typedef ONBMatrix< geometryId, scalar_t > This; typedef Dune::LFEMatrix< scalar_t > Base; public: @@ -102,7 +101,8 @@ explicit ONBMatrix ( unsigned int order ) { // get all multiindecies for monomial basis - const unsigned int dim = Topology::dimension; + constexpr Dune::GeometryType geometry = geometryId; + constexpr unsigned int dim = geometry.dim(); typedef Dune::MultiIndex< dim, scalar_t > MI; Dune::StandardMonomialBasis< dim, MI > basis( order ); const std::size_t size = basis.size(); @@ -123,7 +123,7 @@ { for( std::size_t j = 0; j < size; ++j ) { - Integral< Topology >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q ); + Integral< geometryId >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q ); S( i, j ) = p; S( i, j ) /= q; } @@ -167,6 +167,7 @@ void gramSchmidt () { + using std::sqrt; // setup identity const std::size_t N = Base::rows(); for( std::size_t i = 0; i < N; ++i ) diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/CMakeLists.txt dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/CMakeLists.txt 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -1,4 +1,5 @@ add_subdirectory(raviartthomas02d) +add_subdirectory(raviartthomas03d) add_subdirectory(raviartthomas0cube2d) add_subdirectory(raviartthomas0cube3d) add_subdirectory(raviartthomas12d) @@ -11,6 +12,7 @@ install(FILES raviartthomas02d.hh + raviartthomas03d.hh raviartthomas0cube2d.hh raviartthomas0cube3d.hh raviartthomas12d.hh diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas02d/raviartthomas02dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas02d/raviartthomas02dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas02d/raviartthomas02dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas02d/raviartthomas02dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -19,6 +19,7 @@ //! \brief Constructor with given set of edge orientations RT02DLocalInterpolation (std::bitset<3> s = 0) { + using std::sqrt; for (std::size_t i=0; i + +#include + +#include + +namespace Dune +{ + /**@ingroup LocalBasisImplementation + \brief Lowest order Raviart-Thomas shape functions on the reference tetrahedron. + + \tparam D Type to represent the field in the domain. + \tparam R Type to represent the field in the range. + + \nosubgrouping + */ + template + class RT03DLocalBasis + { + public: + typedef LocalBasisTraits,R,3,Dune::FieldVector, + Dune::FieldMatrix > Traits; + + //! \brief Make set number s, where 0 <= s < 16 + RT03DLocalBasis (std::bitset<4> s = 0) + { + for (int i=0; i<4; i++) + sign_[i] = s[i] ? -1.0 : 1.0; + } + + //! \brief number of shape functions + unsigned int size () const + { + return 4; + } + + //! \brief Evaluate all shape functions + inline void evaluateFunction (const typename Traits::DomainType& in, + std::vector& out) const + { + out.resize(4); + auto c = std::sqrt(2.0); + out[0] = {sign_[0]*c* in[0], sign_[0]*c* in[1], sign_[0]*c*(in[2]-D(1))}; + out[1] = {sign_[1]*c* in[0], sign_[1]*c*(in[1]-D(1)), sign_[1]*c* in[2] }; + out[2] = {sign_[2]*c*(in[0]-D(1)), sign_[2]*c* in[1], sign_[2]*c* in[2] }; + out[3] = {sign_[3]*c* in[0], sign_[3]*c* in[1], sign_[3]*c* in[2] }; + } + + //! \brief Evaluate Jacobian of all shape functions + inline void + evaluateJacobian (const typename Traits::DomainType& in, // position + std::vector& out) const // return value + { + out.resize(4); + for (int i=0; i<4; i++) + { + auto c = std::sqrt(2.0); + out[i][0] = {c*sign_[i], 0, 0}; + out[i][1] = { 0,c*sign_[i], 0}; + out[i][2] = { 0, 0,c*sign_[i]}; + } + } + + //! \brief Evaluate partial derivatives of all shape functions + void partial (const std::array& order, + const typename Traits::DomainType& in, // position + std::vector& out) const // return value + { + auto totalOrder = std::accumulate(order.begin(), order.end(), 0); + if (totalOrder == 0) { + evaluateFunction(in, out); + } else if (totalOrder == 1) { + auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1)); + out.resize(size()); + + for (int i=0; i sign_; + }; +} +#endif diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalcoefficients.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalcoefficients.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalcoefficients.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalcoefficients.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,49 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMAS03D_RAVIARTTHOMAS03DLOCALCOEFFICIENTS_HH +#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMAS03D_RAVIARTTHOMAS03DLOCALCOEFFICIENTS_HH + +#include +#include +#include + +#include + +namespace Dune +{ + + /**@ingroup LocalLayoutImplementation + \brief Layout map for RT0 elements + + \nosubgrouping + \implements Dune::LocalCoefficientsVirtualImp + */ + class RT03DLocalCoefficients + { + public: + //! \brief Standard constructor + RT03DLocalCoefficients () : li(4) + { + for (std::size_t i=0; i<4; i++) + li[i] = LocalKey(i,1,0); + } + + //! number of coefficients + std::size_t size () const + { + return 4; + } + + //! get i'th index + const LocalKey& localKey (std::size_t i) const + { + return li[i]; + } + + private: + std::vector li; + }; + +} + +#endif diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalinterpolation.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas03d/raviartthomas03dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,67 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMAS03D_RAVIARTTHOMAS03DLOCALINTERPOLATION_HH +#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMAS03D_RAVIARTTHOMAS03DLOCALINTERPOLATION_HH + +#include +#include +#include +#include +#include + +namespace Dune +{ + template + class RT03DLocalInterpolation + { + public: + + //! \brief Constructor with given set of face orientations + RT03DLocalInterpolation (std::bitset<4> s = 0) + { + using std::sqrt; + for (std::size_t i=0; i + void interpolate (const F& ff, std::vector& out) const + { + // f gives v*outer normal at a point on the face! + auto&& f = Impl::makeFunctionWithCallOperator(ff); + + out.resize(4); + + for (int i=0; i<4; i++) + { + auto y = f(m_[i]); + out[i] = (y[0]*n_[i][0]+y[1]*n_[i][1]+y[2]*n_[i][2])*sign_[i]/c_[i]; + } + } + + private: + // Face orientations + std::array sign_; + // Face midpoints of the reference tetrahedron + std::array m_; + // Unit outer normals of the reference tetrahedron + std::array n_; + // Inverse triangle face area + std::array c_; + }; +} + +#endif diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas03d.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas03d.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas03d.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas03d.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,79 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMAS03D_HH +#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS_RAVIARTTHOMAS03D_HH + +#include + +#include +#include "raviartthomas03d/raviartthomas03dlocalbasis.hh" +#include "raviartthomas03d/raviartthomas03dlocalcoefficients.hh" +#include "raviartthomas03d/raviartthomas03dlocalinterpolation.hh" + +namespace Dune +{ + + /** + * \brief Zero order Raviart-Thomas shape functions on tetrahedra. + * + * \ingroup RaviartThomas + * + * \tparam D Type to represent the field in the domain. + * \tparam R Type to represent the field in the range. + */ + template + class + RT03DLocalFiniteElement + { + public: + typedef LocalFiniteElementTraits,RT03DLocalCoefficients, + RT03DLocalInterpolation > > Traits; + + //! \brief Standard constructor + RT03DLocalFiniteElement () + {} + + /** + * \brief Constructor with explicitly given face orientations + * + * \param s Face orientation indicator + */ + RT03DLocalFiniteElement (std::bitset<4> s) : + basis(s), + interpolation(s) + {} + + const typename Traits::LocalBasisType& localBasis () const + { + return basis; + } + + const typename Traits::LocalCoefficientsType& localCoefficients () const + { + return coefficients; + } + + const typename Traits::LocalInterpolationType& localInterpolation () const + { + return interpolation; + } + + unsigned int size () const + { + return 4; + } + + static constexpr GeometryType type () + { + return GeometryTypes::tetrahedron; + } + + private: + RT03DLocalBasis basis; + RT03DLocalCoefficients coefficients; + RT03DLocalInterpolation > interpolation; + }; + +} + +#endif diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalbasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -34,21 +34,10 @@ * * \param s Edge orientation indicator */ - RT12DLocalBasis (unsigned int s = 0) + RT12DLocalBasis (std::bitset<3> s = 0) { - sign0 = sign1 = sign2 = 1.0; - if (s & 1) - { - sign0 = -1.0; - } - if (s & 2) - { - sign1 = -1.0; - } - if (s & 4) - { - sign2 = -1.0; - } + for (size_t i=0; i<3; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; } //! \brief number of shape functions @@ -67,12 +56,12 @@ std::vector& out) const { out.resize(8); - out[0][0] = sign0*(in[0] - 4.0*in[0]*in[1]); - out[0][1] = sign0*(-1.0 + 5.0*in[1] - 4.0*in[1]*in[1]); - out[1][0] = sign1*(-1.0 + 5.0*in[0] - 4.0*in[0]*in[0]); - out[1][1] = sign1*(in[1] - 4.0*in[0]*in[1]); - out[2][0] = sign2*(-3.0*in[0] + 4.0*in[0]*in[0] + 4.0*in[1]*in[0]); - out[2][1] = sign2*(-3.0*in[1] + 4.0*in[0]*in[1] + 4.0*in[1]*in[1]); + out[0][0] = sign_[0]*(in[0] - 4.0*in[0]*in[1]); + out[0][1] = sign_[0]*(-1.0 + 5.0*in[1] - 4.0*in[1]*in[1]); + out[1][0] = sign_[1]*(-1.0 + 5.0*in[0] - 4.0*in[0]*in[0]); + out[1][1] = sign_[1]*(in[1] - 4.0*in[0]*in[1]); + out[2][0] = sign_[2]*(-3.0*in[0] + 4.0*in[0]*in[0] + 4.0*in[1]*in[0]); + out[2][1] = sign_[2]*(-3.0*in[1] + 4.0*in[0]*in[1] + 4.0*in[1]*in[1]); out[3][0] = -5.0*in[0] + 8.0*in[0]*in[0] + 4.0*in[1]*in[0]; out[3][1] = 3.0 - 6.0*in[0] - 7.0*in[1] + 8.0*in[0]*in[1] + 4.0*in[1]*in[1]; out[4][0] = -3.0 + 7.0*in[0] + 6.0*in[1] - 4.0*in[0]*in[0] - 8.0*in[1]*in[0]; @@ -96,20 +85,20 @@ { out.resize(8); - out[0][0][0] = sign0*(1.0 - 4.0*in[1]); - out[0][0][1] = sign0*(-4.0*in[0]); + out[0][0][0] = sign_[0]*(1.0 - 4.0*in[1]); + out[0][0][1] = sign_[0]*(-4.0*in[0]); out[0][1][0] = 0.0; - out[0][1][1] = sign0*(5.0 - 8.0*in[1]); + out[0][1][1] = sign_[0]*(5.0 - 8.0*in[1]); - out[1][0][0] = sign1*(5.0 - 8.0*in[0]); + out[1][0][0] = sign_[1]*(5.0 - 8.0*in[0]); out[1][0][1] = 0.0; - out[1][1][0] = sign1*(-4.0*in[1]); - out[1][1][1] = sign1*(1.0 - 4.0*in[0]); + out[1][1][0] = sign_[1]*(-4.0*in[1]); + out[1][1][1] = sign_[1]*(1.0 - 4.0*in[0]); - out[2][0][0] = sign2*(-3.0 + 8.0*in[0] + 4.0*in[1]); - out[2][0][1] = sign2*(4.0*in[0]); - out[2][1][0] = sign2*(4.0*in[1]); - out[2][1][1] = sign2*(-3.0 + 4.0*in[0] + 8.0*in[1]); + out[2][0][0] = sign_[2]*(-3.0 + 8.0*in[0] + 4.0*in[1]); + out[2][0][1] = sign_[2]*(4.0*in[0]); + out[2][1][0] = sign_[2]*(4.0*in[1]); + out[2][1][1] = sign_[2]*(-3.0 + 4.0*in[0] + 8.0*in[1]); out[3][0][0] = -5.0 + 16.0*in[0] + 4.0*in[1]; out[3][0][1] = 4.0*in[0]; @@ -151,12 +140,12 @@ switch (direction) { case 0: - out[0][0] = sign0*(1.0 - 4.0*in[1]); + out[0][0] = sign_[0]*(1.0 - 4.0*in[1]); out[0][1] = 0.0; - out[1][0] = sign1*(5.0 - 8.0*in[0]); - out[1][1] = sign1*(-4.0*in[1]); - out[2][0] = sign2*(-3.0 + 8.0*in[0] + 4.0*in[1]); - out[2][1] = sign2*(4.0*in[1]); + out[1][0] = sign_[1]*(5.0 - 8.0*in[0]); + out[1][1] = sign_[1]*(-4.0*in[1]); + out[2][0] = sign_[2]*(-3.0 + 8.0*in[0] + 4.0*in[1]); + out[2][1] = sign_[2]*(4.0*in[1]); out[3][0] = -5.0 + 16.0*in[0] + 4.0*in[1]; out[3][1] = -6.0 + 8.0*in[1]; out[4][0] = 7.0 - 8.0*in[0] - 8.0*in[1]; @@ -169,12 +158,12 @@ out[7][1] = -8.0*in[1]; break; case 1: - out[2][1] = sign2*(-3.0 + 4.0*in[0] + 8.0*in[1]); - out[2][0] = sign2*(4.0*in[0]); - out[1][1] = sign1*(1.0 - 4.0*in[0]); + out[2][1] = sign_[2]*(-3.0 + 4.0*in[0] + 8.0*in[1]); + out[2][0] = sign_[2]*(4.0*in[0]); + out[1][1] = sign_[1]*(1.0 - 4.0*in[0]); out[1][0] = 0.0; - out[0][0] = sign0*(-4.0*in[0]); - out[0][1] = sign0*(5.0 - 8.0*in[1]); + out[0][0] = sign_[0]*(-4.0*in[0]); + out[0][1] = sign_[0]*(5.0 - 8.0*in[1]); out[3][0] = 4.0*in[0]; out[3][1] = -7.0 + 8.0*in[0] + 8.0*in[1]; out[4][0] = 6.0 - 8.0*in[0]; @@ -201,7 +190,7 @@ } private: - R sign0, sign1, sign2; + std::array sign_; }; } #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALBASIS_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas12d/raviartthomas12dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -30,30 +30,19 @@ * * \param s Edge orientation indicator */ - RT12DLocalInterpolation (unsigned int s = 0) + RT12DLocalInterpolation (std::bitset<3> s = 0) { - sign0 = sign1 = sign2 = 1.0; - if (s & 1) - { - sign0 = -1.0; - } - if (s & 2) - { - sign1 = -1.0; - } - if (s & 4) - { - sign2 = -1.0; - } - n0[0] = 0.0; - n0[1] = -1.0; - n1[0] = -1.0; - n1[1] = 0.0; - n2[0] = 1.0/sqrt(2.0); - n2[1] = 1.0/sqrt(2.0); - c0 = 0.5*n0[0] - 1.0*n0[1]; - c1 = -1.0*n1[0] + 0.5*n1[1]; - c2 = 0.5*n2[0] + 0.5*n2[1]; + using std::sqrt; + for (size_t i=0; i<3; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; + + n_[0] = { 0.0, -1.0}; + n_[1] = {-1.0, 0.0}; + n_[2] = { 1.0/sqrt(2.0), 1.0/sqrt(2.0)}; + + c_ = { 0.5*n_[0][0] - 1.0*n_[0][1], + -1.0*n_[1][0] + 0.5*n_[1][1], + 0.5*n_[2][0] + 0.5*n_[2][1]}; } /** @@ -77,51 +66,50 @@ fill(out.begin(), out.end(), 0.0); const int qOrder1 = 4; - const Dune::QuadratureRule& rule1 = Dune::QuadratureRules::rule(Dune::GeometryTypes::simplex(1), qOrder1); + const auto& rule1 = Dune::QuadratureRules::rule(Dune::GeometryTypes::simplex(1), qOrder1); - for (typename Dune::QuadratureRule::const_iterator it = rule1.begin(); - it != rule1.end(); ++it) + for (auto&& qp : rule1) { - Scalar qPos = it->position(); + Scalar qPos = qp.position(); typename LB::Traits::DomainType localPos; - localPos[0] = qPos; - localPos[1] = 0.0; + localPos = {qPos, 0.0}; auto y = f(localPos); - out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0/c0; - out[3] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight()/c0; + out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]/c_[0]; + out[3] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight()/c_[0]; - localPos[0] = 0.0; - localPos[1] = qPos; + localPos = {0.0, qPos}; y = f(localPos); - out[1] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1/c1; - out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight()/c1; + out[1] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]/c_[1]; + out[4] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight()/c_[1]; - localPos[0] = 1.0 - qPos; - localPos[1] = qPos; + localPos = {1.0 - qPos, qPos}; y = f(localPos); - out[2] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2/c2; - out[5] += (y[0]*n2[0] + y[1]*n2[1])*(2.0*qPos - 1.0)*it->weight()/c2; + out[2] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]/c_[2]; + out[5] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(2.0*qPos - 1.0)*qp.weight()/c_[2]; } const int qOrder2 = 8; - const Dune::QuadratureRule& rule2 = Dune::QuadratureRules::rule(Dune::GeometryTypes::simplex(2), qOrder2); + const auto& rule2 = Dune::QuadratureRules::rule(Dune::GeometryTypes::simplex(2), qOrder2); - for (typename Dune::QuadratureRule::const_iterator it = rule2.begin(); - it != rule2.end(); ++it) + for (auto&& qp : rule2) { - Dune::FieldVector qPos = it->position(); + auto qPos = qp.position(); auto y = f(qPos); - out[6] += y[0]*it->weight(); - out[7] += y[1]*it->weight(); + out[6] += y[0]*qp.weight(); + out[7] += y[1]*qp.weight(); } } private: - typename LB::Traits::RangeFieldType sign0,sign1,sign2; - typename LB::Traits::DomainType n0,n1,n2; - typename LB::Traits::RangeFieldType c0,c1,c2; + // Edge orientations + std::array sign_; + + // Edge normals + std::array n_; + + std::array c_; }; } #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalbasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -34,25 +34,10 @@ * * \param s Edge orientation indicator */ - RT1Cube2DLocalBasis (unsigned int s = 0) + RT1Cube2DLocalBasis (std::bitset<4> s = 0) { - sign0 = sign1 = sign2 = sign3 = 1.0; - if (s & 1) - { - sign0 = -1.0; - } - if (s & 2) - { - sign1 = -1.0; - } - if (s & 4) - { - sign2 = -1.0; - } - if (s & 8) - { - sign3 = -1.0; - } + for (size_t i=0; i<4; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; } //! \brief number of shape functions @@ -72,20 +57,20 @@ { out.resize(12); - out[0][0] = sign0*(-1.0 + 4.0*in[0]-3*in[0]*in[0]); + out[0][0] = sign_[0]*(-1.0 + 4.0*in[0]-3*in[0]*in[0]); out[0][1] = 0.0; out[1][0] = 3.0 - 12.0*in[0] - 6.0*in[1] + 24.0*in[0]*in[1]+9*in[0]*in[0] - 18.0*in[0]*in[0]*in[1]; out[1][1] = 0.0; - out[2][0] = sign1*(-2.0*in[0] + 3.0*in[0]*in[0]); + out[2][0] = sign_[1]*(-2.0*in[0] + 3.0*in[0]*in[0]); out[2][1] = 0.0; out[3][0] = -6.0*in[0] + 12.0*in[0]*in[1] + 9.0*in[0]*in[0] - 18.0*in[0]*in[0]*in[1]; out[3][1] = 0.0; out[4][0] = 0.0; - out[4][1] = sign2*(-1.0 + 4.0*in[1] - 3.0*in[1]*in[1]); + out[4][1] = sign_[2]*(-1.0 + 4.0*in[1] - 3.0*in[1]*in[1]); out[5][0] = 0.0; out[5][1] = -3.0 + 6.0*in[0] + 12.0*in[1] - 24.0*in[0]*in[1] - 9.0*in[1]*in[1] + 18.0*in[0]*in[1]*in[1]; out[6][0] = 0.0; - out[6][1] = sign3*(-2.0*in[1] + 3.0*in[1]*in[1]); + out[6][1] = sign_[3]*(-2.0*in[1] + 3.0*in[1]*in[1]); out[7][0] = 0.0; out[7][1] = 6.0*in[1] - 12.0*in[0]*in[1] - 9.0*in[1]*in[1] + 18.0*in[0]*in[1]*in[1]; out[8][0] = 24.0*in[0] - 36.0*in[0]*in[1] - 24.0*in[0]*in[0] + 36.0*in[0]*in[0]*in[1]; @@ -109,7 +94,7 @@ { out.resize(12); - out[0][0][0] = sign0*(4.0 - 6.0*in[0]); + out[0][0][0] = sign_[0]*(4.0 - 6.0*in[0]); out[0][0][1] = 0.0; out[0][1][0] = 0.0; out[0][1][1] = 0.0; @@ -119,7 +104,7 @@ out[1][1][0] = 0.0; out[1][1][1] = 0.0; - out[2][0][0] = sign1*(-2.0 + 6.0*in[0]); + out[2][0][0] = sign_[1]*(-2.0 + 6.0*in[0]); out[2][0][1] = 0.0; out[2][1][0] = 0.0; out[2][1][1] = 0.0; @@ -132,7 +117,7 @@ out[4][0][0] = 0.0; out[4][0][1] = 0.0; out[4][1][0] = 0.0; - out[4][1][1] = sign2*(4.0 - 6.0*in[1]); + out[4][1][1] = sign_[2]*(4.0 - 6.0*in[1]); out[5][0][0] = 0.0; out[5][0][1] = 0.0; @@ -142,7 +127,7 @@ out[6][0][0] = 0.0; out[6][0][1] = 0.0; out[6][1][0] = 0.0; - out[6][1][1] = sign3*(-2.0 + 6.0*in[1]); + out[6][1][1] = sign_[3]*(-2.0 + 6.0*in[1]); out[7][0][0] = 0.0; out[7][0][1] = 0.0; @@ -190,7 +175,7 @@ } private: - R sign0, sign1, sign2, sign3; + std::array sign_; }; } #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE2D_LOCALBASIS_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas1cube2d/raviartthomas1cube2dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -30,34 +30,15 @@ * * \param s Edge orientation indicator */ - RT1Cube2DLocalInterpolation (unsigned int s = 0) + RT1Cube2DLocalInterpolation (std::bitset<4> s = 0) { - sign0 = sign1 = sign2 = sign3 = 1.0; - if (s & 1) - { - sign0 = -1.0; - } - if (s & 2) - { - sign1 = -1.0; - } - if (s & 4) - { - sign2 = -1.0; - } - if (s & 8) - { - sign3 = -1.0; - } + for (size_t i=0; i<4; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; - n0[0] = -1.0; - n0[1] = 0.0; - n1[0] = 1.0; - n1[1] = 0.0; - n2[0] = 0.0; - n2[1] = -1.0; - n3[0] = 0.0; - n3[1] = 1.0; + n_[0] = {-1.0, 0.0}; + n_[1] = { 1.0, 0.0}; + n_[2] = { 0.0, -1.0}; + n_[3] = { 0.0, 1.0}; } /** @@ -81,56 +62,53 @@ fill(out.begin(), out.end(), 0.0); const int qOrder = 3; - const QuadratureRule& rule1 = QuadratureRules::rule(GeometryTypes::cube(1), qOrder); + const auto& rule1 = QuadratureRules::rule(GeometryTypes::cube(1), qOrder); - for (typename QuadratureRule::const_iterator it = rule1.begin(); - it != rule1.end(); ++it) + for (auto&& qp : rule1) { - Scalar qPos = it->position(); - typename LB::Traits::DomainType localPos; + Scalar qPos = qp.position(); + typename LB::Traits::DomainType localPos = {0.0, qPos}; - localPos[0] = 0.0; - localPos[1] = qPos; auto y = f(localPos); - out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0; - out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight(); + out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]; + out[1] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight(); - localPos[0] = 1.0; - localPos[1] = qPos; + localPos = {1.0, qPos}; y = f(localPos); - out[2] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1; - out[3] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight(); + out[2] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]; + out[3] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight(); - localPos[0] = qPos; - localPos[1] = 0.0; + localPos = {qPos, 0.0}; y = f(localPos); - out[4] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2; - out[5] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight(); + out[4] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]; + out[5] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(1.0 - 2.0*qPos)*qp.weight(); - localPos[0] = qPos; - localPos[1] = 1.0; + localPos = {qPos, 1.0}; y = f(localPos); - out[6] += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3; - out[7] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight(); + out[6] += (y[0]*n_[3][0] + y[1]*n_[3][1])*qp.weight()*sign_[3]; + out[7] += (y[0]*n_[3][0] + y[1]*n_[3][1])*(2.0*qPos - 1.0)*qp.weight(); } - const QuadratureRule& rule2 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); + const auto& rule2 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); - for (typename QuadratureRule::const_iterator it=rule2.begin(); it!=rule2.end(); ++it) + for (auto&& qp : rule2) { - FieldVector qPos = it->position(); + auto qPos = qp.position(); auto y = f(qPos); - out[8] += y[0]*it->weight(); - out[9] += y[1]*it->weight(); - out[10] += y[0]*qPos[1]*it->weight(); - out[11] += y[1]*qPos[0]*it->weight(); + out[8] += y[0]*qp.weight(); + out[9] += y[1]*qp.weight(); + out[10] += y[0]*qPos[1]*qp.weight(); + out[11] += y[1]*qPos[0]*qp.weight(); } } private: - typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3; - typename LB::Traits::DomainType n0, n1, n2, n3; + // Edge orientations + std::array sign_; + + // Edge normals + std::array n_; }; } #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE2D_LOCALINTERPOLATION_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas1cube3d/raviartthomas1cube3dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas1cube3d/raviartthomas1cube3dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas1cube3d/raviartthomas1cube3dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas1cube3d/raviartthomas1cube3dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -29,52 +29,17 @@ * * \param s Edge orientation indicator */ - RT1Cube3DLocalInterpolation (unsigned int s = 0) + RT1Cube3DLocalInterpolation (std::bitset<6> s = 0) { - sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0; - if (s & 1) - { - sign0 = -1.0; - } - if (s & 2) - { - sign1 = -1.0; - } - if (s & 4) - { - sign2 = -1.0; - } - if (s & 8) - { - sign3 = -1.0; - } - if (s & 16) - { - sign4 = -1.0; - } - if (s & 32) - { - sign5 = -1.0; - } + for (size_t i=0; i<6; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; - n0[0] = -1.0; - n0[1] = 0.0; - n0[2] = 0.0; - n1[0] = 1.0; - n1[1] = 0.0; - n1[2] = 0.0; - n2[0] = 0.0; - n2[1] = -1.0; - n2[2] = 0.0; - n3[0] = 0.0; - n3[1] = 1.0; - n3[2] = 0.0; - n4[0] = 0.0; - n4[1] = 0.0; - n4[2] = -1.0; - n5[0] = 0.0; - n5[1] = 0.0; - n5[2] = 1.0; + n_[0] = {-1.0, 0.0, 0.0}; + n_[1] = { 1.0, 0.0, 0.0}; + n_[2] = { 0.0, -1.0, 0.0}; + n_[3] = { 0.0, 1.0, 0.0}; + n_[4] = { 0.0, 0.0, -1.0}; + n_[5] = { 0.0, 0.0, 1.0}; } /** @@ -98,94 +63,83 @@ fill(out.begin(), out.end(), 0.0); const int qOrder = 3; - const QuadratureRule& rule1 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); + const auto& rule1 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); - for (typename QuadratureRule::const_iterator it = rule1.begin(); - it != rule1.end(); ++it) + for (auto&& qp : rule1) { - Dune::FieldVector qPos = it->position(); + Dune::FieldVector qPos = qp.position(); typename LB::Traits::DomainType localPos; - localPos[0] = 0.0; - localPos[1] = qPos[0]; - localPos[2] = qPos[1]; + localPos = {0.0, qPos[0], qPos[1]}; auto y = f(localPos); - out[0] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*it->weight()*sign0; - out[6] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*it->weight(); - out[12] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[1] - 1.0)*it->weight(); - out[18] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight(); - - localPos[0] = 1.0; - localPos[1] = qPos[0]; - localPos[2] = qPos[1]; + out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1] + y[2]*n_[0][2])*qp.weight()*sign_[0]; + out[6] += (y[0]*n_[0][0] + y[1]*n_[0][1] + y[2]*n_[0][2])*(2.0*qPos[0] - 1.0)*qp.weight(); + out[12] += (y[0]*n_[0][0] + y[1]*n_[0][1] + y[2]*n_[0][2])*(2.0*qPos[1] - 1.0)*qp.weight(); + out[18] += (y[0]*n_[0][0] + y[1]*n_[0][1] + y[2]*n_[0][2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*qp.weight(); + + localPos = {1.0, qPos[0], qPos[1]}; y = f(localPos); - out[1] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*it->weight()*sign1; - out[7] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*it->weight(); - out[13] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[1])*it->weight(); - out[19] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight(); - - localPos[0] = qPos[0]; - localPos[1] = 0.0; - localPos[2] = qPos[1]; + out[1] += (y[0]*n_[1][0] + y[1]*n_[1][1] + y[2]*n_[1][2])*qp.weight()*sign_[1]; + out[7] += (y[0]*n_[1][0] + y[1]*n_[1][1] + y[2]*n_[1][2])*(1.0 - 2.0*qPos[0])*qp.weight(); + out[13] += (y[0]*n_[1][0] + y[1]*n_[1][1] + y[2]*n_[1][2])*(1.0 - 2.0*qPos[1])*qp.weight(); + out[19] += (y[0]*n_[1][0] + y[1]*n_[1][1] + y[2]*n_[1][2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*qp.weight(); + + localPos = {qPos[0], 0.0, qPos[1]}; y = f(localPos); - out[2] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*it->weight()*sign2; - out[8] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*it->weight(); - out[14] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(2.0*qPos[1] - 1.0)*it->weight(); - out[20] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight(); - - localPos[0] = qPos[0]; - localPos[1] = 1.0; - localPos[2] = qPos[1]; + out[2] += (y[0]*n_[2][0] + y[1]*n_[2][1] + y[2]*n_[2][2])*qp.weight()*sign_[2]; + out[8] += (y[0]*n_[2][0] + y[1]*n_[2][1] + y[2]*n_[2][2])*(1.0 - 2.0*qPos[0])*qp.weight(); + out[14] += (y[0]*n_[2][0] + y[1]*n_[2][1] + y[2]*n_[2][2])*(2.0*qPos[1] - 1.0)*qp.weight(); + out[20] += (y[0]*n_[2][0] + y[1]*n_[2][1] + y[2]*n_[2][2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*qp.weight(); + + localPos = {qPos[0], 1.0, qPos[1]}; y = f(localPos); - out[3] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*it->weight()*sign3; - out[9] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*it->weight(); - out[15] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(1.0 - 2.0*qPos[1])*it->weight(); - out[21] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight(); - - localPos[0] = qPos[0]; - localPos[1] = qPos[1]; - localPos[2] = 0.0; + out[3] += (y[0]*n_[3][0] + y[1]*n_[3][1] + y[2]*n_[3][2])*qp.weight()*sign_[3]; + out[9] += (y[0]*n_[3][0] + y[1]*n_[3][1] + y[2]*n_[3][2])*(2.0*qPos[0] - 1.0)*qp.weight(); + out[15] += (y[0]*n_[3][0] + y[1]*n_[3][1] + y[2]*n_[3][2])*(1.0 - 2.0*qPos[1])*qp.weight(); + out[21] += (y[0]*n_[3][0] + y[1]*n_[3][1] + y[2]*n_[3][2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*qp.weight(); + + localPos = {qPos[0], qPos[1], 0.0}; y = f(localPos); - out[4] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*it->weight()*sign4; - out[10] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*it->weight(); - out[16] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[1])*it->weight(); - out[22] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight(); - - localPos[0] = qPos[0]; - localPos[1] = qPos[1]; - localPos[2] = 1.0; + out[4] += (y[0]*n_[4][0] + y[1]*n_[4][1] + y[2]*n_[4][2])*qp.weight()*sign_[4]; + out[10] += (y[0]*n_[4][0] + y[1]*n_[4][1] + y[2]*n_[4][2])*(1.0 - 2.0*qPos[0])*qp.weight(); + out[16] += (y[0]*n_[4][0] + y[1]*n_[4][1] + y[2]*n_[4][2])*(1.0 - 2.0*qPos[1])*qp.weight(); + out[22] += (y[0]*n_[4][0] + y[1]*n_[4][1] + y[2]*n_[4][2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*qp.weight(); + + localPos = {qPos[0], qPos[1], 1.0}; y = f(localPos); - out[5] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*it->weight()*sign5; - out[11] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*it->weight(); - out[17] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[1] - 1.0)*it->weight(); - out[23] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight(); + out[5] += (y[0]*n_[5][0] + y[1]*n_[5][1] + y[2]*n_[5][2])*qp.weight()*sign_[5]; + out[11] += (y[0]*n_[5][0] + y[1]*n_[5][1] + y[2]*n_[5][2])*(2.0*qPos[0] - 1.0)*qp.weight(); + out[17] += (y[0]*n_[5][0] + y[1]*n_[5][1] + y[2]*n_[5][2])*(2.0*qPos[1] - 1.0)*qp.weight(); + out[23] += (y[0]*n_[5][0] + y[1]*n_[5][1] + y[2]*n_[5][2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*qp.weight(); } - const QuadratureRule& rule2 = QuadratureRules::rule(GeometryTypes::cube(3), qOrder); - for (typename QuadratureRule::const_iterator it = rule2.begin(); - it != rule2.end(); ++it) + const auto& rule2 = QuadratureRules::rule(GeometryTypes::cube(3), qOrder); + for (auto&& qp : rule2) { - FieldVector qPos = it->position(); + FieldVector qPos = qp.position(); auto y = f(qPos); - out[24] += y[0]*it->weight(); - out[25] += y[1]*it->weight(); - out[26] += y[2]*it->weight(); - out[27] += y[0]*qPos[1]*it->weight(); - out[28] += y[0]*qPos[2]*it->weight(); - out[29] += y[1]*qPos[0]*it->weight(); - out[30] += y[1]*qPos[2]*it->weight(); - out[31] += y[2]*qPos[0]*it->weight(); - out[32] += y[2]*qPos[1]*it->weight(); - out[33] += y[0]*qPos[1]*qPos[2]*it->weight(); - out[34] += y[1]*qPos[0]*qPos[2]*it->weight(); - out[35] += y[2]*qPos[0]*qPos[1]*it->weight(); + out[24] += y[0]*qp.weight(); + out[25] += y[1]*qp.weight(); + out[26] += y[2]*qp.weight(); + out[27] += y[0]*qPos[1]*qp.weight(); + out[28] += y[0]*qPos[2]*qp.weight(); + out[29] += y[1]*qPos[0]*qp.weight(); + out[30] += y[1]*qPos[2]*qp.weight(); + out[31] += y[2]*qPos[0]*qp.weight(); + out[32] += y[2]*qPos[1]*qp.weight(); + out[33] += y[0]*qPos[1]*qPos[2]*qp.weight(); + out[34] += y[1]*qPos[0]*qPos[2]*qp.weight(); + out[35] += y[2]*qPos[0]*qPos[1]*qp.weight(); } } private: - typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3, sign4, sign5; - typename LB::Traits::DomainType n0, n1, n2, n3, n4, n5; + // Facet orientations + std::array sign_; + + // Facet normals + std::array n_; }; } #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALINTERPOLATION_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas2cube2d/raviartthomas2cube2dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas2cube2d/raviartthomas2cube2dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas2cube2d/raviartthomas2cube2dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas2cube2d/raviartthomas2cube2dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -30,34 +30,15 @@ * * \param s Edge orientation indicator */ - RT2Cube2DLocalInterpolation (unsigned int s = 0) + RT2Cube2DLocalInterpolation (std::bitset<4> s = 0) { - sign0 = sign1 = sign2 = sign3 = 1.0; - if (s & 1) - { - sign0 *= -1.0; - } - if (s & 2) - { - sign1 *= -1.0; - } - if (s & 4) - { - sign2 *= -1.0; - } - if (s & 8) - { - sign3 *= -1.0; - } + for (size_t i=0; i<4; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; - n0[0] = -1.0; - n0[1] = 0.0; - n1[0] = 1.0; - n1[1] = 0.0; - n2[0] = 0.0; - n2[1] = -1.0; - n3[0] = 0.0; - n3[1] = 1.0; + n_[0] = {-1.0, 0.0}; + n_[1] = { 1.0, 0.0}; + n_[2] = { 0.0, -1.0}; + n_[3] = { 0.0, 1.0}; } /** @@ -81,68 +62,66 @@ fill(out.begin(), out.end(), 0.0); const int qOrder = 6; - const QuadratureRule& rule = QuadratureRules::rule(GeometryTypes::cube(1), qOrder); + const auto& rule1 = QuadratureRules::rule(GeometryTypes::cube(1), qOrder); - for (typename QuadratureRule::const_iterator it=rule.begin(); it!=rule.end(); ++it) + for (auto&& qp : rule1) { - Scalar qPos = it->position(); + Scalar qPos = qp.position(); typename LB::Traits::DomainType localPos; - localPos[0] = 0.0; - localPos[1] = qPos; + localPos = {0.0, qPos}; auto y = f(localPos); - out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0; - out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight(); - out[2] += (y[0]*n0[0] + y[1]*n0[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign0; + out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]; + out[1] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight(); + out[2] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[0]; - localPos[0] = 1.0; - localPos[1] = qPos; + localPos = {1.0, qPos}; y = f(localPos); - out[3] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1; - out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight(); - out[5] += (y[0]*n1[0] + y[1]*n1[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign1; + out[3] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]; + out[4] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight(); + out[5] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[1]; - localPos[0] = qPos; - localPos[1] = 0.0; + localPos = {qPos, 0.0}; y = f(localPos); - out[6] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2; - out[7] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight(); - out[8] += (y[0]*n2[0] + y[1]*n2[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign2; + out[6] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]; + out[7] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(1.0 - 2.0*qPos)*qp.weight(); + out[8] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[2]; - localPos[0] = qPos; - localPos[1] = 1.0; + localPos = {qPos, 1.0}; y = f(localPos); - out[9] += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3; - out[10] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight(); - out[11] += (y[0]*n3[0] + y[1]*n3[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign3; + out[9] += (y[0]*n_[3][0] + y[1]*n_[3][1])*qp.weight()*sign_[3]; + out[10] += (y[0]*n_[3][0] + y[1]*n_[3][1])*(2.0*qPos - 1.0)*qp.weight(); + out[11] += (y[0]*n_[3][0] + y[1]*n_[3][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[3]; } - const QuadratureRule& rule2 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); + const auto& rule2 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); - for (typename QuadratureRule::const_iterator it = rule2.begin(); - it != rule2.end(); ++it) + for (auto&& qp : rule2) { - FieldVector qPos = it->position(); + FieldVector qPos = qp.position(); auto y = f(qPos); - out[12] += y[0]*it->weight(); - out[13] += y[1]*it->weight(); - out[14] += y[0]*qPos[0]*it->weight(); - out[15] += y[1]*qPos[0]*it->weight(); - out[16] += y[0]*qPos[1]*it->weight(); - out[17] += y[1]*qPos[1]*it->weight(); - out[18] += y[0]*qPos[0]*qPos[1]*it->weight(); - out[19] += y[1]*qPos[0]*qPos[1]*it->weight(); - out[20] += y[0]*qPos[1]*qPos[1]*it->weight(); - out[21] += y[1]*qPos[0]*qPos[0]*it->weight(); - out[22] += y[0]*qPos[0]*qPos[1]*qPos[1]*it->weight(); - out[23] += y[1]*qPos[0]*qPos[0]*qPos[1]*it->weight(); + out[12] += y[0]*qp.weight(); + out[13] += y[1]*qp.weight(); + out[14] += y[0]*qPos[0]*qp.weight(); + out[15] += y[1]*qPos[0]*qp.weight(); + out[16] += y[0]*qPos[1]*qp.weight(); + out[17] += y[1]*qPos[1]*qp.weight(); + out[18] += y[0]*qPos[0]*qPos[1]*qp.weight(); + out[19] += y[1]*qPos[0]*qPos[1]*qp.weight(); + out[20] += y[0]*qPos[1]*qPos[1]*qp.weight(); + out[21] += y[1]*qPos[0]*qPos[0]*qp.weight(); + out[22] += y[0]*qPos[0]*qPos[1]*qPos[1]*qp.weight(); + out[23] += y[1]*qPos[0]*qPos[0]*qPos[1]*qp.weight(); } } private: - typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3; - typename LB::Traits::DomainType n0, n1, n2, n3; + // Edge orientations + std::array sign_; + + // Edge normals + std::array n_; }; } #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2_CUBE2D_LOCALINTERPOLATION_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas3cube2d/raviartthomas3cube2dlocalinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas3cube2d/raviartthomas3cube2dlocalinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomas3cube2d/raviartthomas3cube2dlocalinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomas3cube2d/raviartthomas3cube2dlocalinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -30,34 +30,15 @@ * * \param s Edge orientation indicator */ - RT3Cube2DLocalInterpolation (unsigned int s = 0) + RT3Cube2DLocalInterpolation (std::bitset<4> s = 0) { - sign0 = sign1 = sign2 = sign3 = 1.0; - if (s & 1) - { - sign0 *= -1.0; - } - if (s & 2) - { - sign1 *= -1.0; - } - if (s & 4) - { - sign2 *= -1.0; - } - if (s & 8) - { - sign3 *= -1.0; - } + for (size_t i=0; i<4; i++) + sign_[i] = (s[i]) ? -1.0 : 1.0; - n0[0] = -1.0; - n0[1] = 0.0; - n1[0] = 1.0; - n1[1] = 0.0; - n2[0] = 0.0; - n2[1] = -1.0; - n3[0] = 0.0; - n3[1] = 1.0; + n_[0] = {-1.0, 0.0}; + n_[1] = { 1.0, 0.0}; + n_[2] = { 0.0, -1.0}; + n_[3] = { 0.0, 1.0}; } /** @@ -81,52 +62,47 @@ fill(out.begin(), out.end(), 0.0); const int qOrder = 9; - const QuadratureRule& rule = QuadratureRules::rule(GeometryTypes::cube(1), qOrder); + const auto& rule1 = QuadratureRules::rule(GeometryTypes::cube(1), qOrder); - for (typename QuadratureRule::const_iterator it=rule.begin(); it!=rule.end(); ++it) + for (auto&& qp : rule1) { - Scalar qPos = it->position(); + Scalar qPos = qp.position(); typename LB::Traits::DomainType localPos; - localPos[0] = 0.0; - localPos[1] = qPos; + localPos = {0.0, qPos}; auto y = f(localPos); - out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0; - out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight(); - out[2] += (y[0]*n0[0] + y[1]*n0[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign0; - out[3] += (y[0]*n0[0] + y[1]*n0[1])*(20.0*qPos*qPos*qPos - 30.0*qPos*qPos + 12.0*qPos - 1.0)*it->weight(); + out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]; + out[1] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight(); + out[2] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[0]; + out[3] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(20.0*qPos*qPos*qPos - 30.0*qPos*qPos + 12.0*qPos - 1.0)*qp.weight(); - localPos[0] = 1.0; - localPos[1] = qPos; + localPos = {1.0, qPos}; y = f(localPos); - out[4] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1; - out[5] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight(); - out[6] += (y[0]*n1[0] + y[1]*n1[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign1; - out[7] += (y[0]*n1[0] + y[1]*n1[1])*(-20.0*qPos*qPos*qPos + 30.0*qPos*qPos - 12.0*qPos + 1.0)*it->weight(); + out[4] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]; + out[5] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight(); + out[6] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[1]; + out[7] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(-20.0*qPos*qPos*qPos + 30.0*qPos*qPos - 12.0*qPos + 1.0)*qp.weight(); - localPos[0] = qPos; - localPos[1] = 0.0; + localPos = {qPos, 0.0}; y = f(localPos); - out[8] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2; - out[9] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight(); - out[10] += (y[0]*n2[0] + y[1]*n2[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign2; - out[11] += (y[0]*n2[0] + y[1]*n2[1])*(-20.0*qPos*qPos*qPos + 30.0*qPos*qPos - 12.0*qPos + 1.0)*it->weight(); + out[8] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]; + out[9] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(1.0 - 2.0*qPos)*qp.weight(); + out[10] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[2]; + out[11] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(-20.0*qPos*qPos*qPos + 30.0*qPos*qPos - 12.0*qPos + 1.0)*qp.weight(); - localPos[0] = qPos; - localPos[1] = 1.0; + localPos = {qPos, 1.0}; y = f(localPos); - out[12] += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3; - out[13] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight(); - out[14] += (y[0]*n3[0] + y[1]*n3[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign3; - out[15] += (y[0]*n3[0] + y[1]*n3[1])*(20.0*qPos*qPos*qPos - 30.0*qPos*qPos + 12.0*qPos - 1.0)*it->weight(); + out[12] += (y[0]*n_[3][0] + y[1]*n_[3][1])*qp.weight()*sign_[3]; + out[13] += (y[0]*n_[3][0] + y[1]*n_[3][1])*(2.0*qPos - 1.0)*qp.weight(); + out[14] += (y[0]*n_[3][0] + y[1]*n_[3][1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*qp.weight()*sign_[3]; + out[15] += (y[0]*n_[3][0] + y[1]*n_[3][1])*(20.0*qPos*qPos*qPos - 30.0*qPos*qPos + 12.0*qPos - 1.0)*qp.weight(); } - const QuadratureRule& rule2 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); + const auto& rule2 = QuadratureRules::rule(GeometryTypes::cube(2), qOrder); - for (typename QuadratureRule::const_iterator it = rule2.begin(); - it != rule2.end(); ++it) + for (auto&& qp : rule2) { - FieldVector qPos = it->position(); + auto qPos = qp.position(); auto y = f(qPos); double l0_x=1.0; @@ -138,37 +114,40 @@ double l2_y=6.0*qPos[1]*qPos[1]-6.0*qPos[1]+1.0; double l3_y=20.0*qPos[1]*qPos[1]*qPos[1] - 30.0*qPos[1]*qPos[1] + 12.0*qPos[1] - 1.0; - out[16] += y[0]*l0_x*l0_y*it->weight(); - out[17] += y[0]*l0_x*l1_y*it->weight(); - out[18] += y[0]*l0_x*l2_y*it->weight(); - out[19] += y[0]*l0_x*l3_y*it->weight(); - out[20] += y[0]*l1_x*l0_y*it->weight(); - out[21] += y[0]*l1_x*l1_y*it->weight(); - out[22] += y[0]*l1_x*l2_y*it->weight(); - out[23] += y[0]*l1_x*l3_y*it->weight(); - out[24] += y[0]*l2_x*l0_y*it->weight(); - out[25] += y[0]*l2_x*l1_y*it->weight(); - out[26] += y[0]*l2_x*l2_y*it->weight(); - out[27] += y[0]*l2_x*l3_y*it->weight(); - - out[28] += y[1]*l0_x*l0_y*it->weight(); - out[29] += y[1]*l0_x*l1_y*it->weight(); - out[30] += y[1]*l0_x*l2_y*it->weight(); - out[31] += y[1]*l1_x*l0_y*it->weight(); - out[32] += y[1]*l1_x*l1_y*it->weight(); - out[33] += y[1]*l1_x*l2_y*it->weight(); - out[34] += y[1]*l2_x*l0_y*it->weight(); - out[35] += y[1]*l2_x*l1_y*it->weight(); - out[36] += y[1]*l2_x*l2_y*it->weight(); - out[37] += y[1]*l3_x*l0_y*it->weight(); - out[38] += y[1]*l3_x*l1_y*it->weight(); - out[39] += y[1]*l3_x*l2_y*it->weight(); + out[16] += y[0]*l0_x*l0_y*qp.weight(); + out[17] += y[0]*l0_x*l1_y*qp.weight(); + out[18] += y[0]*l0_x*l2_y*qp.weight(); + out[19] += y[0]*l0_x*l3_y*qp.weight(); + out[20] += y[0]*l1_x*l0_y*qp.weight(); + out[21] += y[0]*l1_x*l1_y*qp.weight(); + out[22] += y[0]*l1_x*l2_y*qp.weight(); + out[23] += y[0]*l1_x*l3_y*qp.weight(); + out[24] += y[0]*l2_x*l0_y*qp.weight(); + out[25] += y[0]*l2_x*l1_y*qp.weight(); + out[26] += y[0]*l2_x*l2_y*qp.weight(); + out[27] += y[0]*l2_x*l3_y*qp.weight(); + + out[28] += y[1]*l0_x*l0_y*qp.weight(); + out[29] += y[1]*l0_x*l1_y*qp.weight(); + out[30] += y[1]*l0_x*l2_y*qp.weight(); + out[31] += y[1]*l1_x*l0_y*qp.weight(); + out[32] += y[1]*l1_x*l1_y*qp.weight(); + out[33] += y[1]*l1_x*l2_y*qp.weight(); + out[34] += y[1]*l2_x*l0_y*qp.weight(); + out[35] += y[1]*l2_x*l1_y*qp.weight(); + out[36] += y[1]*l2_x*l2_y*qp.weight(); + out[37] += y[1]*l3_x*l0_y*qp.weight(); + out[38] += y[1]*l3_x*l1_y*qp.weight(); + out[39] += y[1]*l3_x*l2_y*qp.weight(); } } private: - typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3; - typename LB::Traits::DomainType n0, n1, n2, n3; + // Edge orientations + std::array sign_; + + // Edge normals + std::array n_; }; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomaslfecache.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomaslfecache.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomaslfecache.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomaslfecache.hh 2021-08-31 14:03:38.000000000 +0000 @@ -67,6 +67,7 @@ static auto getImplementations() { return std::make_tuple( + std::make_pair(index(GeometryTypes::tetrahedron), []() { return RT03DLocalFiniteElement(); }), std::make_pair(index(GeometryTypes::hexahedron), []() { return RT0Cube3DLocalFiniteElement(); }) ); } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexbasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -12,6 +12,24 @@ namespace Dune { + /* + * `RTPreBasisFactory` provides a basis for the Raviart-Thomas function space. + * `RaviartThomasL2InterpolationFactory` provides the linear functionals. + * + * `Defaultbasisfactory::create` first builds the function space and the linear functionals. + * Then the constructor of `BasisMatrix` gets called. There the matrix + * + * \begin{equation} + * A_{i,j} := N_j(\phi_i) + * \end{equation} + * + * with linear functionals $N_j$ and basisfunctions $\phi_i$ gets assembled. + * Then the matrix gets inverted and is then used as a coefficent matrix for the standard monomial basis. + * + * For more details on the theory see the first chapter "Construction of Local Finite Element Spaces Using the Generic Reference Elements" + * of the book "Advances in Dune" by Dedner, Flemisch and Klöfkorn published in 2012. + */ + template< unsigned int dim, class SF, class CF > struct RaviartThomasBasisFactory : public DefaultBasisFactory< RTPreBasisFactory, diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexinterpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexinterpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexinterpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexinterpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -68,22 +69,22 @@ typedef std::size_t Key; typedef const LocalCoefficientsContainer Object; - template< class Topology > + template< GeometryType::Id geometryId > static Object *create( const Key &key ) { typedef RaviartThomasL2InterpolationFactory< dim, double > InterpolationFactory; - if( !supports< Topology >( key ) ) + if( !supports< geometryId >( key ) ) return nullptr; - typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< Topology >( key ); + typename InterpolationFactory::Object *interpolation = InterpolationFactory::template create< geometryId >( key ); Object *localKeys = new Object( *interpolation ); InterpolationFactory::release( interpolation ); return localKeys; } - template< class Topology > + template< GeometryType::Id geometryId > static bool supports ( const Key &key ) { - return Impl::IsSimplex< Topology >::value; + return GeometryType(geometryId).isSimplex(); } static void release( Object *object ) { delete object; } }; @@ -103,10 +104,16 @@ struct RTL2InterpolationBuilder { static const unsigned int dimension = dim; + + // for the dofs associated to the element typedef OrthonormalBasisFactory< dimension, Field > TestBasisFactory; typedef typename TestBasisFactory::Object TestBasis; + + // for the dofs associated to the faces typedef OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory; typedef typename TestFaceBasisFactory::Object TestFaceBasis; + + // the normals of the faces typedef FieldVector< Field, dimension > Normal; RTL2InterpolationBuilder () = default; @@ -121,33 +128,53 @@ TestFaceBasisFactory::release( f.basis_ ); } - unsigned int topologyId () const { return topologyId_; } + [[deprecated("Use type().id() instead.")]] + unsigned int topologyId () const { return type().id(); } - GeometryType type () const { return GeometryType( topologyId(), dimension ); } + GeometryType type () const { return geometry_; } std::size_t order () const { return order_; } + // number of faces unsigned int faceSize () const { return faceSize_; } + // basis associated to the element TestBasis *testBasis () const { return testBasis_; } + + // basis associated to face f TestFaceBasis *testFaceBasis ( unsigned int f ) const { assert( f < faceSize() ); return faceStructure_[ f ].basis_; } + // normal of face f const Normal &normal ( unsigned int f ) const { assert( f < faceSize() ); return *(faceStructure_[ f ].normal_); } - template< class Topology > + template< GeometryType::Id geometryId > void build ( std::size_t order ) { + constexpr GeometryType geometry = geometryId; + geometry_ = geometry; order_ = order; - topologyId_ = Topology::id; - testBasis_ = (order > 0 ? TestBasisFactory::template create< Topology >( order-1 ) : nullptr); + testBasis_ = (order > 0 ? TestBasisFactory::template create< geometry >( order-1 ) : nullptr); const auto &refElement = ReferenceElements< Field, dimension >::general( type() ); faceSize_ = refElement.size( 1 ); faceStructure_.reserve( faceSize_ ); for( unsigned int face = 0; face < faceSize_; ++face ) { - TestFaceBasis *faceBasis = Impl::IfTopology< CreateFaceBasis, dimension-1 >::apply( refElement.type( face, 1 ).id(), order ); + /* For simplices or cubes of arbitrary dimension you could just use + * + * ``` + * GeometryType faceGeometry = Impl::getBase(geometry_); + * TestFaceBasis *faceBasis = TestFaceBasisFactory::template create< faceGeometry >( order ); + * ``` + * + * For i.e. Prisms and Pyramids in 3d this does not work because they contain squares and triangles as faces. + * And depending on the dynamic face index a different face geometry is needed. + * + */ + TestFaceBasis *faceBasis = Impl::toGeometryTypeIdConstant(refElement.type( face, 1 ), [&](auto faceGeometryTypeId) { + return TestFaceBasisFactory::template create< decltype(faceGeometryTypeId)::value >( order ); + }); faceStructure_.emplace_back( faceBasis, refElement.integrationOuterNormal( face ) ); } assert( faceStructure_.size() == faceSize_ ); @@ -164,15 +191,10 @@ const Dune::FieldVector< Field, dimension > *normal_; }; - template< class FaceTopology > - struct CreateFaceBasis - { - static TestFaceBasis *apply ( std::size_t order ) { return TestFaceBasisFactory::template create< FaceTopology >( order ); } - }; - std::vector< FaceStructure > faceStructure_; TestBasis *testBasis_ = nullptr; - unsigned int topologyId_, faceSize_; + GeometryType geometry_; + unsigned int faceSize_; std::size_t order_; }; @@ -228,12 +250,12 @@ { return size_; } - template + template void build( std::size_t order ) { size_ = 0; order_ = order; - builder_.template build(order_); + builder_.template build(order_); if (builder_.testBasis()) size_ += dimension*builder_.testBasis()->size(); for ( unsigned int f=0; f void interpolate ( typename Base::template Helper &func ) const { - const Dune::GeometryType geoType( builder_.topologyId(), dimension ); + const Dune::GeometryType geoType = builder_.type(); std::vector< Field > testBasisVal; @@ -327,7 +349,15 @@ } private: - /** /brief evaluate boundary functionals **/ + /** \brief evaluate boundary functionals + * + * \param startRow row of matrix to start + * \param mVal value of the testBasis at a quadrature point on a face + * \param rtVal value of the RaviartThomasBasis at a quadrature point on a face + * \param normal the normal of the face + * \param weight quadrature weight + * \param matrix result gets written into matrix starting with row: row + */ template void fillBnd (unsigned int startRow, const MVal &mVal, @@ -350,6 +380,14 @@ assert( miter == mVal.end() ); } } + /** \brief evaluate interior functionals + * + * \param startRow row of matrix to start + * \param mVal value of the testBasis at a quadrature point in the interior of the ReferenceElement + * \param rtVal value of the RaviartThomasBasis at a quadrature point in the interior of the ReferenceElement + * \param weight quadrature weight + * \param matrix result gets written into matrix starting with row: row + */ template void fillInterior (unsigned int startRow, const MVal &mVal, @@ -386,19 +424,20 @@ typedef const RaviartThomasL2Interpolation Object; typedef std::size_t Key; typedef typename std::remove_const::type NonConstObject; - template + + template static Object *create( const Key &key ) { - if ( !supports(key) ) + if ( !supports(key) ) return 0; NonConstObject *interpol = new NonConstObject(); - interpol->template build(key); + interpol->template build(key); return interpol; } - template< class Topology > + template< GeometryType::Id geometryId > static bool supports ( const Key &key ) { - return Impl::IsSimplex::value; + return GeometryType(geometryId).isSimplex(); } static void release( Object *object ) { delete object; } }; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexprebasis.hh dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexprebasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexprebasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/raviartthomas/raviartthomassimplex/raviartthomassimplexprebasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -12,7 +12,7 @@ namespace Dune { - template < class Topology, class Field > + template < GeometryType::Id geometryId, class Field > struct RTVecMatrix; template @@ -31,48 +31,105 @@ { typedef MonomialBasisProvider Type; }; - template< class Topology > + template< GeometryType::Id geometryId > static Object *create ( const Key &order ) { - RTVecMatrix vecMatrix(order); - MBasis *mbasis = MBasisFactory::template create(order+1); - typename std::remove_const::type *tmBasis = - new typename std::remove_const::type(*mbasis); + RTVecMatrix vecMatrix(order); + MBasis *mbasis = MBasisFactory::template create(order+1); + typename std::remove_const::type *tmBasis = new typename std::remove_const::type(*mbasis); tmBasis->fill(vecMatrix); return tmBasis; } static void release( Object *object ) { delete object; } }; - template + + template struct RTVecMatrix { - static const unsigned int dim = Topology::dimension; + static constexpr GeometryType geometry = geometryId; + static const unsigned int dim = geometry.dim(); typedef MultiIndex MI; - typedef MonomialBasis MIBasis; + typedef MonomialBasis MIBasis; RTVecMatrix(std::size_t order) { + /* + * Construction of Raviart-Thomas elements in high dimensions see "Mixed Finite Elements in \R^3" by Nedelec, 1980. + * + * Let $\P_{n,k}$ be the space of polynomials in $n$ variables with degree $\leq k$. + * The space of Raviart-Thomas functions in $n$ dimensions with index $k$ is defined as + * + * \begin{equation*} + * RT_k := (\P_{k-1})^n \oplus \widetilde \P_k x + * \end{equation*} + * with $x=(x_1,x_2,\dots, x_n)$ in $n$ dimensions and $\widetilde \P_k$ the homogeneous polynomials of degree $k$. + * + * For $RT_k$ holds + * \begin{equation*} + * (\P_{k-1})^n \subset RT_k \subset (\P_k)^n. + * \end{equation*} + * + * We construct $(\P_k)^n$ and and only use the monomials contained in $RT_k$. + * + */ + MIBasis basis(order+1); FieldVector< MI, dim > x; + /* + * Init MultiIndices + * x[0]=(1,0,0) x + * x[1]=(0,1,0) y + * x[2]=(0,0,1) z + */ for( unsigned int i = 0; i < dim; ++i ) x[ i ].set( i, 1 ); std::vector< MI > val( basis.size() ); + + // val now contains all monomials in $n$ dimensions with degree $\leq order+1$ basis.evaluate( x, val ); col_ = basis.size(); + + // get $\dim (\P_{order-1})$ unsigned int notHomogen = 0; if (order>0) notHomogen = basis.sizes()[order-1]; + + // get $\dim \widetilde (\P_order)$ unsigned int homogen = basis.sizes()[order]-notHomogen; + + /* + * + * The set $RT_k$ is defined as + * + * \begin{equation} + * RT_k := (\P_k)^dim + \widetilde \P_k x with x\in \R^n. + * \end{equation} + * + * The space $\P_k$ is split in $\P_k = \P_{k-1} + \widetilde \P_k$. + * + * \begin{align} + * RT_k &= (\P_{k-1} + \widetilde \P_k)^dim + \widetilde \P_k x with x\in \R^n + * &= (\P_{k-1})^n + (\widetilde \P_k)^n + \widetilde \P_k x with x\in \R^n + * \end{align} + * + * Thus $\dim RT_k = n * \dim \P_{k-1} + (n+1)*\dim (\widetilde \P_k)$ + */ + + // row_ = \dim RT_k *dim row_ = (notHomogen*dim+homogen*(dim+1))*dim; - row1_ = basis.sizes()[order]*dim*dim; mat_ = new Field*[row_]; int row = 0; + + /* Assign the correct values for the nonhomogeneous polymonials $p\in (\P_{oder-1})^dim$ + * A basis function is represented by $dim$ rows. + */ for (unsigned int i=0; i void row( const unsigned int row, Vector &vec ) const { @@ -130,7 +198,7 @@ for (unsigned int i=0; i #include +#include #include #include #include diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/brezzidouglasmarinielementtest.cc dune-localfunctions-2.8.0/dune/localfunctions/test/brezzidouglasmarinielementtest.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/brezzidouglasmarinielementtest.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/brezzidouglasmarinielementtest.cc 2021-08-31 14:03:38.000000000 +0000 @@ -17,7 +17,9 @@ TEST_FE(bdm1cube2dlfem); Dune::BrezziDouglasMariniCubeLocalFiniteElement bdm1cube3dlfem(1); - TEST_FE2(bdm1cube3dlfem, DisableLocalInterpolation); + // \todo Implement the missing LocalInterpolation + // DisableRepresentConstants is only set because the test also uses DisableLocalInterpolation internally. + TEST_FE2(bdm1cube3dlfem, DisableLocalInterpolation + DisableRepresentConstants); Dune::BrezziDouglasMariniCubeLocalFiniteElement bdm2cube2dlfem(1); TEST_FE(bdm2cube2dlfem); diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/CMakeLists.txt dune-localfunctions-2.8.0/dune/localfunctions/test/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/localfunctions/test/CMakeLists.txt 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -1,3 +1,5 @@ +add_definitions(-DDUNE_DEPRECATED_INTERPOLATE_CHECK=1) + dune_add_test(SOURCES bdfmelementtest.cc) dune_add_test(SOURCES brezzidouglasmarinielementtest.cc) @@ -14,6 +16,8 @@ dune_add_test(SOURCES monomialshapefunctiontest.cc) +dune_add_test(SOURCES nedelec1stkindelementtest.cc) + dune_add_test(SOURCES rannacherturekelementtest.cc) dune_add_test(SOURCES raviartthomaselementtest.cc) @@ -84,3 +88,11 @@ dune_add_test(NAME test-raviartthomassimplex4 SOURCES test-raviartthomassimplex.cc COMPILE_DEFINITIONS "CHECKDIM=4") + +dune_add_test(NAME test-nedelecsimplex2 + SOURCES test-nedelecsimplex.cc + COMPILE_DEFINITIONS "CHECKDIM=2") + +dune_add_test(NAME test-nedelecsimplex3 + SOURCES test-nedelecsimplex.cc + COMPILE_DEFINITIONS "CHECKDIM=3") diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/globalmonomialfunctionstest.cc dune-localfunctions-2.8.0/dune/localfunctions/test/globalmonomialfunctionstest.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/globalmonomialfunctionstest.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/globalmonomialfunctionstest.cc 2021-08-31 14:03:38.000000000 +0000 @@ -12,7 +12,6 @@ #include #include #include -#include #include #include @@ -56,14 +55,14 @@ template static void Dim(int &result) { - Dune::Hybrid::forEach(Dune::Std::make_index_sequence<4>{},[&](auto i){Order(result);}); + Dune::Hybrid::forEach(std::make_index_sequence<4>{},[&](auto i){Order(result);}); } int main(int argc, char** argv) { try { int result = 77; - Dune::Hybrid::forEach(Dune::Std::make_index_sequence<3>{},[&](auto i){Dim(result);}); + Dune::Hybrid::forEach(std::make_index_sequence<3>{},[&](auto i){Dim(result);}); return result; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/lagrangeshapefunctiontest.cc dune-localfunctions-2.8.0/dune/localfunctions/test/lagrangeshapefunctiontest.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/lagrangeshapefunctiontest.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/lagrangeshapefunctiontest.cc 2021-08-31 14:03:38.000000000 +0000 @@ -185,11 +185,15 @@ PrismP2LocalFiniteElement prismp2fem; TEST_FE3(prismp2fem, DisableNone, 1); + // Pyramid shapefunctions are not differentiable on the plane where xi[0]=xi[1]. + // So let's skip test points on this plane + auto xySkip = [](const FieldVector& xi){return std::abs(xi[0]-xi[1]) < 1e-8;}; + PyramidP1LocalFiniteElement pyramidp1fem; - TEST_FE2(pyramidp1fem, DisableJacobian); + TEST_FE4(pyramidp1fem, DisableNone, 1, xySkip); PyramidP2LocalFiniteElement pyramidp2fem; - TEST_FE2(pyramidp2fem, DisableJacobian); + TEST_FE4(pyramidp2fem, DisableNone, 1, xySkip); Hybrid::forEach(std::make_index_sequence<4>{},[&success](auto i) { diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/nedelec1stkindelementtest.cc dune-localfunctions-2.8.0/dune/localfunctions/test/nedelec1stkindelementtest.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/nedelec1stkindelementtest.cc 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/nedelec1stkindelementtest.cc 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,58 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: + +#include "config.h" + +#include +#include + +#include + +using namespace Dune; + +int main(int argc, char** argv) +{ + bool success = true; + + // First order on a triangle + Nedelec1stKindSimplexLocalFiniteElement nedelecLFEMTriangle1stOrder; + TEST_FE3(nedelecLFEMTriangle1stOrder, DisableNone, 2); + + for (unsigned int s = 0; s < 8; s++) + { + Nedelec1stKindSimplexLocalFiniteElement nedelecLFEMTriangle1stOrder(s); + TEST_FE3(nedelecLFEMTriangle1stOrder, DisableNone, 2); + } + + // First order on a tetrahedron + Nedelec1stKindSimplexLocalFiniteElement nedelecLFEMTetrahedron1stOrder; + TEST_FE3(nedelecLFEMTetrahedron1stOrder, DisableNone, 2); + + for (unsigned int s = 0; s < 64; s++) + { + Nedelec1stKindSimplexLocalFiniteElement nedelecLFEMTetrahedron1stOrder(s); + TEST_FE3(nedelecLFEMTetrahedron1stOrder, DisableNone, 2); + } + + // First order on a square + Nedelec1stKindCubeLocalFiniteElement nedelecLFEMSquare1stOrder; + TEST_FE3(nedelecLFEMSquare1stOrder, DisableNone, 2); + + for (unsigned int s = 0; s < 16; s++) + { + Nedelec1stKindCubeLocalFiniteElement nedelecLFEMSquare1stOrder(s); + TEST_FE3(nedelecLFEMSquare1stOrder, DisableNone, 2); + } + + // First order on a cube + Nedelec1stKindCubeLocalFiniteElement nedelecLFEMCube1stOrder; + TEST_FE3(nedelecLFEMCube1stOrder, DisableNone, 2); + + for (unsigned int s = 0; s < 4096; s++) + { + Nedelec1stKindCubeLocalFiniteElement nedelecLFEMCube1stOrder(s); + TEST_FE3(nedelecLFEMCube1stOrder, DisableNone, 2); + } + + return success ? 0 : 1; +} diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/raviartthomaselementtest.cc dune-localfunctions-2.8.0/dune/localfunctions/test/raviartthomaselementtest.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/raviartthomaselementtest.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/raviartthomaselementtest.cc 2021-08-31 14:03:38.000000000 +0000 @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -23,6 +24,9 @@ Dune::RaviartThomasSimplexLocalFiniteElement<2,double,double> rt1simplex2dlfem(Dune::GeometryTypes::simplex(2),1); TEST_FE(rt1simplex2dlfem); + Dune::RaviartThomasSimplexLocalFiniteElement<3,double,double> rt0simplex3dlfem(Dune::GeometryTypes::simplex(3),0); + TEST_FE(rt0simplex3dlfem); + Dune::RaviartThomasCubeLocalFiniteElement rt0cube2dlfem; TEST_FE(rt0cube2dlfem); for (unsigned int s = 0; s < 16; s++) @@ -111,6 +115,14 @@ TEST_FE(rt12dlfemDedicated); } + Dune::RT03DLocalFiniteElement rt03dlfemDedicated; + TEST_FE(rt03dlfemDedicated); + for (unsigned int s = 0; s < 16; s++) + { + Dune::RT03DLocalFiniteElement rt03dlfemDedicated(s); + TEST_FE(rt03dlfemDedicated); + } + Dune::RT0Cube3DLocalFiniteElement rt0cube3dlfemDedicated; TEST_FE(rt0cube3dlfemDedicated); for (unsigned int s = 0; s < 64; s++) diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-fe.hh dune-localfunctions-2.8.0/dune/localfunctions/test/test-fe.hh --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-fe.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-fe.hh 2021-08-31 14:03:38.000000000 +0000 @@ -20,7 +20,6 @@ #include #include -#include #include @@ -30,15 +29,19 @@ // This provides the evaluate method needed by the interpolate() // method. template -class FEFunction : - public Dune::Function +class FEFunction { - typedef typename FE::Traits::Basis::Traits::DomainLocal DomainLocal; - typedef typename FE::Traits::Basis::Traits::Range Range; - const FE& fe; public: + using RangeType = typename FE::Traits::Basis::Traits::Range; + using DomainType = typename FE::Traits::Basis::Traits::DomainLocal; + + struct Traits { + using RangeType = typename FE::Traits::Basis::Traits::Range; + using DomainType = typename FE::Traits::Basis::Traits::DomainLocal; + }; + typedef typename FE::Traits::Basis::Traits::RangeField CT; std::vector coeff; @@ -57,10 +60,9 @@ coeff[i] = ((1.0*std::rand()) / RAND_MAX - 0.5)*2.0*max; } - void evaluate (const DomainLocal& x, Range& y) const { - std::vector yy; + void evaluate (const DomainType& x, RangeType& y) const { + std::vector yy; fe.basis().evaluateFunction(x, yy); - y = 0.0; for (std::size_t i=0; i f(fe); std::vector::CT> coeff; + std::vector::CT> coeff2; for(int i=0; i() << ":" << std::endl; + success = false; + } + // Check size of weight vector if (coeff.size() != fe.basis().size()) { std::cout << "Bug in LocalInterpolation for finite element type " diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-finiteelementcache.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-finiteelementcache.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-finiteelementcache.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-finiteelementcache.cc 2021-08-31 14:03:38.000000000 +0000 @@ -5,13 +5,15 @@ #include "config.h" #endif +#include + #include -#include #include #include #include +#include template static void test(Dune::GeometryType type) @@ -19,12 +21,12 @@ FiniteElementCache cache; using FiniteElement = typename FiniteElementCache::FiniteElementType; - DUNE_UNUSED const FiniteElement& finiteElement = cache.get(type); + [[maybe_unused]] const FiniteElement& finiteElement = cache.get(type); } int main() { static constexpr std::size_t max_k = 3; - Dune::Hybrid::forEach(Dune::Std::make_index_sequence{},[&](auto k) + Dune::Hybrid::forEach(std::make_index_sequence{},[&](auto k) { constexpr int dim = 2; using FiniteElementCache = typename @@ -39,6 +41,33 @@ test(Dune::GeometryTypes::simplex(dim)); test(Dune::GeometryTypes::cube(dim)); } + + { + constexpr int dim = 2; + constexpr int order = 0; + using FiniteElementCache = typename + Dune::RaviartThomasLocalFiniteElementCache; + test(Dune::GeometryTypes::simplex(dim)); + test(Dune::GeometryTypes::cube(dim)); + } + + { + constexpr int dim = 2; + constexpr int order = 1; + using FiniteElementCache = typename + Dune::RaviartThomasLocalFiniteElementCache; + test(Dune::GeometryTypes::simplex(dim)); + test(Dune::GeometryTypes::cube(dim)); + } + + { + constexpr int dim = 3; + constexpr int order = 0; + using FiniteElementCache = typename + Dune::RaviartThomasLocalFiniteElementCache; + test(Dune::GeometryTypes::simplex(dim)); + test(Dune::GeometryTypes::cube(dim)); + } return 0; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/testgenericfem.cc dune-localfunctions-2.8.0/dune/localfunctions/test/testgenericfem.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/testgenericfem.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/testgenericfem.cc 2021-08-31 14:03:38.000000000 +0000 @@ -11,8 +11,6 @@ #include #include -#include - #include // Lagrange type elements diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-lagrange.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-lagrange.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-lagrange.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-lagrange.cc 2021-08-31 14:03:38.000000000 +0000 @@ -2,6 +2,8 @@ // vi: set et ts=4 sw=2 sts=2: #include +#include + #include #include @@ -14,9 +16,9 @@ * shape functions on simplices. * * The topology can be chosen at compile time by setting TOPOLOGY - * to a string like + * to a Dune::GeometryType like * \code - * Pyramid > > + * GeometryTypes::simplex(2) * \endcode * which generates a 2d simplex. If TOPOLOGY is not set, all * topologies up to 4d are tested. Note, this may lead to prolonged @@ -72,24 +74,25 @@ return ret; } -template +template bool test(unsigned int order, bool verbose = false) { - typedef Dune::LagrangeBasisFactory BasisFactory; - typedef Dune::LagrangeCoefficientsFactory< Dune::EquidistantPointSet, Topology::dimension,double > LagrangeCoefficientsFactory; + constexpr Dune::GeometryType geometry = geometryId; + typedef Dune::LagrangeBasisFactory BasisFactory; + typedef Dune::LagrangeCoefficientsFactory< Dune::EquidistantPointSet, geometry.dim(), double > LagrangeCoefficientsFactory; bool ret = true; for (unsigned int o = 0; o <= order; ++o) { - const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< Topology >( o ); + const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< geometry >( o ); if ( pointsPtr == 0) continue; - std::cout << "# Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl; + std::cout << "Testing " << geometry << " with order " << o << std::endl; - typename BasisFactory::Object &basis = *BasisFactory::template create(o); + typename BasisFactory::Object &basis = *BasisFactory::template create(o); ret |= test(basis,*pointsPtr,verbose); @@ -97,10 +100,10 @@ // derivatives in a human readabible form (aka LaTeX source) #ifdef TEST_OUTPUT_FUNCTIONS std::stringstream name; - name << "lagrange_" << Topology::name() << "_p" << o << ".basis"; + name << "lagrange_" << geometry << "_p" << o << ".basis"; std::ofstream out(name.str().c_str()); - Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField>(out,basis); - Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField>(out,basis); + Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis); + Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis); #endif // TEST_OUTPUT_FUNCTIONS LagrangeCoefficientsFactory::release( pointsPtr ); @@ -152,31 +155,29 @@ bool tests = true; #ifdef CHECKDIM1 - tests &= test > (order); - tests &= test > (order); + tests &= test (order); + tests &= test (order); #endif #ifdef CHECKDIM2 - tests &= test > > (order); - tests &= test > >(order); + tests &= test (order); + tests &= test (order); #endif #ifdef CHECKDIM3 - tests &= test > > >(order); - tests &= test > > >(order); - tests &= test > > >(order); + tests &= test (order); + tests &= test (order); + tests &= test (order); + tests &= test (order); #endif - // tests &= test > > >(order); - std::cout << "NOT CHECKING PYRAMID!" << std::endl; - // reduce tested order to 4 in 4d unless explicitly asked for more if (argc < 2) order = 4; #ifdef CHECKDIM4 - tests &= test > > > >(order); - tests &= test > > > >(order); + tests &= test (order); + tests &= test (order); #endif return (tests ? 0 : 1); diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-localfe.hh dune-localfunctions-2.8.0/dune/localfunctions/test/test-localfe.hh --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-localfe.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-localfe.hh 2021-08-31 14:03:38.000000000 +0000 @@ -30,15 +30,16 @@ // This class wraps one shape function of a local finite element as a function // that can be feed to the LocalInterpolation::interpolate method. template -class ShapeFunctionAsFunction : - // public Dune::LocalFiniteElementFunctionBase::type - public Dune::LocalFiniteElementFunctionBase::FunctionBase - // public Dune::LocalFiniteElementFunctionBase::VirtualFunctionBase +class ShapeFunctionAsFunction { public: typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType; typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType; - typedef typename Dune::Function Base; + + struct Traits { + typedef typename FE::Traits::LocalBasisType::Traits::DomainType DomainType; + typedef typename FE::Traits::LocalBasisType::Traits::RangeType RangeType; + }; typedef typename FE::Traits::LocalBasisType::Traits::RangeFieldType CT; @@ -98,7 +99,7 @@ { ////////////////////////////////////////////////////////////////////////////// // Part A: Feed the shape functions to the 'interpolate' method in form of - // a class derived from FunctionBase. + // a class providing an evaluate() method. // This way is deprecated since dune-localfunctions 2.7. ////////////////////////////////////////////////////////////////////////////// @@ -178,9 +179,63 @@ } +// Check whether the space spanned by the shape functions +// contains the constant functions +template +bool testCanRepresentConstants(const FE& fe, + unsigned order = 5) +{ + typedef typename FE::Traits::LocalBasisType LB; + using RangeType = typename LB::Traits::RangeType; + + bool success = true; + + // Construct the constant '1' function + auto constantOne = [](const typename LB::Traits::DomainType& xi) { return RangeType(1.0); }; + + // Project the constant function onto the FE space + std::vector coefficients; + fe.localInterpolation().interpolate(constantOne, coefficients); + + // A set of test points + const auto& quad = Dune::QuadratureRules::rule(fe.type(),order); + + // Loop over all quadrature points + for (size_t i=0; i values; + fe.localBasis().evaluateFunction(testPoint, values); + + RangeType sum(0); + for (size_t j=0; j TOL) + { + std::cout << "Finite element type " << Dune::className(fe) + << " cannot represent constant functions!" << std::endl; + std::cout << " At position: " << testPoint << "," + << std::endl; + std::cout << " discrete approximation of the '1' function has value " << sum + << std::endl; + std::cout << std::endl; + success = false; + } + + } // Loop over all quadrature points + + return success; +} + // check whether Jacobian agrees with FD approximation template -bool testJacobian(const FE& fe, unsigned order = 2) +bool testJacobian(const FE& fe, + unsigned order = 2, + const std::function derivativePointSkip = nullptr) { typedef typename FE::Traits::LocalBasisType LB; @@ -216,6 +271,10 @@ return false; } + // Skip this test point if we are supposed to + if (derivativePointSkip && derivativePointSkip(quad[i].position())) + continue; + // Loop over all directions for (int k=0; k static bool test(const FE& fe, - double eps, double delta, unsigned int diffOrder, std::size_t order = 2) + double eps, double delta, unsigned int diffOrder, std::size_t order = 2, + const std::function derivativePointSkip = nullptr) { bool success = true; @@ -282,10 +342,10 @@ std::cout << "No test for differentiability orders larger than 2!" << std::endl; if (diffOrder >= 2) - success = success and testOrder2(fe, eps, delta, order); + success = success and testOrder2(fe, eps, delta, order, derivativePointSkip); if (diffOrder >= 1) - success = success and testOrder1(fe, eps, delta, order); + success = success and testOrder1(fe, eps, delta, order, derivativePointSkip); success = success and testOrder0(fe, eps, delta, order); @@ -377,7 +437,8 @@ static bool testOrder1(const FE& fe, double eps, double delta, - std::size_t order = 2) + std::size_t order = 2, + const std::function derivativePointSkip = nullptr) { typedef typename FE::Traits::LocalBasisType LB; typedef typename LB::Traits::RangeFieldType RangeField; @@ -400,6 +461,10 @@ // Get a test point const Dune::FieldVector& testPoint = quad[i].position(); + // Skip the test points we are supposed to skip + if (derivativePointSkip && derivativePointSkip(testPoint)) + continue; + // Loop over all directions for (int k = 0; k < LB::Traits::dimDomain; k++) { @@ -474,7 +539,8 @@ static bool testOrder2(const FE& fe, double eps, double delta, - std::size_t order = 2) + std::size_t order = 2, + const std::function derivativePointSkip = nullptr) { typedef typename FE::Traits::LocalBasisType LocalBasis; typedef typename LocalBasis::Traits::DomainFieldType DF; @@ -502,6 +568,10 @@ // Get a test point const Domain& testPoint = quad[i].position(); + // Skip the test points we are supposed to skip + if (derivativePointSkip && derivativePointSkip(testPoint)) + continue; + // For testing the 'partial' method std::array >, dimR> partialHessians; for (size_t k = 0; k < dimR; k++) @@ -611,12 +681,25 @@ DisableLocalInterpolation = 1, DisableVirtualInterface = 2, DisableJacobian = 4, - DisableEvaluate = 8 + DisableEvaluate = 8, + DisableRepresentConstants = 16 }; -// call tests for given finite element +/** \brief Call tests for given finite element + * + * \param derivativePointSkip This is a small predicate class that allows to skip certain + * points when testing the derivative implementations. It exists because some + * finite elements are not everywhere differentiable, but we still want to run + * the tests for derivatives. Rather than constructing special sets of test + * points that avoid the problematic parts of the domain, we simply skip + * all test points that happen to be somewhere where the shape functions are + * not differentiable. + */ template -bool testFE(const FE& fe, char disabledTests = DisableNone, unsigned int diffOrder = 0) +bool testFE(const FE& fe, + char disabledTests = DisableNone, + unsigned int diffOrder = 0, + const std::function derivativePointSkip = nullptr) { // Order of the quadrature rule used to generate test points unsigned int quadOrder = 2; @@ -681,9 +764,15 @@ { success = testLocalInterpolation(fe) and success; } + + if (not (disabledTests & DisableRepresentConstants)) + { + success = testCanRepresentConstants(fe) and success; + } + if (not (disabledTests & DisableJacobian)) { - success = testJacobian(fe, quadOrder) and success; + success = testJacobian(fe, quadOrder, derivativePointSkip) and success; } else { @@ -693,7 +782,7 @@ if (not (disabledTests & DisableEvaluate)) { - success = TestPartial::test(fe, TOL, jacobianTOL, diffOrder, quadOrder) and success; + success = TestPartial::test(fe, TOL, jacobianTOL, diffOrder, quadOrder, derivativePointSkip) and success; } if (not (disabledTests & DisableVirtualInterface)) @@ -707,7 +796,7 @@ success = testLocalInterpolation(virtualFE) and success; if (not (disabledTests & DisableJacobian)) { - success = testJacobian(virtualFE) and success; + success = testJacobian(virtualFE, quadOrder, derivativePointSkip) and success; } else { @@ -722,5 +811,6 @@ #define TEST_FE(A) { bool b = testFE(A); std::cout << "testFE(" #A ") " << (b?"succeeded\n":"failed\n"); success &= b; } #define TEST_FE2(A,B) { bool b = testFE(A, B); std::cout << "testFE(" #A ", " #B ") " << (b?"succeeded\n":"failed\n"); success &= b; } #define TEST_FE3(A,B,C) { bool b = testFE(A, B, C); std::cout << "testFE(" #A ", " #B ", " #C ") " << (b?"succeeded\n":"failed\n"); success &= b; } +#define TEST_FE4(A,B,C,D) { bool b = testFE(A, B, C, D); std::cout << "testFE(" #A ", " #B ", " #C ", " #D ") " << (b?"succeeded\n":"failed\n"); success &= b; } #endif // DUNE_LOCALFUNCTIONS_TEST_TEST_LOCALFE_HH diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-nedelecsimplex.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-nedelecsimplex.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-nedelecsimplex.cc 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-nedelecsimplex.cc 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,120 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#include +#include +#include +#include + +/** + * \file + * \brief Performs some tests for the generic Nedelec + * shape functions on simplices. + * + * The topology can be chosen at compile time by setting TOPOLOGY + * to a Dune::GeometryType like + * \code + * GeometryTypes::simplex(2) + * \endcode + * which generates a 2d simplex. If TOPOLOGY is not set, triangles and tetrahedra + * are tested. Note, this may lead to prolonged compiler runs. + * + * For debugging purpuse the functions and the derivatives can be + * printed. You have to define the macro TEST_OUTPUT_FUNCTIONS to + * activate this function. + */ + +#if HAVE_GMP +typedef Dune::GMPField< 128 > StorageField; +typedef Dune::GMPField< 512 > ComputeField; +#else +typedef double StorageField; +typedef double ComputeField; +#endif + +template< Dune::GeometryType::Id geometryId > +bool test(unsigned int order) +{ + bool ret = true; + constexpr Dune::GeometryType geometry = geometryId; + + for (unsigned int o = 1; o <= order; ++o) + { + std::cout << "Testing " << geometry << " of the " << o <<"-th kind"<< std::endl; + typedef Dune::NedelecBasisFactory BasisFactory; + const typename BasisFactory::Object &basis = *BasisFactory::template create(o); + + // define the macro TEST_OUTPUT_FUNCTIONS to output files containing functions and + // derivatives in a human readabible form (aka LaTeX source) +#ifdef TEST_OUTPUT_FUNCTIONS + std::stringstream name; + name << "ned_" << geometry << "_p" << o << ".basis"; + std::ofstream out(name.str().c_str()); + Dune::basisPrint<0, BasisFactory, typename BasisFactory::StorageField, geometry>(out,basis); + Dune::basisPrint<1, BasisFactory, typename BasisFactory::StorageField, geometry>(out,basis); +#endif // TEST_OUTPUT_FUNCTIONS + + // test interpolation + typedef Dune::NedelecL2InterpolationFactory InterpolationFactory; + const typename InterpolationFactory::Object &interpol = *InterpolationFactory::template create(o); + Dune::LFEMatrix matrix; + interpol.interpolate(basis,matrix); + for (unsigned int i=0; i 1000.*Dune::Zero::epsilon() ) + std::cout << " non-zero entry in interpolation matrix: " + << "(" << i << "," << j << ") = " << Dune::field_cast(matrix(i,j)) + << std::endl; + + BasisFactory::release(&basis); + std::cout<<"----------------------------------------------------------------------------------------------------------------\n"; + } + if (!ret) { + std::cout << " FAILED !" << std::endl; + } + return ret; +} + +#ifdef CHECKDIM + #if CHECKDIM==2 + #define CHECKDIM2 + #elif CHECKDIM==3 + #define CHECKDIM3 + #endif +#else + #define CHECKDIM2 + #define CHECKDIM3 +#endif + + + +int main ( int argc, char **argv ) +{ + using namespace Dune; + using namespace Impl; + + unsigned int order = (argc < 2) ? 5 : atoi(argv[1]); + + if (argc < 2) + { + std::cerr << "Usage: " << argv[ 0 ] << "

" << std::endl + << "Using default kind of " << order << std::endl; + } +#ifdef TOPOLOGY + return (test(order) ? 0 : 1 ); +#else + bool tests = true; + + +#ifdef CHECKDIM2 + tests &= test(order); +#endif + +#ifdef CHECKDIM3 + tests &= test(order); +#endif + + return (tests ? 0 : 1); +#endif // TOPOLOGY +} diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-orthonormal.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-orthonormal.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-orthonormal.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-orthonormal.cc 2021-08-31 14:03:38.000000000 +0000 @@ -15,9 +15,9 @@ * shape functions on simplices. * * The topology can be chosen at compile time by setting TOPOLOGY - * to a string like + * to a Dune::GeometryType like * \code - * Pyramid > > + * GeometryTypes::simplex(2) * \endcode * which generates a 2d simplex. If TOPOLOGY is not set, all * topologies up to 4d are tested. Note, this may lead to prolonged @@ -36,16 +36,16 @@ typedef double ComputeField; #endif -template +template< Dune::GeometryType::Id geometryId > bool test(unsigned int order) { bool ret = true; - Dune::GeometryType gt(Topology::id, Topology::dimension); + constexpr Dune::GeometryType geometry = geometryId; for (unsigned int o = 0; o <= order; ++o) { - std::cout << "Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl; - typedef Dune::OrthonormalBasisFactory BasisFactory; - const typename BasisFactory::Object &basis = *BasisFactory::template create(o); + std::cout << "Testing " << geometry << " with order " << o << std::endl; + typedef Dune::OrthonormalBasisFactory BasisFactory; + const typename BasisFactory::Object &basis = *BasisFactory::template create(o); const unsigned int size = basis.size( ); @@ -55,8 +55,8 @@ for( unsigned int i = 0; i < size * size; ++i ) m[ i ] = 0; - const Dune::QuadratureRule &quadrature = - Dune::QuadratureRules::rule(gt,2*order+1); + const Dune::QuadratureRule &quadrature = + Dune::QuadratureRules::rule(geometry,2*order+1); const unsigned int quadratureSize = quadrature.size(); for( unsigned int qi = 0; qi < quadratureSize; ++qi ) { @@ -84,10 +84,10 @@ // derivatives in a human readabible form (aka LaTeX source) #ifdef TEST_OUTPUT_FUNCTIONS std::stringstream name; - name << "orthonormal_" << Topology::name() << "_p" << o << ".basis"; + name << "orthonormal_" << geometry << "_p" << o << ".basis"; std::ofstream out(name.str().c_str()); - Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField>(out,basis); - Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField>(out,basis); + Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis); + Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis); #endif // TEST_OUTPUT_FUNCTIONS BasisFactory::release(&basis); @@ -136,25 +136,25 @@ bool tests = true; #ifdef CHECKDIM1 - tests &= test > (order); - tests &= test > (order); + tests &= test (order); + tests &= test (order); #endif #ifdef CHECKDIM2 - tests &= test > > (order); - tests &= test > >(order); + tests &= test (order); + tests &= test (order); #endif #ifdef CHECKDIM3 - tests &= test > > >(order); - tests &= test > > >(order); - tests &= test > > >(order); - tests &= test > > >(order); + tests &= test (order); + tests &= test (order); + tests &= test (order); + tests &= test (order); #endif #ifdef CHECKDIM4 - tests &= test > > > >(order); - tests &= test > > > >(order); + tests &= test (order); + tests &= test (order); #endif return (tests ? 0 : 1); diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-pk2d.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-pk2d.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-pk2d.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-pk2d.cc 2021-08-31 14:03:38.000000000 +0000 @@ -8,11 +8,11 @@ #include #include #include +#include #include #include #include -#include #include #include @@ -59,7 +59,7 @@ int result = 77; static constexpr std::size_t max_k = 20; - Dune::Hybrid::forEach(Dune::Std::make_index_sequence{},[&](auto i){test(result);}); + Dune::Hybrid::forEach(std::make_index_sequence{},[&](auto i){test(result);}); return result; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-power-monomial.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-power-monomial.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-power-monomial.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-power-monomial.cc 2021-08-31 14:03:38.000000000 +0000 @@ -8,10 +8,10 @@ #include #include #include +#include #include #include -#include #include #include @@ -64,20 +64,20 @@ template static void DimR(int &result) { - Dune::Hybrid::forEach(Dune::Std::make_index_sequence<4>{},[&](auto i){Order(result);}); + Dune::Hybrid::forEach(std::make_index_sequence<4>{},[&](auto i){Order(result);}); } template static void DimD(int &result) { - Dune::Hybrid::forEach(Dune::Std::make_index_sequence<4>{},[&](auto i){DimR(result);}); + Dune::Hybrid::forEach(std::make_index_sequence<4>{},[&](auto i){DimR(result);}); } int main(int argc, char** argv) { try { int result = 77; - Dune::Hybrid::forEach(Dune::Std::make_index_sequence<3>{},[&](auto i){DimD(result);}); + Dune::Hybrid::forEach(std::make_index_sequence<3>{},[&](auto i){DimD(result);}); return result; } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/test-raviartthomassimplex.cc dune-localfunctions-2.8.0/dune/localfunctions/test/test-raviartthomassimplex.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/test-raviartthomassimplex.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/test-raviartthomassimplex.cc 2021-08-31 14:03:38.000000000 +0000 @@ -11,9 +11,9 @@ * shape functions on simplices. * * The topology can be chosen at compile time by setting TOPOLOGY - * to a string like + * to a Dune::GeometryType like * \code - * Pyramid > > + * GeometryTypes::simplex(2) * \endcode * which generates a 2d simplex. If TOPOLOGY is not set, all * topologies up to 4d are tested. Note, this may lead to prolonged @@ -32,30 +32,31 @@ typedef double ComputeField; #endif -template +template< Dune::GeometryType::Id geometryId > bool test(unsigned int order) { bool ret = true; + constexpr Dune::GeometryType geometry = geometryId; for (unsigned int o = 0; o <= order; ++o) { - std::cout << "Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl; - typedef Dune::RaviartThomasBasisFactory BasisFactory; - const typename BasisFactory::Object &basis = *BasisFactory::template create(o); + std::cout << "Testing " << geometry << " with order " << o << std::endl; + typedef Dune::RaviartThomasBasisFactory BasisFactory; + const typename BasisFactory::Object &basis = *BasisFactory::template create(o); // define the macro TEST_OUTPUT_FUNCTIONS to output files containing functions and // derivatives in a human readabible form (aka LaTeX source) #ifdef TEST_OUTPUT_FUNCTIONS std::stringstream name; - name << "rt_" << Topology::name() << "_p" << o << ".basis"; + name << "rt_" << geometry << "_p" << o << ".basis"; std::ofstream out(name.str().c_str()); - Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField>(out,basis); - Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField>(out,basis); + Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis); + Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis); #endif // TEST_OUTPUT_FUNCTIONS // test interpolation - typedef Dune::RaviartThomasL2InterpolationFactory InterpolationFactory; - const typename InterpolationFactory::Object &interpol = *InterpolationFactory::template create(o); + typedef Dune::RaviartThomasL2InterpolationFactory InterpolationFactory; + const typename InterpolationFactory::Object &interpol = *InterpolationFactory::template create(o); Dune::LFEMatrix matrix; interpol.interpolate(basis,matrix); for (unsigned int i=0; i > (order); + tests &= test (order); #endif #ifdef CHECKDIM2 - tests &= test > >(order); + tests &= test (order); #endif #ifdef CHECKDIM3 - tests &= test > > >(order); + tests &= test (order); #endif // reduce tested order to 4 in 4d unless explicitly asked for more @@ -128,7 +129,7 @@ order = 4; #ifdef CHECKDIM4 - tests &= test > > > >(order); + tests &= test (order); #endif return (tests ? 0 : 1); diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/test/virtualshapefunctiontest.cc dune-localfunctions-2.8.0/dune/localfunctions/test/virtualshapefunctiontest.cc --- dune-localfunctions-2.7.1/dune/localfunctions/test/virtualshapefunctiontest.cc 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/test/virtualshapefunctiontest.cc 2021-08-31 14:03:38.000000000 +0000 @@ -32,11 +32,17 @@ {} // A test function to test the local interpolation -template +template struct TestFunction -// : public VirtualFunction - : public Function { + using DomainType = D; + using RangeType = R; + + struct Traits { + using DomainType = D; + using RangeType = R; + }; + void evaluate(const DomainType& in, RangeType& out) const { // May not be flexible enough to compile for all range types out = 1; @@ -49,7 +55,7 @@ { // call each method once to test that it's there syntax_check( localBasis->order() ); - unsigned int size DUNE_UNUSED = localBasis->size(); + [[maybe_unused]] unsigned int size = localBasis->size(); // evaluate the local basis at (0,...,0) typename T::DomainType in(0); @@ -87,7 +93,24 @@ // Test interpolation of a function object derived from VirtualFunction TestFunction testFunction; std::vector coefficients; + + ////////////////////////////////////////////////////////////////////////////// + // Part A: Feed the function to the 'interpolate' method in form of + // a class providing an evaluate() method. + // This way is deprecated since dune-localfunctions 2.7. + ////////////////////////////////////////////////////////////////////////////// localInterpolation->interpolate(testFunction, coefficients); + + ////////////////////////////////////////////////////////////////////////////// + // Part B: Redo the same test, but feed the function to the + // 'interpolate' method in form of a callable. + ////////////////////////////////////////////////////////////////////////////// + auto callableTestFunction = [&](const auto& x) { + RangeType y(0); + testFunction.evaluate(x,y); + return y; + }; + localInterpolation->interpolate(callableTestFunction, coefficients); } // Test all methods of a local finite element given as a pointer to the abstract base class diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/basismatrix.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/basismatrix.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/basismatrix.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/basismatrix.hh 2021-08-31 14:03:38.000000000 +0000 @@ -53,13 +53,13 @@ unsigned int cols_; }; - template< class Topology, class F, + template< GeometryType::Id geometryId, class F, class Interpolation, class Field > - struct BasisMatrix< const MonomialBasis< Topology, F >, Interpolation, Field > - : public BasisMatrixBase< const MonomialBasis< Topology, F >, Interpolation, Field > + struct BasisMatrix< const MonomialBasis< geometryId, F >, Interpolation, Field > + : public BasisMatrixBase< const MonomialBasis< geometryId, F >, Interpolation, Field > { - typedef const MonomialBasis< Topology, F > PreBasis; + typedef const MonomialBasis< geometryId, F > PreBasis; typedef BasisMatrixBase Base; typedef typename Base::Matrix Matrix; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/basisprint.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/basisprint.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/basisprint.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/basisprint.hh 2021-08-31 14:03:38.000000000 +0000 @@ -17,7 +17,7 @@ **********************************************/ // default argument does not work for gcc 4.1.2 // template - template + template void basisPrint(std::ostream &out, typename BasisFactory::Object &basis) { @@ -31,7 +31,7 @@ typedef typename Basis::CoefficientMatrix CMatrix; typedef PolynomialBasis, CMatrix > PrintBasis; - MIBasis *miBasis = MIBasisFactory::create( Dune::GeometryType( basis.basis().topologyId(),dimension ),basis.basis().order()); + MIBasis *miBasis = MIBasisFactory::template create( basis.basis().order()); PrintBasis printBasis(*miBasis,basis.matrix(),basis.size()); unsigned int size = printBasis.size(); @@ -39,11 +39,6 @@ out << "% Number of base functions: " << size << std::endl; out << "% Derivative order: " << deriv << std::endl; - /* - std::vector< FieldVector< - LFETensor,PrintBasis::dimRange> > - y( size ); - */ std::vector< FieldVector< FieldVector::size>, PrintBasis::dimRange> > y( size ); @@ -62,8 +57,8 @@ } MIBasisFactory::release(miBasis); } - // template - template + + template void basisPrint(std::ostream &out, typename BasisFactory::Key &key) { diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/defaultbasisfactory.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/defaultbasisfactory.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/defaultbasisfactory.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/defaultbasisfactory.hh 2021-08-31 14:03:38.000000000 +0000 @@ -56,15 +56,15 @@ Type; }; - template< class Topology > + template< GeometryType::Id geometryId > static Object *create ( const Key &key ) { const typename PreBasisFactory::Key preBasisKey = PreBasisKeyExtractor::apply(key); - const PreBasis *preBasis = PreBasisFactory::template create( preBasisKey ); - const Interpolation *interpol = InterpolationFactory::template create( key ); + const PreBasis *preBasis = PreBasisFactory::template create( preBasisKey ); + const Interpolation *interpol = InterpolationFactory::template create( key ); BasisMatrix< PreBasis, Interpolation, ComputeField > matrix( *preBasis, *interpol ); - const MonomialBasis *monomialBasis = MonomialBasisFactory::template create< Topology >( preBasis->order() ); + const MonomialBasis *monomialBasis = MonomialBasisFactory::template create< geometryId >( preBasis->order() ); Basis *basis = new Basis( *monomialBasis ); diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/dglocalcoefficients.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/dglocalcoefficients.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/dglocalcoefficients.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/dglocalcoefficients.hh 2021-08-31 14:03:38.000000000 +0000 @@ -59,11 +59,11 @@ typedef typename BasisFactory::Key Key; typedef const DGLocalCoefficients Object; - template< class Topology > + template< GeometryType::Id geometryId > static Object *create ( const Key &key ) { const typename BasisFactory::Object *basis - = BasisFactory::template create< Topology >( key ); + = BasisFactory::template create< geometryId >( key ); Object *coefficients = new Object( basis->size() ); BasisFactory::release( basis ); return coefficients; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/l2interpolation.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/l2interpolation.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/l2interpolation.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/l2interpolation.hh 2021-08-31 14:03:38.000000000 +0000 @@ -204,12 +204,12 @@ typedef LocalL2Interpolation< Basis, Quadrature, onb > LocalInterpolation; typedef const LocalInterpolation Object; - template< class Topology > + template< GeometryType::Id geometryId > static Object *create ( const Key &key ) { - Dune::GeometryType gt(Topology::id, Topology::dimension); - const Basis *basis = BasisFactory::template create< Topology >( key ); - const Quadrature & quadrature = QuadratureProvider::rule(gt, 2*basis->order()+1); + constexpr Dune::GeometryType geometry = geometryId; + const Basis *basis = BasisFactory::template create< geometry >( key ); + const Quadrature & quadrature = QuadratureProvider::rule(geometry, 2*basis->order()+1); return new Object( *basis, quadrature ); } static void release ( Object *object ) diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/lfematrix.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/lfematrix.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/lfematrix.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/lfematrix.hh 2021-08-31 14:03:38.000000000 +0000 @@ -86,6 +86,7 @@ bool invert () { + using std::abs; assert( rows() == cols() ); std::vector p(rows()); for (unsigned int j=0; j max ) + if ( abs( (*this)(i,j) ) > max ) { - max = std::abs( (*this)(i,j) ); + max = abs( (*this)(i,j) ); r = i; } } diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/localfiniteelement.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/localfiniteelement.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/localfiniteelement.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/localfiniteelement.hh 2021-08-31 14:03:38.000000000 +0000 @@ -4,6 +4,7 @@ #define DUNE_GENERIC_LOCALFINITEELEMENT_HH #include +#include #include #include @@ -39,20 +40,24 @@ /** \todo Please doc me */ GenericLocalFiniteElement ( const GeometryType >, const Key &key ) - : topologyId_( gt.id() ), + : geometry_( gt ), key_( key ), finiteElement_() { - Impl::IfTopology< FiniteElement::template Maker, dimDomain >::apply( topologyId_, key_, finiteElement_ ); + Impl::toGeometryTypeIdConstant(type(), [&](auto geometryTypeId) { + finiteElement_.template create(key_); + }); } /** \todo Please doc me */ GenericLocalFiniteElement ( const GenericLocalFiniteElement &other ) - : topologyId_( other.topologyId_ ), + : geometry_( other.type() ), key_( other.key_ ), finiteElement_() { - Impl::IfTopology< FiniteElement::template Maker, dimDomain >::apply( topologyId_, key_, finiteElement_ ); + Impl::toGeometryTypeIdConstant(type(), [&](auto geometryTypeId) { + finiteElement_.template create(key_); + }); } ~GenericLocalFiniteElement() @@ -91,28 +96,20 @@ */ GeometryType type () const { - return GeometryType(topologyId_,dimDomain); - } - - /** \todo Please doc me ! - * \deprecated Deprecated in dune-localfunctions 2.7. Use type().id() instead! - */ - DUNE_DEPRECATED_MSG("Use type().id() instead!") - unsigned int topologyId () const - { - return topologyId_; + return geometry_; } private: struct FiniteElement { FiniteElement() : basis_(0), coeff_(0), interpol_(0) {} - template + + template < GeometryType::Id geometryId > void create( const Key &key ) { release(); - basis_ = BasisF::template create(key); - coeff_ = CoeffF::template create(key); - interpol_ = InterpolF::template create(key); + basis_ = BasisF::template create(key); + coeff_ = CoeffF::template create(key); + interpol_ = InterpolF::template create(key); } void release() { @@ -126,19 +123,11 @@ coeff_=0; interpol_=0; } - template< class Topology > - struct Maker - { - static void apply ( const Key &key, FiniteElement &finiteElement ) - { - finiteElement.template create(key); - }; - }; typename Traits::LocalBasisType *basis_; typename Traits::LocalCoefficientsType *coeff_; typename Traits::LocalInterpolationType *interpol_; }; - unsigned int topologyId_; + GeometryType geometry_; Key key_; FiniteElement finiteElement_; }; diff -Nru dune-localfunctions-2.7.1/dune/localfunctions/utility/monomialbasis.hh dune-localfunctions-2.8.0/dune/localfunctions/utility/monomialbasis.hh --- dune-localfunctions-2.7.1/dune/localfunctions/utility/monomialbasis.hh 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/localfunctions/utility/monomialbasis.hh 2021-08-31 14:03:38.000000000 +0000 @@ -8,8 +8,8 @@ #include #include -#include #include +#include #include #include @@ -41,12 +41,12 @@ * build over a non simplex base geometry. * * Central classes: - * 1) template< class Topology, class F > + * 1) template< GeometryType::Id geometryId, class F > * class MonomialBasisImpl; * Implementation of the monomial evaluation for * a given topology and field type. * The method evaluate fills a F* vector - * 2) template< class Topology, class F > + * 2) template< GeometryType::Id geometryId, class F > * class MonomialBasis * The base class for the static monomial evaluation * providing addiional evaluate methods including @@ -66,10 +66,10 @@ // Internal Forward Declarations // ----------------------------- - template< class Topology > + template< GeometryType::Id geometryId > class MonomialBasisSize; - template< class Topology, class F > + template< GeometryType::Id geometryId, class F > class MonomialBasis; @@ -77,10 +77,10 @@ // MonomialBasisSize // ----------------- - template< class TopologyType > + template< GeometryType::Id geometryId > class MonomialBasisSize { - typedef MonomialBasisSize< TopologyType > This; + typedef MonomialBasisSize< geometryId > This; public: static This &instance () @@ -133,7 +133,8 @@ sizes_ = new unsigned int[ order+1 ]; numBaseFunctions_ = new unsigned int[ order+1 ]; - constexpr auto dim = TopologyType::dimension; + constexpr GeometryType geometry = geometryId; + constexpr auto dim = geometry.dim(); sizes_[ 0 ] = 1; for( unsigned int k = 1; k <= order; ++k ) @@ -143,7 +144,7 @@ for( int codim=dim-1; codim>=0; codim--) { - if (Impl::isPrism(TopologyType::id,dim,codim)) + if (Impl::isPrism(geometry.id(),dim,codim)) { for( unsigned int k = 1; k <= order; ++k ) { @@ -172,8 +173,8 @@ template< int mydim, int dim, class F > struct MonomialBasisHelper { - typedef MonomialBasisSize< typename Impl::SimplexTopology< mydim >::type > MySize; - typedef MonomialBasisSize< typename Impl::SimplexTopology< dim >::type > Size; + typedef MonomialBasisSize< GeometryTypes::simplex(mydim).toId() > MySize; + typedef MonomialBasisSize< GeometryTypes::simplex(dim).toId() > Size; static void copy ( const unsigned int deriv, F *&wit, F *&rit, const unsigned int numBaseFunctions, const F &z ) @@ -226,26 +227,23 @@ // MonomialBasisImpl // ----------------- - template< class Topology, class F > - class MonomialBasisImpl; - - template< class F > - class MonomialBasisImpl< Impl::Point, F > + template< GeometryType::Id geometryId, class F> + class MonomialBasisImpl { - typedef MonomialBasisImpl< Impl::Point, F > This; - public: - typedef Impl::Point Topology; typedef F Field; - static const unsigned int dimDomain = Topology::dimension; + static constexpr GeometryType geometry = geometryId; + + static const unsigned int dimDomain = geometry.dim(); typedef FieldVector< Field, dimDomain > DomainVector; private: - friend class MonomialBasis< Topology, Field >; - friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >; - friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >; + friend class MonomialBasis< geometryId, Field >; + + MonomialBasisImpl () + {} template< int dimD > void evaluate ( const unsigned int deriv, const unsigned int order, @@ -253,82 +251,120 @@ const unsigned int block, const unsigned int *const offsets, Field *const values ) const { + //start with vertex *values = Unity< F >(); F *const end = values + block; for( Field *it = values+1 ; it != end; ++it ) *it = Zero< F >(); - } - void integrate ( const unsigned int order, - const unsigned int *const offsets, - Field *const values ) const - { - values[ 0 ] = Unity< Field >(); - } - }; - - template< class BaseTopology, class F > - class MonomialBasisImpl< Impl::Prism< BaseTopology >, F > - { - typedef MonomialBasisImpl< Impl::Prism< BaseTopology >, F > This; + constexpr GeometryType gt = GeometryTypes::vertex; - public: - typedef Impl::Prism< BaseTopology > Topology; - typedef F Field; - - static const unsigned int dimDomain = Topology::dimension; - - typedef FieldVector< Field, dimDomain > DomainVector; - - private: - friend class MonomialBasis< Topology, Field >; - friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >; - friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >; - - typedef MonomialBasisSize< BaseTopology > BaseSize; - typedef MonomialBasisSize< Topology > Size; - - MonomialBasisImpl< BaseTopology, Field > baseBasis_; - - MonomialBasisImpl () - {} + if constexpr ( geometry == gt) + return; + else + evaluate(deriv, order, x, block, offsets, values ); + } - template< int dimD > + template void evaluate ( const unsigned int deriv, const unsigned int order, const FieldVector< Field, dimD > &x, const unsigned int block, const unsigned int *const offsets, Field *const values ) const { - typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper; + + static constexpr GeometryType baseGeometry = baseGeometryId; + + auto constexpr isPrismatic = geometry.isPrismatic(baseGeometry.dim()); + + // compute + typedef MonomialBasisHelper< baseGeometry.dim() + 1, dimD, Field > Helper; + typedef MonomialBasisSize BaseSize; + const BaseSize &size = BaseSize::instance(); const_cast(size).computeSizes(order); - const Field &z = x[ dimDomain-1 ]; - - // fill first column - baseBasis_.evaluate( deriv, order, x, block, offsets, values ); + const Field &z = x[ baseGeometry.dim() ]; Field *row0 = values; for( unsigned int k = 1; k <= order; ++k ) { Field *row1 = values + block*offsets[ k-1 ]; Field *wit = row1 + block*size.sizes_[ k ]; - Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z ); + if constexpr ( isPrismatic ) + Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z ); Helper::copy( deriv, wit, row0, size( k-1 ), z ); row0 = row1; } + + // stop if desired dimension is reached + if constexpr( baseGeometry.dim() == dimDomain-1) + return; + else + { + constexpr GeometryType nextGeometry = isPrismatic ? GeometryTypes::prismaticExtension(baseGeometry) + : GeometryTypes::conicalExtension(baseGeometry); + + evaluate(deriv, order, x, block, offsets, values ); + } } void integrate ( const unsigned int order, const unsigned int *const offsets, Field *const values ) const { - const BaseSize &size = BaseSize::instance(); - const Size &mySize = Size::instance(); - // fill first column - baseBasis_.integrate( order, offsets, values ); + //start with vertex + values[ 0 ] = Unity< Field >(); + static constexpr GeometryType gt = GeometryTypes::vertex; + + if constexpr ( geometry == gt) + return; + else + integrate(order, offsets, values); + } + + template + void integrate ( const unsigned int order, + const unsigned int *const offsets, + Field *const values) const + { + static constexpr GeometryType baseGeometry = baseGeometryId; + + auto constexpr isPrismatic = geometry.isPrismatic(baseGeometry.dim()); + + // decide which kind of integration should be performed + if constexpr ( isPrismatic ) + integratePrismatic(order, offsets, values); + else + integrateConical(order, offsets, values); + + // stop if the desired dimension is reached + if constexpr( baseGeometry.dim() == dimDomain-1) + return; + else + { + static constexpr GeometryType nextGeometry = (isPrismatic ? GeometryTypes::prismaticExtension(baseGeometry) + : GeometryTypes::conicalExtension(baseGeometry)); + + integrate(order, offsets, values); + } + + } + + template + void integratePrismatic ( const unsigned int order, + const unsigned int *const offsets, + Field *const values ) const + { + typedef MonomialBasisSize BaseSize; + static const BaseSize &size = BaseSize::instance(); const unsigned int *const baseSizes = size.sizes_; + static constexpr GeometryType baseGeometry = baseGeometryId; + static constexpr GeometryType nextGeometry = GeometryTypes::prismaticExtension(baseGeometry); + + typedef MonomialBasisSize Size; + static const Size &mySize = Size::instance(); + Field *row0 = values; for( unsigned int k = 1; k <= order; ++k ) { @@ -350,133 +386,29 @@ row0 = row1; } } - }; - - template< class BaseTopology, class F > - class MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F > - { - typedef MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F > This; - - public: - typedef Impl::Pyramid< BaseTopology > Topology; - typedef F Field; - - static const unsigned int dimDomain = Topology::dimension; - - typedef FieldVector< Field, dimDomain > DomainVector; - private: - friend class MonomialBasis< Topology, Field >; - friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >; - friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >; - - typedef MonomialBasisSize< BaseTopology > BaseSize; - typedef MonomialBasisSize< Topology > Size; - - MonomialBasisImpl< BaseTopology, Field > baseBasis_; - MonomialBasisImpl () - {} - - template< int dimD > - void evaluateSimplexBase ( const unsigned int deriv, const unsigned int order, - const FieldVector< Field, dimD > &x, - const unsigned int block, const unsigned int *const offsets, - Field *const values, - const BaseSize &size ) const + template + void integrateConical ( const unsigned int order, + const unsigned int *const offsets, + Field *const values) const { - baseBasis_.evaluate( deriv, order, x, block, offsets, values ); - } - - template< int dimD > - void evaluatePyramidBase ( const unsigned int deriv, const unsigned int order, - const FieldVector< Field, dimD > &x, - const unsigned int block, const unsigned int *const offsets, - Field *const values, - const BaseSize &size ) const - { - Field omz = Unity< Field >() - x[ dimDomain-1 ]; - - if( Zero< Field >() < omz ) - { - const Field invomz = Unity< Field >() / omz; - FieldVector< Field, dimD > y; - for( unsigned int i = 0; i < dimDomain-1; ++i ) - y[ i ] = x[ i ] * invomz; - - // fill first column - baseBasis_.evaluate( deriv, order, y, block, offsets, values ); - - Field omzk = omz; - for( unsigned int k = 1; k <= order; ++k ) - { - Field *it = values + block*offsets[ k-1 ]; - Field *const end = it + block*size.sizes_[ k ]; - for( ; it != end; ++it ) - *it *= omzk; - omzk *= omz; - } - } - else - { - assert( deriv==0 ); - *values = Unity< Field >(); - for( unsigned int k = 1; k <= order; ++k ) - { - Field *it = values + block*offsets[ k-1 ]; - Field *const end = it + block*size.sizes_[ k ]; - for( ; it != end; ++it ) - *it = Zero< Field >(); - } - } - } - - template< int dimD > - void evaluate ( const unsigned int deriv, const unsigned int order, - const FieldVector< Field, dimD > &x, - const unsigned int block, const unsigned int *const offsets, - Field *const values ) const - { - typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper; - const BaseSize &size = BaseSize::instance(); - const_cast(size).computeSizes(order); - - if( Impl::IsSimplex< Topology >::value ) - evaluateSimplexBase( deriv, order, x, block, offsets, values, size ); - else - evaluatePyramidBase( deriv, order, x, block, offsets, values, size ); - - Field *row0 = values; - for( unsigned int k = 1; k <= order; ++k ) - { - Field *row1 = values + block*offsets[ k-1 ]; - Field *wit = row1 + block*size.sizes_[ k ]; - Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] ); - row0 = row1; - } - } - - void integrate ( const unsigned int order, - const unsigned int *const offsets, - Field *const values ) const - { - const BaseSize &size = BaseSize::instance(); - - // fill first column - baseBasis_.integrate( order, offsets, values ); - + typedef MonomialBasisSize BaseSize; + static const BaseSize &size = BaseSize::instance(); const unsigned int *const baseSizes = size.sizes_; + static constexpr GeometryType baseGeometry = baseGeometryId; + { Field *const col0End = values + baseSizes[ 0 ]; for( Field *it = values; it != col0End; ++it ) - *it *= Field( 1 ) / Field( int(dimDomain) ); + *it *= Field( 1 ) / Field( int(baseGeometry.dim()+1) ); } Field *row0 = values; for( unsigned int k = 1; k <= order; ++k ) { - const Field factor = (Field( 1 ) / Field( k + dimDomain )); + const Field factor = (Field( 1 ) / Field( k + baseGeometry.dim()+1)); Field *const row1 = values+offsets[ k-1 ]; Field *const col0End = row1 + baseSizes[ k ]; @@ -493,19 +425,20 @@ row0 = row1; } } - }; + }; // MonomialBasis // ------------- - template< class Topology, class F > + template< GeometryType::Id geometryId, class F > class MonomialBasis - : public MonomialBasisImpl< Topology, F > + : public MonomialBasisImpl< geometryId, F > { - typedef MonomialBasis< Topology, F > This; - typedef MonomialBasisImpl< Topology, F > Base; + static constexpr GeometryType geometry = geometryId; + typedef MonomialBasis< geometryId, F > This; + typedef MonomialBasisImpl< geometryId, F > Base; public: static const unsigned int dimension = Base::dimDomain; @@ -517,7 +450,7 @@ typedef Dune::FieldVector RangeVector; - typedef MonomialBasisSize Size; + typedef MonomialBasisSize Size; MonomialBasis (unsigned int order) : Base(), @@ -546,9 +479,8 @@ unsigned int derivSize ( const unsigned int deriv ) const { - typedef typename Impl::SimplexTopology< dimension >::type SimplexTopology; - MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv ); - return MonomialBasisSize< SimplexTopology >::instance() ( deriv ); + MonomialBasisSize< GeometryTypes::simplex(dimension).toId() >::instance().computeSizes( deriv ); + return MonomialBasisSize< GeometryTypes::simplex(dimension).toId() >::instance() ( deriv ); } unsigned int order () const @@ -558,7 +490,7 @@ unsigned int topologyId ( ) const { - return Topology::id; + return geometry.id(); } void evaluate ( const unsigned int deriv, const DomainVector &x, @@ -633,13 +565,13 @@ template< int dim,class F > class StandardMonomialBasis - : public MonomialBasis< typename Impl::SimplexTopology< dim >::type, F > + : public MonomialBasis< GeometryTypes::simplex(dim).toId() , F > { typedef StandardMonomialBasis< dim, F > This; - typedef MonomialBasis< typename Impl::SimplexTopology< dim >::type, F > Base; + typedef MonomialBasis< GeometryTypes::simplex(dim).toId(), F > Base; public: - typedef typename Impl::SimplexTopology< dim >::type Topology; + static constexpr GeometryType geometry = GeometryTypes::simplex(dim); static const int dimension = dim; StandardMonomialBasis ( unsigned int order ) @@ -654,13 +586,13 @@ template< int dim, class F > class StandardBiMonomialBasis - : public MonomialBasis< typename Impl::CubeTopology< dim >::type, F > + : public MonomialBasis< GeometryTypes::cube(dim).toId() , F > { typedef StandardBiMonomialBasis< dim, F > This; - typedef MonomialBasis< typename Impl::CubeTopology< dim >::type, F > Base; + typedef MonomialBasis< GeometryTypes::cube(dim).toId() , F > Base; public: - typedef typename Impl::CubeTopology< dim >::type Topology; + static constexpr GeometryType geometry = GeometryTypes::cube(dim); static const int dimension = dim; StandardBiMonomialBasis ( unsigned int order ) @@ -687,9 +619,14 @@ typedef FieldVector DomainVector; typedef FieldVector RangeVector; + [[deprecated("Use VirtualMonomialBasis(GeometryType gt, unsigned int order) instead.")]] explicit VirtualMonomialBasis(unsigned int topologyId, unsigned int order) - : order_(order), topologyId_(topologyId) {} + : order_(order), geometry_(GeometryType(topologyId,dim)) {} + + explicit VirtualMonomialBasis(const GeometryType& gt, + unsigned int order) + : order_(order), geometry_(gt) {} virtual ~VirtualMonomialBasis() {} @@ -705,9 +642,15 @@ return order_; } + [[deprecated("Use type().id() instead.")]] unsigned int topologyId ( ) const { - return topologyId_; + return type().id(); + } + + GeometryType type() const + { + return geometry_; } virtual void evaluate ( const unsigned int deriv, const DomainVector &x, @@ -769,22 +712,23 @@ } protected: unsigned int order_; - unsigned int topologyId_; + GeometryType geometry_; }; - template< class Topology, class F > + template< GeometryType::Id geometryId, class F > class VirtualMonomialBasisImpl - : public VirtualMonomialBasis< Topology::dimension, F > + : public VirtualMonomialBasis< static_cast(geometryId).dim(), F > { - typedef VirtualMonomialBasis< Topology::dimension, F > Base; - typedef VirtualMonomialBasisImpl< Topology, F > This; + static constexpr GeometryType geometry = geometryId; + typedef VirtualMonomialBasis< geometry.dim(), F > Base; + typedef VirtualMonomialBasisImpl< geometryId, F > This; public: typedef typename Base::Field Field; typedef typename Base::DomainVector DomainVector; VirtualMonomialBasisImpl(unsigned int order) - : Base(Topology::id,order), basis_(order) + : Base(geometry,order), basis_(order) {} const unsigned int *sizes ( ) const @@ -804,7 +748,7 @@ } private: - MonomialBasis basis_; + MonomialBasis basis_; using Base::order_; }; @@ -826,10 +770,10 @@ typedef MonomialBasisFactory Type; }; - template< class Topology > + template< GeometryType::Id geometryId > static Object* create ( const Key &order ) { - return new VirtualMonomialBasisImpl< Topology, StorageField >( order ); + return new VirtualMonomialBasisImpl< geometryId, StorageField >( order ); } static void release( Object *object ) { delete object; } }; diff -Nru dune-localfunctions-2.7.1/dune/python/CMakeLists.txt dune-localfunctions-2.8.0/dune/python/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/python/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/python/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1 @@ +add_subdirectory(localfunctions) diff -Nru dune-localfunctions-2.7.1/dune/python/localfunctions/CMakeLists.txt dune-localfunctions-2.8.0/dune/python/localfunctions/CMakeLists.txt --- dune-localfunctions-2.7.1/dune/python/localfunctions/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/python/localfunctions/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,4 @@ +install(FILES + localfiniteelement.hh + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/python/localfunctions +) diff -Nru dune-localfunctions-2.7.1/dune/python/localfunctions/localfiniteelement.hh dune-localfunctions-2.8.0/dune/python/localfunctions/localfiniteelement.hh --- dune-localfunctions-2.7.1/dune/python/localfunctions/localfiniteelement.hh 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/dune/python/localfunctions/localfiniteelement.hh 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,73 @@ +#ifndef DUNE_PYTHON_LOCALFUNCTIONS_LOCALFINITEELEMENT_HH +#define DUNE_PYTHON_LOCALFUNCTIONS_LOCALFINITEELEMENT_HH + +#include + +#include +#include +#include + +namespace Dune { +namespace Python { + +namespace detail { + +template +DUNE_EXPORT auto registerLocalBasis(pybind11::handle scope) +{ + static auto cls = pybind11::class_(scope, "LocalBasis"); + + cls.def("__len__", [](const LocalBasis& basis) { return basis.size(); }); + cls.def_property_readonly("order", [](const LocalBasis& basis) { return basis.order(); }); + cls.def("evaluateFunction", + [](const LocalBasis& basis, const typename LocalBasis::Traits::DomainType& in) { + std::vector out; + basis.evaluateFunction(in, out); + return out; + }); + cls.def("evaluateJacobian", + [](const LocalBasis& basis, const typename LocalBasis::Traits::DomainType& in) { + std::vector out; + basis.evaluateJacobian(in, out); + return out; + }); + return cls; +} + +DUNE_EXPORT auto registerLocalKey(pybind11::handle scope) +{ + static auto cls = pybind11::class_(scope, "LocalKey"); + + cls.def_property_readonly("subEntity", &LocalKey::subEntity); + cls.def_property_readonly("codim", &LocalKey::codim); + cls.def_property("index", + [](const LocalKey& key) { return key.index(); }, + [](LocalKey& key, unsigned int index) { key.index(index); }); + cls.def("__lt__", &LocalKey::operator<); + + return cls; +} + +} /* namespace detail */ + +template +DUNE_EXPORT auto registerLocalFiniteElement(pybind11::handle scope, const char* name = "LocalFiniteElement") +{ + static auto cls = pybind11::class_(scope, name); + + detail::registerLocalBasis(cls); + + cls.def_property_readonly("localBasis", &LocalFiniteElement::localBasis, pybind11::return_value_policy::reference_internal); + // cls.def_property_readonly("localCoefficients", &LocalFiniteElement::localCoefficients); + // cls.def_property_readonly("localInterpolation", &LocalFiniteElement::localInterpolation); + cls.def("__len__", &LocalFiniteElement::size); + cls.def_property_readonly("type", &LocalFiniteElement::type); + + return cls; +} + + +} /* namespace Python */ +} /* namespace Dune */ + +#endif diff -Nru dune-localfunctions-2.7.1/dune-localfunctions.pc.in dune-localfunctions-2.8.0/dune-localfunctions.pc.in --- dune-localfunctions-2.7.1/dune-localfunctions.pc.in 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune-localfunctions.pc.in 2021-08-31 14:03:38.000000000 +0000 @@ -8,8 +8,8 @@ Name: @PACKAGE_NAME@ Version: @VERSION@ -Description: Dune (Distributed and Unified Numerics Environment) local functions module -URL: http://dune-project.org/ +Description: @DESCRIPTION@ +URL: @URL@ Requires: ${DEPENDENCIES} Libs: -L${libdir} Cflags: -I${includedir} diff -Nru dune-localfunctions-2.7.1/dune.module dune-localfunctions-2.8.0/dune.module --- dune-localfunctions-2.7.1/dune.module 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/dune.module 2021-08-31 14:03:38.000000000 +0000 @@ -1,5 +1,9 @@ Module: dune-localfunctions -Version: 2.7.1 +Version: 2.8.0 +Author: The Dune Core developers Maintainer: dune-devel@lists.dune-project.org -Depends: dune-geometry (>= 2.7) +Description: Provides interface and implementation for shape functions defined on the DUNE reference elements +URL: https://gitlab.dune-project.org/core/dune-localfunctions +Python-Requires: +Depends: dune-geometry (>= 2.8) Whitespace-Hook: Yes diff -Nru dune-localfunctions-2.7.1/.gitignore dune-localfunctions-2.8.0/.gitignore --- dune-localfunctions-2.7.1/.gitignore 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/.gitignore 2021-08-31 14:03:38.000000000 +0000 @@ -1,2 +1,8 @@ *~ /build*/ +# ignore files generated during python setup.py sdist +MANIFEST +_skbuild/ +dist +# ignore macOS filesystem +.DS_Store diff -Nru dune-localfunctions-2.7.1/.gitlab-ci.yml dune-localfunctions-2.8.0/.gitlab-ci.yml --- dune-localfunctions-2.7.1/.gitlab-ci.yml 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/.gitlab-ci.yml 2021-08-31 14:03:38.000000000 +0000 @@ -2,12 +2,31 @@ include: - project: 'core/ci-config' ref: master - file: 'config/common/releases/2.7.yml' + file: 'config/common/master.yml' - project: 'core/ci-config' ref: master - file: 'jobs/common/releases/2.7.yml' + file: 'jobs/common/master.yml' before_script: - . /duneci/bin/duneci-init-job - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git - duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git + +debian-11-gcc-9-17-python: + image: registry.dune-project.org/docker/ci/debian:11 + script: duneci-standard-test + variables: + DUNECI_TOOLCHAIN: gcc-9-17 + # so we need some variables to build the dune-py module during execution of the first python test: + # we need to find the dune mdoule + DUNE_CONTROL_PATH: /duneci/modules:$CI_PROJECT_DIR + # the position for the dune-py module + DUNE_PY_DIR: /duneci/modules/dune-py + # during dune-py build this variable is used - do know a way to access + # the CMAKE_FLAGS used to build the modules... + DUNE_CMAKE_FLAGS: "-DCMAKE_CXX_COMPILER=g++-9 -DCMAKE_C_COMPILER=gcc-9 -DCMAKE_CXX_FLAGS='-std=c++17 -O2 -g -Wall -fdiagnostics-color=always' -DDUNE_ENABLE_PYTHONBINDINGS=ON -DDUNE_MAX_TEST_CORES=4 -DBUILD_SHARED_LIBS=TRUE -DDUNE_PYTHON_INSTALL_LOCATION=none -DCMAKE_POSITION_INDEPENDENT_CODE=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Vc=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE" + # cmake flags we use for all dune moudle - important is that build shared libs is set (need some better way of doing this) + DUNECI_CMAKE_FLAGS: $DUNE_CMAKE_FLAGS + # finally set the python path to all modules + PYTHONPATH: /duneci/modules/dune-common/build-cmake/python:/duneci/modules/dune-geometry/build-cmake/python:$CI_PROJECT_DIR/build-cmake/python + tags: [duneci] diff -Nru dune-localfunctions-2.7.1/LICENSE.md dune-localfunctions-2.8.0/LICENSE.md --- dune-localfunctions-2.7.1/LICENSE.md 2020-11-26 23:45:36.000000000 +0000 +++ dune-localfunctions-2.8.0/LICENSE.md 2021-08-31 14:03:38.000000000 +0000 @@ -7,33 +7,38 @@ 2016 Lukas Böger 2014 Andreas Buhr 2014--2017 Ansgar Burchardt -2009--2019 Andreas Dedner +2009--2021 Andreas Dedner +2020 Nils-Arne Dreier 2010 Martin Drohmann -2008--2018 Christian Engwer -2008--2017 Jorrit Fahlke +2008--2020 Christian Engwer +2008--2020 Jorrit Fahlke 2011--2013 Bernd Flemisch 2015 Elisa Friebel 2012 Christoph Gersbacher 2019 Janick Gerstenberger -2009--2019 Carsten Gräser +2009--2021 Carsten Gräser 2015--2018 Felix Gruber -2011--2019 Christoph Grüninger +2011--2021 Christoph Grüninger +2020 René Heß 2017 Patrick Jaap 2013 Guillaume Jouvet -2015--2019 Dominic Kempf +2015--2020 Dominic Kempf 2015 Angela Klewinghaus +2020 Robert Klöfkorn 2014 Arne Morten Kvarving 2015 Jizhou Li 2013--2017 Tobias Malkmus 2009 Sven Marnach 2013--2018 Steffen Müthing +2020 Lisa Julia Nebel 2012 Rebecca Neumann 2008--2018 Martin Nolte 2011--2016 Elias Pipping -2016--2019 Simon Praetorius +2016--2021 Simon Praetorius 2012--2013 Human Rezaijafari 2011--2012 Uli Sack -2008--2019 Oliver Sander +2008--2021 Oliver Sander +2020--2021 Henrik Stolzmann 2011--2012 Matthias Wohlmuth 2010--2015 Jonathan Youett diff -Nru dune-localfunctions-2.7.1/pyproject.toml dune-localfunctions-2.8.0/pyproject.toml --- dune-localfunctions-2.7.1/pyproject.toml 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/pyproject.toml 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,3 @@ +[build-system] +requires = ['setuptools', 'wheel', 'scikit-build', 'cmake', 'ninja', 'requests', 'dune-geometry>=2.8.0.dev0'] +build-backend = 'setuptools.build_meta' diff -Nru dune-localfunctions-2.7.1/python/CMakeLists.txt dune-localfunctions-2.8.0/python/CMakeLists.txt --- dune-localfunctions-2.7.1/python/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,3 @@ +add_subdirectory(dune) + +configure_file(setup.py.in setup.py) diff -Nru dune-localfunctions-2.7.1/python/dune/CMakeLists.txt dune-localfunctions-2.8.0/python/dune/CMakeLists.txt --- dune-localfunctions-2.7.1/python/dune/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/dune/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,5 @@ +add_subdirectory(localfunctions) + +add_python_targets(dune + __init__ +) diff -Nru dune-localfunctions-2.7.1/python/dune/__init__.py dune-localfunctions-2.8.0/python/dune/__init__.py --- dune-localfunctions-2.7.1/python/dune/__init__.py 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/dune/__init__.py 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1 @@ +__import__('pkg_resources').declare_namespace(__name__) diff -Nru dune-localfunctions-2.7.1/python/dune/localfunctions/CMakeLists.txt dune-localfunctions-2.8.0/python/dune/localfunctions/CMakeLists.txt --- dune-localfunctions-2.7.1/python/dune/localfunctions/CMakeLists.txt 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/dune/localfunctions/CMakeLists.txt 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,5 @@ +add_python_targets(localfunctions + __init__ +) + +dune_add_pybind11_module(NAME _localfunctions) diff -Nru dune-localfunctions-2.7.1/python/dune/localfunctions/__init__.py dune-localfunctions-2.8.0/python/dune/localfunctions/__init__.py --- dune-localfunctions-2.7.1/python/dune/localfunctions/__init__.py 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/dune/localfunctions/__init__.py 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,3 @@ +from __future__ import absolute_import, unicode_literals + +from ._localfunctions import * diff -Nru dune-localfunctions-2.7.1/python/dune/localfunctions/_localfunctions.cc dune-localfunctions-2.8.0/python/dune/localfunctions/_localfunctions.cc --- dune-localfunctions-2.7.1/python/dune/localfunctions/_localfunctions.cc 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/dune/localfunctions/_localfunctions.cc 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,7 @@ +#include +#include + +PYBIND11_MODULE(_localfunctions, module) +{ + Dune::Python::detail::registerLocalKey(module); +} diff -Nru dune-localfunctions-2.7.1/python/setup.py.in dune-localfunctions-2.8.0/python/setup.py.in --- dune-localfunctions-2.7.1/python/setup.py.in 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/python/setup.py.in 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,15 @@ +from setuptools import setup, find_packages + +pkg = [m for m in "${ProjectPythonRequires}".split(' ') if "dune" not in m] + +setup(name="${ProjectName}", + namespace_packages=['dune'], + description="${ProjectDescription}", + version="${ProjectVersionString}", + author="${ProjectAuthor}", + author_email="${ProjectMaintainerEmail}", + packages = find_packages(), + zip_safe = 0, + package_data = {'': ['*.so']}, + install_requires = pkg + ) diff -Nru dune-localfunctions-2.7.1/setup.py dune-localfunctions-2.8.0/setup.py --- dune-localfunctions-2.7.1/setup.py 1970-01-01 00:00:00.000000000 +0000 +++ dune-localfunctions-2.8.0/setup.py 2021-08-31 14:03:38.000000000 +0000 @@ -0,0 +1,6 @@ +try: + from dune.packagemetadata import metaData +except ImportError: + from packagemetadata import metaData +from skbuild import setup +setup(**metaData()[1])