From 42d61610ff22ba113928bcca106d41c9e2f42e19 Mon Sep 17 00:00:00 2001 From: Carlos Brito Date: Mon, 8 Jul 2024 00:29:05 +0200 Subject: [PATCH 1/4] ComplexFunction --- CMakeLists.txt | 1 - examples/PDEs/Helmholtz.cpp | 7 +- src/Rodin/Variational.h | 1 - src/Rodin/Variational/BoundaryNormal.h | 20 +- src/Rodin/Variational/ComplexFunction.h | 308 +++++++++++++++ src/Rodin/Variational/Div.h | 10 +- src/Rodin/Variational/Division.h | 2 +- src/Rodin/Variational/Dot.h | 6 +- src/Rodin/Variational/FaceNormal.h | 18 +- src/Rodin/Variational/ForwardDecls.h | 356 +++++++++-------- src/Rodin/Variational/Function.h | 138 +++---- src/Rodin/Variational/Grad.h | 16 +- src/Rodin/Variational/GridFunction.h | 454 +++++++++++++--------- src/Rodin/Variational/IdentityMatrix.h | 6 +- src/Rodin/Variational/Jacobian.h | 17 +- src/Rodin/Variational/LazyEvaluator.h | 8 + src/Rodin/Variational/MatrixFunction.h | 53 ++- src/Rodin/Variational/Mult.h | 4 +- src/Rodin/Variational/Normal.h | 63 --- src/Rodin/Variational/P0/GridFunction.h | 9 +- src/Rodin/Variational/P1/Div.h | 45 +-- src/Rodin/Variational/P1/Grad.h | 16 +- src/Rodin/Variational/P1/GridFunction.h | 9 +- src/Rodin/Variational/P1/Jacobian.h | 9 +- src/Rodin/Variational/P1/P1Element.cpp | 74 ++-- src/Rodin/Variational/P1/QuadratureRule.h | 2 +- src/Rodin/Variational/RealFunction.h | 29 +- src/Rodin/Variational/ScalarFunction.h | 282 +++++++++++++- src/Rodin/Variational/VectorFunction.h | 124 +++--- src/Rodin/Variational/Zero.h | 33 +- 30 files changed, 1434 insertions(+), 686 deletions(-) create mode 100644 src/Rodin/Variational/ComplexFunction.h delete mode 100644 src/Rodin/Variational/Normal.h diff --git a/CMakeLists.txt b/CMakeLists.txt index d8e5d6166..fc106b4b3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -79,7 +79,6 @@ rodin_add_compile_options(LANG CXX C OPTIONS -Wno-error=old-style-cast -Wno-erro add_compile_definitions(_LIBCPP_ENABLE_CXX17_REMOVED_UNARY_BINARY_FUNCTION) add_compile_definitions(__MAC_OS_X_VERSION_MIN_REQUIRED=__MAC_10_8) - if (RODIN_LTO) rodin_add_compile_options(LANG CXX C OPTIONS -flto) endif() diff --git a/examples/PDEs/Helmholtz.cpp b/examples/PDEs/Helmholtz.cpp index b2d5ef43b..e6d94a64b 100644 --- a/examples/PDEs/Helmholtz.cpp +++ b/examples/PDEs/Helmholtz.cpp @@ -21,12 +21,11 @@ int main(int, char**) // Functions P1 vh(mesh); + GridFunction gf(vh); Math::Vector m(2); - m << Complex(3, 1), Complex(2, 9); - - std::cout << m << std::endl; - std::cout << m.dot(m) << std::endl; + // m << Complex(3, 1), Complex(2, 9); + m.setOnes(); // TrialFunction uRe(vh), uIm(vh); // TestFunction vRe(vh), vIm(vh); diff --git a/src/Rodin/Variational.h b/src/Rodin/Variational.h index 460f15827..abe0336bf 100644 --- a/src/Rodin/Variational.h +++ b/src/Rodin/Variational.h @@ -48,7 +48,6 @@ #include "Variational/Div.h" #include "Variational/Grad.h" -#include "Variational/Normal.h" #include "Variational/FaceNormal.h" #include "Variational/BoundaryNormal.h" #include "Variational/Jacobian.h" diff --git a/src/Rodin/Variational/BoundaryNormal.h b/src/Rodin/Variational/BoundaryNormal.h index 62a96cb50..84ca9dd35 100644 --- a/src/Rodin/Variational/BoundaryNormal.h +++ b/src/Rodin/Variational/BoundaryNormal.h @@ -14,10 +14,14 @@ namespace Rodin::Variational /** * @brief Outward unit normal. */ - class BoundaryNormal final : public VectorFunctionBase + class BoundaryNormal final : public VectorFunctionBase { public: - using Parent = VectorFunctionBase; + using ScalarType = Real; + + using VectorType = Math::SpatialVector; + + using Parent = VectorFunctionBase; /** * @brief Constructs the outward unit normal. @@ -65,14 +69,14 @@ namespace Rodin::Variational } inline - Math::SpatialVector getValue(const Geometry::Point& p) const + VectorType getValue(const Geometry::Point& p) const { - Math::SpatialVector res; + VectorType res; getValue(res, p); return res; } - void interpolate(Math::SpatialVector& res, const Geometry::Point& p) const + void interpolate(VectorType& res, const Geometry::Point& p) const { const auto& polytope = p.getPolytope(); const auto& d = polytope.getDimension(); @@ -90,7 +94,7 @@ namespace Rodin::Variational if (inc.size() == 1) { const auto& tracePolytope = mesh.getPolytope(meshDim - 1, *inc.begin()); - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const VectorType rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(res, np); return; @@ -115,7 +119,7 @@ namespace Rodin::Variational const auto& tracePolytope = mesh.getPolytope(meshDim - 1, idx); if (traceDomain.count(tracePolytope->getAttribute())) { - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const VectorType rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(res, np); return; @@ -167,7 +171,7 @@ namespace Rodin::Variational } } - void getValue(Math::SpatialVector& out, const Geometry::Point& p) const + void getValue(VectorType& out, const Geometry::Point& p) const { out.setConstant(NAN); const auto& polytope = p.getPolytope(); diff --git a/src/Rodin/Variational/ComplexFunction.h b/src/Rodin/Variational/ComplexFunction.h new file mode 100644 index 000000000..e6e036349 --- /dev/null +++ b/src/Rodin/Variational/ComplexFunction.h @@ -0,0 +1,308 @@ +/* + * Copyright Carlos BRITO PACHECO 2021 - 2022. + * Distributed under the Boost Software License, Version 1.0. + * (See accompanying file LICENSE or copy at + * https://www.boost.org/LICENSE_1_0.txt) + */ +#ifndef RODIN_VARIATIONAL_COMPLEXFUNCTION_H +#define RODIN_VARIATIONAL_COMPLEXFUNCTION_H + +#include +#include +#include +#include +#include + +#include "Rodin/Geometry/Polytope.h" + +#include "ForwardDecls.h" + +#include "RangeShape.h" + +#include "ScalarFunction.h" + + +namespace Rodin::Variational +{ + /** + * @defgroup RealFunctionSpecializations RealFunction Template Specializations + * @brief Template specializations of the RealFunction class. + * @see RealFunction + */ + + template + class ComplexFunctionBase : public ScalarFunctionBase> + { + public: + using ScalarType = Complex; + + using Parent = ScalarFunctionBase>; + + using Parent::traceOf; + + ComplexFunctionBase() = default; + + ComplexFunctionBase(const ComplexFunctionBase& other) + : Parent(other) + {} + + ComplexFunctionBase(ComplexFunctionBase&& other) + : Parent(std::move(other)) + {} + + virtual ~ComplexFunctionBase() = default; + + inline + const Derived& getDerived() const + { + return static_cast(*this); + } + + inline + constexpr + auto getValue(const Geometry::Point& p) const + { + return static_cast(*this).getValue(p); + } + + virtual ComplexFunctionBase* copy() const noexcept override = 0; + }; + + /** + * @ingroup ComplexFunctionSpecializations + */ + template + class ComplexFunction> final + : public ComplexFunctionBase> + { + public: + using Parent = ComplexFunctionBase>; + + using NestedRangeType = typename FormLanguage::Traits>::RangeType; + + static_assert(std::is_same_v); + + ComplexFunction(const FunctionBase& nested) + : m_nested(nested.copy()) + {} + + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_nested(other.m_nested->copy()) + {} + + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_nested(std::move(other.m_nested)) + {} + + inline + constexpr + auto getValue(const Geometry::Point& v) const + { + return m_nested->getValue(v); + } + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute attrs) + { + m_nested->traceOf(attrs); + return *this; + } + + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + std::unique_ptr> m_nested; + }; + + /** + * @brief CTAD for ComplexFunction. + */ + template + ComplexFunction(const FunctionBase&) -> ComplexFunction>; + + /** + * @ingroup ComplexFunctionSpecializations + * @brief Represents a constant scalar function with type Complex. + */ + template <> + class ComplexFunction final + : public ComplexFunctionBase> + { + public: + using Parent = ComplexFunctionBase>; + + /** + * @brief Constructs a ComplexFunction from a constant scalar value. + * @param[in] x Constant scalar value + */ + ComplexFunction(Complex x) + : m_x(x) + {} + + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_x(other.m_x) + {} + + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_x(other.m_x) + {} + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + Complex getValue() const + { + return m_x; + } + + inline + constexpr + Complex getValue(const Geometry::Point&) const + { + return m_x; + } + + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + const Complex m_x; + }; + + /** + * @brief CTAD for ComplexFunction. + */ + ComplexFunction(Complex) -> ComplexFunction; + + template <> + class ComplexFunction final + : public ComplexFunctionBase> + { + public: + using Parent = ComplexFunctionBase>; + + /** + * @brief Constructs a ComplexFunction from an integer value. + * @param[in] x Constant integer value + */ + ComplexFunction(Integer x) + : m_x(x) + {} + + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_x(other.m_x) + {} + + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_x(other.m_x) + {} + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + Complex getValue() const + { + return m_x; + } + + inline + constexpr + Complex getValue(const Geometry::Point&) const + { + return static_cast(m_x); + } + + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + const Integer m_x; + }; + + ComplexFunction(Integer) -> ComplexFunction; + + /** + * @ingroup ComplexFunctionSpecializations + * @brief Represents a scalar function given by an arbitrary scalar function. + */ + template + class ComplexFunction final : public ComplexFunctionBase> + { + static_assert(std::is_invocable_r_v); + + public: + using Parent = ComplexFunctionBase>; + + ComplexFunction(F f) + : m_f(f) + {} + + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_f(other.m_f) + {} + + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_f(std::move(other.m_f)) + {} + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + Complex getValue(const Geometry::Point& v) const + { + return m_f(v); + } + + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + const F m_f; + }; + + /** + * @brief CTAD for ComplexFunction. + */ + template >> + ComplexFunction(F) -> ComplexFunction; +} + +#endif + diff --git a/src/Rodin/Variational/Div.h b/src/Rodin/Variational/Div.h index 99335665c..ced60d2f6 100644 --- a/src/Rodin/Variational/Div.h +++ b/src/Rodin/Variational/Div.h @@ -20,16 +20,16 @@ namespace Rodin::Variational /** * @brief Base class for Div classes. */ - template + template class DivBase; /** * @ingroup DivSpecializations * @brief Divergence of a P1 GridFunction */ - template - class DivBase> - : public RealFunctionBase>> + template + class DivBase, Derived> + : public RealFunctionBase, Derived>> { public: using FESType = FES; @@ -37,7 +37,7 @@ namespace Rodin::Variational using OperandType = GridFunction; /// Parent class - using Parent = RealFunctionBase>; + using Parent = RealFunctionBase>; /** * @brief Constructs the Div of a @f$ \mathbb{P}_1 @f$ function @f$ u diff --git a/src/Rodin/Variational/Division.h b/src/Rodin/Variational/Division.h index 43e188bcd..41dd6d73d 100644 --- a/src/Rodin/Variational/Division.h +++ b/src/Rodin/Variational/Division.h @@ -96,7 +96,7 @@ namespace Rodin::Variational constexpr void getValue(Math::Vector& res, const Geometry::Point& p) const { - getLHS().getValue(res, p); + getLHS().getDerived().getValue(res, p); res /= getRHS().getValue(p); } diff --git a/src/Rodin/Variational/Dot.h b/src/Rodin/Variational/Dot.h index d0eb0d6be..a1db812c6 100644 --- a/src/Rodin/Variational/Dot.h +++ b/src/Rodin/Variational/Dot.h @@ -182,14 +182,14 @@ namespace Rodin::Variational assert(getLHS().getRangeShape() == getRHS().getRangeShape()); const auto& lhs = this->object(getLHS().getValue(p)); const auto& rhs = this->object(getRHS().getValue(p)); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { return lhs * rhs; - } else if constexpr (std::is_same_v>) + } else if constexpr (std::is_same_v>) { return lhs.dot(rhs); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { return (lhs.array() * rhs.array()).rowwise().sum().colwise().sum().value(); } diff --git a/src/Rodin/Variational/FaceNormal.h b/src/Rodin/Variational/FaceNormal.h index bfdeb3d76..19ca0ef97 100644 --- a/src/Rodin/Variational/FaceNormal.h +++ b/src/Rodin/Variational/FaceNormal.h @@ -15,10 +15,14 @@ namespace Rodin::Variational /** * @brief Outward unit normal on a face. */ - class FaceNormal : public VectorFunctionBase + class FaceNormal : public VectorFunctionBase { public: - using Parent = VectorFunctionBase; + using ScalarType = Real; + + using VectorType = Math::SpatialVector; + + using Parent = VectorFunctionBase; /** * @brief Constructs the outward unit on a face. @@ -47,14 +51,14 @@ namespace Rodin::Variational } inline - Math::SpatialVector getValue(const Geometry::Point& p) const + VectorType getValue(const Geometry::Point& p) const { - Math::SpatialVector res; + VectorType res; getValue(res, p); return res; } - void getValue(Math::SpatialVector& res, const Geometry::Point& p) const + void getValue(VectorType& res, const Geometry::Point& p) const { const auto& polytope = p.getPolytope(); const auto& vs = p.getPolytope().getVertices(); @@ -74,9 +78,9 @@ namespace Rodin::Variational { const Index v1 = vs[0]; const Index v2 = vs[1]; - Eigen::Vector3 a = + Eigen::Vector3 a = mesh.getVertexCoordinates(v1) - mesh.getVertexCoordinates(v2); - Eigen::Vector3 n; + Eigen::Vector3 n; n << jacobian(1, 0), -jacobian(0, 0), jacobian(2, 0); n = n.cross(a); n.stableNormalize(); diff --git a/src/Rodin/Variational/ForwardDecls.h b/src/Rodin/Variational/ForwardDecls.h index df5f90207..8b0b7ca7e 100644 --- a/src/Rodin/Variational/ForwardDecls.h +++ b/src/Rodin/Variational/ForwardDecls.h @@ -19,152 +19,6 @@ namespace Rodin::Variational */ enum class RangeType; - /** - * @brief Base class for linear form objects. - * @tparam VectorType Type of vector which will be assembled - */ - template - class LinearFormBase; - - /** - * @brief Represents a linear form on some finite element space. - * @tparam FES Type of finite element space - * @tparam VectorType Type of vector which will be assembled - * - * Represents a linear form @f$ \ell : V_h \rightarrow \mathbb{R} @f$ on a given - * finite element space @f$ V_h @f$. - */ - template - class LinearForm; - - class Integrator; - - /** - * @brief Base class for linear form integrators. - * - * An instance of LinearFormIntegratorBase performs the assembly of the - * element vector for each finite element. - */ - template - class LinearFormIntegratorBase; - - /** - * @brief Abstract base class for objects of type BilinearForm. - * @tparam OperatorType Type of operator which will be assembled - * - * Represents a bilinear form: - * @f[ - * a : V_h \times U_h \rightarrow \mathbb{R} - * @f] - * on given finite element spaces @f$ V_h @f$ and @f$ U_h @f$. - */ - template - class BilinearFormBase; - - /** - * @brief Represents a bilinear form on a trial and test space - * @tparam TrialFES Type of trial finite element space - * @tparam TestFES Type of test finite element space - * - * An object of type BilinearForm represents a bilinear map: - * @f[ - * \begin{aligned} - * a : U \times V &\rightarrow \mathbb{R}\\ - * (u, v) &\mapsto a(u, v) - * \end{aligned} - * @f] - * where @f$ U @f$ and @f$ V @f$ are finite element spaces. - * - * @see BilinearFormSpecializations - */ - template - class BilinearForm; - - /** - * @brief Base class for bilinear form integrators. - */ - template - class BilinearFormIntegratorBase; - - template - class LocalBilinearFormIntegratorBase; - - template - class GlobalBilinearFormIntegratorBase; - - template - class FiniteElementBase; - - template - class FiniteElement; - - class FiniteElementCollectionBase; - - /** - * @brief Base class for finite element spaces. - */ - class FiniteElementSpaceBase; - - /** - * @brief Arbitrary order @f$ H^1(\mathcal{T}_h)^d \subset L^2 (\Omega) @f$ - * broken Sobolev space. - * @tparam Context Type of context for the finite element space - * - * Given some triangulation @f$ \mathcal{T}_h @f$ of @f$ \Omega @f$, - * instances of this class will represent the finite element space: - * @f[ - * H^1(\mathcal{T}_h)^d := \left\{ v \in L^2(\Omega)^d : \forall \tau \in - * \mathcal{T}_h, \ v |_\tau \in H^1 (\tau)^d \right\} - * @f] - */ - template - class L2; - - /** - * @brief Represents the lazy evaluation of a mesh function. - * - * The main objective of this class is to wrap the reference of a data-full - * object into a light object which permits evaluation. This way one is able - * to call the `copy()` method without actually copying the underlying data. - */ - template - class LazyEvaluator; - - /** - * @brief Base class for grid function objects. - */ - template - class GridFunctionBase; - - /** - * @brief Represents a grid function belonging to some finite element space. - * @tparam FES Type of finite element space - * - * A GridFunction object represents a function whose values are known at the - * "grid points". These grid points are the node coordinates of the global - * degrees of freedom in the finite element space. - * - * @section gridfunction-interpolation Interpolation - * - * In general, the GridFunction class represents the global interpolation - * operator on the @f$ \mathcal{T}_h @f$-based family of finite elements - * @f$ \{ K, P_K, \Sigma_K \}_{K \in \mathcal{T}_h } @f$. More precisely, it - * defines the global interpolation operator @f$ \mathcal{I}_h : - * D(\mathcal{I}_h) \rightarrow V_h @f$, where: - * @f[ - * D(\mathcal{I}_h) := \{ v \in V(\mathcal{T}_h) : \forall K \in - * \mathcal{T}_h, \: v|_K \in V(K) \} \: . - * @f] - * - * @note For an overview of all the possible specializations of the - * GridFunction class, please see @ref GridFunctionSpecializations. - * - * @see GridFunctionBase - * @see GridFunctionSpecializations - */ - template - class GridFunction; - /** * @brief Enumeration class to indicate whether a derived instance of * ShapeFunctionBase belongs to either a trial or test space. @@ -233,7 +87,43 @@ namespace Rodin::Variational class Function; /** - * @brief Base class for scalar-valued functions defined on a mesh. + * @brief Base class for objects representing boolean functions. + */ + template + class BooleanFunctionBase; + + /** + * @note For an overview of all the possible specializations of the + * BooleanFunction class, please see @ref BooleanFunctionSpecializations. + * + * @see BooleanFunctionSpecializations + */ + template + class BooleanFunction; + + /** + * @brief Base class for objects representing integer functions. + */ + template + class IntegerFunctionBase; + + /** + * @note For an overview of all the possible specializations of the + * IntegerFunction class, please see @ref IntegerFunctionSpecializations. + * + * @see IntegerFunctionSpecializations + */ + template + class IntegerFunction; + + template + class ScalarFunctionBase; + + template + class ScalarFunction; + + /** + * @brief Base class for real-valued functions defined on a mesh. */ template class RealFunctionBase; @@ -247,8 +137,20 @@ namespace Rodin::Variational template class RealFunction; - template - class Zero; + /** + * @brief Base class for scalar-valued functions defined on a mesh. + */ + template + class ComplexFunctionBase; + + /** + * @note For an overview of all the possible specializations of the + * ComplexFunction class, please see @ref ComplexFunctionSpecializations. + * + * @see ComplexFunctionSpecializations + */ + template + class ComplexFunction; /** * @brief Base class for vector-valued functions defined on a mesh. @@ -256,7 +158,7 @@ namespace Rodin::Variational * @note Vectors are zero indexed. This means that the 0-index corresponds * to the 1st entry of the vector. */ - template + template class VectorFunctionBase; /** @@ -271,7 +173,7 @@ namespace Rodin::Variational /** * @brief Base class for matrix-valued functions defined on a mesh. */ - template + template class MatrixFunctionBase; /** @@ -284,19 +186,151 @@ namespace Rodin::Variational class MatrixFunction; /** - * @brief Base class for objects representing boolean functions. + * @brief Base class for grid function objects. */ - template - class BooleanFunctionBase; + template + class GridFunctionBase; /** + * @brief Represents a grid function belonging to some finite element space. + * @tparam FES Type of finite element space + * + * A GridFunction object represents a function whose values are known at the + * "grid points". These grid points are the node coordinates of the global + * degrees of freedom in the finite element space. + * + * @section gridfunction-interpolation Interpolation + * + * In general, the GridFunction class represents the global interpolation + * operator on the @f$ \mathcal{T}_h @f$-based family of finite elements + * @f$ \{ K, P_K, \Sigma_K \}_{K \in \mathcal{T}_h } @f$. More precisely, it + * defines the global interpolation operator @f$ \mathcal{I}_h : + * D(\mathcal{I}_h) \rightarrow V_h @f$, where: + * @f[ + * D(\mathcal{I}_h) := \{ v \in V(\mathcal{T}_h) : \forall K \in + * \mathcal{T}_h, \: v|_K \in V(K) \} \: . + * @f] + * * @note For an overview of all the possible specializations of the - * BooleanFunction class, please see @ref BooleanFunctionSpecializations. + * GridFunction class, please see @ref GridFunctionSpecializations. * - * @see BooleanFunctionSpecializations + * @see GridFunctionBase + * @see GridFunctionSpecializations */ - template - class BooleanFunction; + template + class GridFunction; + + /** + * @brief Base class for linear form objects. + * @tparam VectorType Type of vector which will be assembled + */ + template + class LinearFormBase; + + /** + * @brief Represents a linear form on some finite element space. + * @tparam FES Type of finite element space + * @tparam VectorType Type of vector which will be assembled + * + * Represents a linear form @f$ \ell : V_h \rightarrow \mathbb{R} @f$ on a given + * finite element space @f$ V_h @f$. + */ + template + class LinearForm; + + class Integrator; + + /** + * @brief Base class for linear form integrators. + * + * An instance of LinearFormIntegratorBase performs the assembly of the + * element vector for each finite element. + */ + template + class LinearFormIntegratorBase; + + /** + * @brief Abstract base class for objects of type BilinearForm. + * @tparam OperatorType Type of operator which will be assembled + * + * Represents a bilinear form: + * @f[ + * a : V_h \times U_h \rightarrow \mathbb{R} + * @f] + * on given finite element spaces @f$ V_h @f$ and @f$ U_h @f$. + */ + template + class BilinearFormBase; + + /** + * @brief Represents a bilinear form on a trial and test space + * @tparam TrialFES Type of trial finite element space + * @tparam TestFES Type of test finite element space + * + * An object of type BilinearForm represents a bilinear map: + * @f[ + * \begin{aligned} + * a : U \times V &\rightarrow \mathbb{R}\\ + * (u, v) &\mapsto a(u, v) + * \end{aligned} + * @f] + * where @f$ U @f$ and @f$ V @f$ are finite element spaces. + * + * @see BilinearFormSpecializations + */ + template + class BilinearForm; + + /** + * @brief Base class for bilinear form integrators. + */ + template + class BilinearFormIntegratorBase; + + template + class LocalBilinearFormIntegratorBase; + + template + class GlobalBilinearFormIntegratorBase; + + template + class FiniteElementBase; + + template + class FiniteElement; + + /** + * @brief Base class for finite element spaces. + */ + class FiniteElementSpaceBase; + + /** + * @brief Arbitrary order @f$ H^1(\mathcal{T}_h)^d \subset L^2 (\Omega) @f$ + * broken Sobolev space. + * @tparam Context Type of context for the finite element space + * + * Given some triangulation @f$ \mathcal{T}_h @f$ of @f$ \Omega @f$, + * instances of this class will represent the finite element space: + * @f[ + * H^1(\mathcal{T}_h)^d := \left\{ v \in L^2(\Omega)^d : \forall \tau \in + * \mathcal{T}_h, \ v |_\tau \in H^1 (\tau)^d \right\} + * @f] + */ + template + class L2; + + /** + * @brief Represents the lazy evaluation of a mesh function. + * + * The main objective of this class is to wrap the reference of a data-full + * object into a light object which permits evaluation. This way one is able + * to call the `copy()` method without actually copying the underlying data. + */ + template + class LazyEvaluator; + + template + class Zero; template class Jump; @@ -716,6 +750,9 @@ namespace Rodin::Variational template class TraceOperator; + template + class Potential; + template class Max; @@ -1036,9 +1073,6 @@ namespace Rodin::Variational template class DenseProblem; - - template - class Potential; } #endif diff --git a/src/Rodin/Variational/Function.h b/src/Rodin/Variational/Function.h index 40fef34e3..f61c50957 100644 --- a/src/Rodin/Variational/Function.h +++ b/src/Rodin/Variational/Function.h @@ -158,37 +158,37 @@ namespace Rodin::Variational return static_cast(*this).getRangeShape(); } - inline - constexpr - RangeType getRangeType() const - { - using R = typename FormLanguage::Traits>::RangeType; - if constexpr (std::is_same_v) - { - return RangeType::Boolean; - } - else if constexpr (std::is_same_v) - { - return RangeType::Integer; - } - else if constexpr (std::is_same_v) - { - return RangeType::Real; - } - else if constexpr (std::is_same_v>) - { - return RangeType::Vector; - } - else if constexpr (std::is_same_v>) - { - return RangeType::Matrix; - } - else - { - assert(false); - static_assert(Utility::DependentFalse::Value); - } - } + // inline + // constexpr + // RangeType getRangeType() const + // { + // using R = typename FormLanguage::Traits>::RangeType; + // if constexpr (std::is_same_v) + // { + // return RangeType::Boolean; + // } + // else if constexpr (std::is_same_v) + // { + // return RangeType::Integer; + // } + // else if constexpr (std::is_same_v) + // { + // return RangeType::Real; + // } + // else if constexpr (Utility::IsSpecialization) + // { + // return RangeType::Vector; + // } + // else if constexpr (std::is_same_v>) + // { + // return RangeType::Matrix; + // } + // else + // { + // assert(false); + // static_assert(Utility::DependentFalse::Value); + // } + // } /** * @brief Evaluates the function on a Point belonging to the mesh. @@ -201,33 +201,33 @@ namespace Rodin::Variational return static_cast(*this).getValue(p); } - inline - constexpr - void getValue(Math::Vector& res, const Geometry::Point& p) const - { - if constexpr (Internal::HasGetValueMethod&, const Geometry::Point&>::Value) - { - return static_cast(*this).getValue(res, p); - } - else - { - res = getValue(p); - } - } - - inline - constexpr - void getValue(Math::Matrix& res, const Geometry::Point& p) const - { - if constexpr (Internal::HasGetValueMethod&, const Geometry::Point&>::Value) - { - return static_cast(*this).getValue(res, p); - } - else - { - res = getValue(p); - } - } + // inline + // constexpr + // void getValue(Math::Vector& res, const Geometry::Point& p) const + // { + // if constexpr (Internal::HasGetValueMethod&, const Geometry::Point&>::Value) + // { + // return static_cast(*this).getValue(res, p); + // } + // else + // { + // res = getValue(p); + // } + // } + + // inline + // constexpr + // void getValue(Math::Matrix& res, const Geometry::Point& p) const + // { + // if constexpr (Internal::HasGetValueMethod&, const Geometry::Point&>::Value) + // { + // return static_cast(*this).getValue(res, p); + // } + // else + // { + // res = getValue(p); + // } + // } inline auto coeff(size_t i, size_t j) const @@ -241,17 +241,17 @@ namespace Rodin::Variational return Component(*this, i); } - inline - void operator()(Math::Vector& res, const Geometry::Point& p) const - { - return getValue(res, p); - } - - inline - void operator()(Math::Matrix& res, const Geometry::Point& p) const - { - return getValue(res, p); - } + // inline + // void operator()(Math::Vector& res, const Geometry::Point& p) const + // { + // return getValue(res, p); + // } + + // inline + // void operator()(Math::Matrix& res, const Geometry::Point& p) const + // { + // return getValue(res, p); + // } /** * @brief Evaluates the function on a Point belonging to the mesh. diff --git a/src/Rodin/Variational/Grad.h b/src/Rodin/Variational/Grad.h index 98d197306..75b3c562d 100644 --- a/src/Rodin/Variational/Grad.h +++ b/src/Rodin/Variational/Grad.h @@ -28,24 +28,28 @@ namespace Rodin::Variational /** * @brief Base class for Grad classes. */ - template + template class GradBase; /** * @ingroup GradSpecializations * @brief Gradient of a P1 GridFunction */ - template - class GradBase> - : public VectorFunctionBase>> + template + class GradBase, Derived> + : public VectorFunctionBase< + typename FormLanguage::Traits::ScalarType, GradBase, Derived>> { public: using FESType = FES; + using ScalarType = typename FormLanguage::Traits::ScalarType; + + using VectorType = Math::SpatialVector; + using OperandType = GridFunction; - /// Parent class - using Parent = VectorFunctionBase>; + using Parent = VectorFunctionBase>; /** * @brief Constructs the gradient of a @f$ \mathbb{P}_1 @f$ function @f$ diff --git a/src/Rodin/Variational/GridFunction.h b/src/Rodin/Variational/GridFunction.h index 708911a35..1701acfc6 100644 --- a/src/Rodin/Variational/GridFunction.h +++ b/src/Rodin/Variational/GridFunction.h @@ -32,15 +32,15 @@ #include "Component.h" #include "Restriction.h" #include "LazyEvaluator.h" -#include "RealFunction.h" +#include "ScalarFunction.h" #include "VectorFunction.h" #include "MatrixFunction.h" #include "FiniteElementSpace.h" namespace Rodin::FormLanguage { - template - struct Traits> + template + struct Traits> { using FESType = FES; using MeshType = typename Traits::MeshType; @@ -70,18 +70,18 @@ namespace Rodin::Variational * @section gridfunction-data-layout Data layout * * The data of the GridFunctionBase object can be accessed via a call to @ref - * getData(). The i-th column of the returned Math::Matrix object corresponds + * getData(). The i-th column of the returned Math::Matrix object corresponds * to the value of the grid function at the global i-th degree of freedom in * the finite element space. Furthermore, the following conditions are * satisfied: * ``` - * const Math::Matrix& data = gf.getData(); + * const auto& data = gf.getData(); * assert(data.rows() == gf.getFiniteElementSpace().getVectorDimension()); * assert(data.cols() == gf.getFiniteElementSpace().getSize()); * ``` */ - template ::FESType> - class GridFunctionBase : public LazyEvaluator> + template + class GridFunctionBase : public LazyEvaluator> { public: using FESType = FES; @@ -100,12 +100,15 @@ namespace Rodin::Variational /// Type of finite element using ElementType = typename FormLanguage::Traits::ElementType; + using DataType = Math::Matrix; + + using WeightVectorType = Math::Vector; + /// Parent class - using Parent = LazyEvaluator>; + using Parent = LazyEvaluator>; - static_assert( - std::is_same_v || - std::is_same_v>); + static_assert(std::is_same_v || + std::is_same_v>); GridFunctionBase(const FES& fes) : Parent(std::cref(*this)), @@ -152,19 +155,15 @@ namespace Rodin::Variational */ inline constexpr - Real max() const + ScalarType max() const { - static_assert(std::is_same_v, - "GridFunction must be scalar valued."); return m_data.maxCoeff(); } inline constexpr - Real max(Index& idx) const + ScalarType max(Index& idx) const { - static_assert(std::is_same_v, - "GridFunction must be scalar valued."); return m_data.maxCoeff(&idx); } @@ -181,19 +180,15 @@ namespace Rodin::Variational */ inline constexpr - Real min() const + ScalarType min() const { - static_assert(std::is_same_v, - "GridFunction must be scalar valued."); return m_data.minCoeff(); } inline constexpr - Real min(Index& idx) const + ScalarType min(Index& idx) const { - static_assert(std::is_same_v, - "GridFunction must be scalar valued."); return m_data.minCoeff(&idx); } @@ -201,8 +196,6 @@ namespace Rodin::Variational constexpr Index argmax() const { - static_assert(std::is_same_v, - "GridFunction must be scalar valued."); Index idx = 0; m_data.maxCoeff(&idx); return idx; @@ -212,8 +205,6 @@ namespace Rodin::Variational constexpr Index argmin() const { - static_assert(std::is_same_v, - "GridFunction must be scalar valued."); Index idx = 0; m_data.minCoeff(&idx); return idx; @@ -222,7 +213,7 @@ namespace Rodin::Variational inline Derived& normalize() { - static_assert(std::is_same_v>, + static_assert(std::is_same_v>, "GridFunction must be vector valued."); for (size_t i = 0; i < getSize(); i++) getData().col(i).normalize(); @@ -232,7 +223,7 @@ namespace Rodin::Variational inline Derived& stableNormalize() { - static_assert(std::is_same_v>, + static_assert(std::is_same_v>, "GridFunction must be vector valued."); for (size_t i = 0; i < getSize(); i++) getData().col(i).stableNormalize(); @@ -250,7 +241,7 @@ namespace Rodin::Variational constexpr auto x() const { - static_assert(std::is_same_v>); + static_assert(std::is_same_v>); assert(getFiniteElementSpace().getVectorDimension() >= 1); return Component(*this, 0); } @@ -259,7 +250,7 @@ namespace Rodin::Variational constexpr auto y() const { - static_assert(std::is_same_v>); + static_assert(std::is_same_v>); assert(getFiniteElementSpace().getVectorDimension() >= 2); return Component(*this, 1); } @@ -268,7 +259,7 @@ namespace Rodin::Variational constexpr auto z() const { - static_assert(std::is_same_v>); + static_assert(std::is_same_v>); assert(getFiniteElementSpace().getVectorDimension() >= 3); return Component(*this, 2); } @@ -294,41 +285,35 @@ namespace Rodin::Variational * @brief Bulk assigns the value to the whole data array. */ inline - Derived& operator=(Real v) + Derived& operator=(const RangeType& v) { - static_assert(std::is_same_v); - m_data.setConstant(v); + if constexpr (std::is_same_v) + { + m_data.setConstant(v); + } + else if constexpr (std::is_same_v>) + { + assert(m_data.cols() >= 0); + for (size_t i = 0; i < static_cast(m_data.cols()); i++) + m_data.col(i) = v; + } + else + { + assert(false); + } return static_cast(*this); } inline - Derived& operator=(const Math::Vector& v) + Derived& operator=(std::function fn) { - static_assert(std::is_same_v>); - Math::Matrix& data = m_data; - assert(data.cols() >= 0); - for (size_t i = 0; i < static_cast(data.cols()); i++) - data.col(i) = v; - return static_cast(*this); + return project(fn); } inline - Derived& operator=(std::function fn) + Derived& operator=(std::function fn) { - if constexpr (std::is_same_v) - { - assert(getFiniteElementSpace().getVectorDimension() == 1); - return project(RealFunction(fn)); - } - else if constexpr (std::is_same_v>) - { - return project(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); - } - else - { - assert(false); - return static_cast(*this); - } + return project(fn); } /** @@ -345,9 +330,9 @@ namespace Rodin::Variational * @brief Addition of a scalar value. */ inline - Derived& operator+=(Real rhs) + Derived& operator+=(const ScalarType& rhs) { - static_assert(std::is_same_v); + static_assert(std::is_same_v); m_data = m_data.array() + rhs; return static_cast(*this); } @@ -356,9 +341,9 @@ namespace Rodin::Variational * @brief Substraction of a scalar value. */ inline - Derived& operator-=(Real rhs) + Derived& operator-=(const ScalarType& rhs) { - static_assert(std::is_same_v); + static_assert(std::is_same_v); m_data = m_data.array() - rhs; return static_cast(*this); } @@ -367,7 +352,7 @@ namespace Rodin::Variational * @brief Multiplication by a scalar value. */ inline - Derived& operator*=(Real rhs) + Derived& operator*=(const ScalarType& rhs) { m_data = m_data.array() * rhs; return static_cast(*this); @@ -377,7 +362,7 @@ namespace Rodin::Variational * @brief Division by a scalar value. */ inline - Derived& operator/=(Real rhs) + Derived& operator/=(const ScalarType& rhs) { m_data = m_data.array() / rhs; return static_cast(*this); @@ -388,7 +373,7 @@ namespace Rodin::Variational { if (this == &rhs) { - operator*=(Real(2)); + operator*=(2); } else { @@ -403,7 +388,7 @@ namespace Rodin::Variational { if (this == &rhs) { - operator=(Real(0)); + m_data.setZero(); } else { @@ -433,7 +418,7 @@ namespace Rodin::Variational { if (this == &rhs) { - operator=(Real(1)); + m_data.setOnes(); } else { @@ -446,25 +431,63 @@ namespace Rodin::Variational /** * @brief Projects a scalar valued function on the region of the mesh * with the given attribute. - * @param[in] fn Real valued function + * @param[in] fn Scalar valued function * @param[in] attr Attribute */ inline - auto& project(std::function fn, Geometry::Attribute attr) + auto& project( + std::function fn, Geometry::Attribute attr) { - return project(fn, FlatSet{attr}); + return project(fn, FlatSet{ attr }); + } + + inline + auto& project( + std::function fn, Geometry::Attribute attr) + { + return project(fn, FlatSet{ attr }); + } + + inline + auto& project( + std::function fn, + const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return project(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return project(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } } - /** - * @brief Projects a scalar valued function on the region of the mesh - * with the given attributes. - * @param[in] fn Real valued function - * @param[in] attrs Set of attributes - */ inline - auto& project(std::function fn, const FlatSet& attrs = {}) + auto& project( + std::function fn, + const FlatSet& attrs = {}) { - return project(RealFunction(fn), attrs); + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return project(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return project(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } } template @@ -508,7 +531,7 @@ namespace Rodin::Variational const auto& mesh = fes.getMesh(); const size_t d = mesh.getDimension(); std::vector ns(fes.getSize(), 0); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { #ifdef RODIN_MULTITHREADED auto& threadPool = Threads::getGlobalThreadPool(); @@ -518,7 +541,7 @@ namespace Rodin::Variational const size_t capacity = fes.getSize() / threadPool.getThreadCount(); std::vector is; is.reserve(capacity); - std::vector vs; + std::vector vs; vs.reserve(capacity); std::unique_ptr> fnt(fn.copy()); for (Index i = start; i < end; ++i) @@ -574,9 +597,9 @@ namespace Rodin::Variational } #endif } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { - Math::Vector value; + Math::Vector value; for (auto it = mesh.getCell(); !it.end(); ++it) { const auto& polytope = *it; @@ -589,7 +612,7 @@ namespace Rodin::Variational { const Geometry::Point p(polytope, trans, fe.getNode(local)); const Index global = fes.getGlobalIndex({ d, i }, local); - fn.getValue(value, p); + fn.getDerived().getValue(value, p); m_data.col(global) = (value + ns[global] * m_data.col(global)) / (ns[global] + 1.0); ns[global] += 1; @@ -605,15 +628,59 @@ namespace Rodin::Variational } inline - auto& projectOnBoundary(std::function fn, Geometry::Attribute attr) + auto& projectOnBoundary( + std::function fn, Geometry::Attribute attr) { return projectOnBoundary(fn, FlatSet{attr}); } inline - auto& projectOnBoundary(std::function fn, const FlatSet& attrs = {}) + auto& projectOnBoundary( + std::function fn, Geometry::Attribute attr) { - return projectOnBoundary(RealFunction(fn), attrs); + return projectOnBoundary(fn, FlatSet{attr}); + } + + inline + auto& projectOnBoundary( + std::function fn, + const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return projectOnBoundary(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return projectOnBoundary(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } + } + + inline + auto& projectOnBoundary( + std::function fn, + const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return projectOnBoundary(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return projectOnBoundary(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } } template @@ -652,13 +719,13 @@ namespace Rodin::Variational { const Geometry::Point p(polytope, trans, fe.getNode(local)); const Index global = fes.getGlobalIndex({ d, i }, local); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(m_data.rows() == 1); m_data(global) = (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1.0); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { m_data.col(global) = (fn.getValue(p) + ns[global] * m_data.col(global)) / (ns[global] + 1.0); @@ -675,15 +742,57 @@ namespace Rodin::Variational } inline - auto& projectOnFaces(std::function fn, Geometry::Attribute attr) + auto& projectOnFaces( + std::function fn, Geometry::Attribute attr) { return projectOnFaces(fn, FlatSet{attr}); } inline - auto& projectOnFaces(std::function fn, const FlatSet& attrs = {}) + auto& projectOnFaces( + std::function fn, Geometry::Attribute attr) { - return projectOnFaces(RealFunction(fn), attrs); + return projectOnFaces(fn, FlatSet{attr}); + } + + inline + auto& projectOnFaces( + std::function fn, const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return projectOnFaces(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return projectOnFaces(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } + } + + inline + auto& projectOnFaces( + std::function fn, const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return projectOnFaces(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return projectOnFaces(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } } template @@ -722,13 +831,13 @@ namespace Rodin::Variational { const Geometry::Point p(polytope, trans, fe.getNode(local)); const Index global = fes.getGlobalIndex({ d, i }, local); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(m_data.rows() == 1); m_data(global) = (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1.0); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { m_data.col(global) = (fn.getValue(p) + ns[global] * m_data.col(global)) / (ns[global] + 1.0); @@ -745,15 +854,57 @@ namespace Rodin::Variational } inline - auto& projectOnInterfaces(std::function fn, Geometry::Attribute attr) + auto& projectOnInterfaces( + std::function fn, Geometry::Attribute attr) { return projectOnInterfaces(fn, FlatSet{attr}); } inline - auto& projectOnInterfaces(std::function fn, const FlatSet& attrs = {}) + auto& projectOnInterfaces( + std::function fn, Geometry::Attribute attr) { - return projectOnInterfaces(RealFunction(fn), attrs); + return projectOnInterfaces(fn, FlatSet{ attr }); + } + + inline + auto& projectOnInterfaces( + std::function fn, const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return projectOnInterfaces(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return projectOnInterfaces(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } + } + + inline + auto& projectOnInterfaces( + std::function fn, const FlatSet& attrs = {}) + { + if constexpr (std::is_same_v) + { + assert(getFiniteElementSpace().getVectorDimension() == 1); + return projectOnInterfaces(ScalarFunction(fn)); + } + else if constexpr (std::is_same_v>) + { + return projectOnInterfaces(VectorFunction(getFiniteElementSpace().getVectorDimension(), fn)); + } + else + { + assert(false); + return static_cast(*this); + } } template @@ -771,7 +922,8 @@ namespace Rodin::Variational } template - Derived& projectOnInterfaces(const FunctionBase& fn, const FlatSet& attrs) + Derived& projectOnInterfaces( + const FunctionBase& fn, const FlatSet& attrs) { using FunctionType = FunctionBase; using FunctionRangeType = typename FormLanguage::Traits::RangeType; @@ -790,12 +942,12 @@ namespace Rodin::Variational for (size_t local = 0; local < fe.getCount(); local++) { const Geometry::Point p(polytope, trans, fe.getNode(local)); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(m_data.rows() == 1); m_data(fes.getGlobalIndex({ d, i }, local)) = fn.getValue(p); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { m_data.col(fes.getGlobalIndex({ d, i }, local)) = fn.getValue(p); } @@ -905,7 +1057,7 @@ namespace Rodin::Variational */ inline constexpr - Math::Matrix& getData() + auto& getData() { return m_data; } @@ -915,21 +1067,21 @@ namespace Rodin::Variational */ inline constexpr - const Math::Matrix& getData() const + const DataType& getData() const { return m_data; } inline constexpr - std::optional>& getWeights() + std::optional& getWeights() { return m_weights; } inline constexpr - const std::optional>& getWeights() const + const std::optional& getWeights() const { return m_weights; } @@ -987,13 +1139,13 @@ namespace Rodin::Variational inline Derived& setValue(Index global, Value&& v) { - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(m_data.size() >= 0); assert(global < static_cast(m_data.size())); m_data.coeffRef(global) = std::forward(v); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { assert(m_data.cols() >= 0); assert(global < static_cast(m_data.cols())); @@ -1024,13 +1176,13 @@ namespace Rodin::Variational inline auto getValue(Index global) const { - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(m_data.size() >= 0); assert(global < static_cast(m_data.size())); return m_data.coeff(global); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { assert(m_data.cols() >= 0); assert(global < static_cast(m_data.cols())); @@ -1043,31 +1195,32 @@ namespace Rodin::Variational } } + inline + constexpr + RangeType getValue(const Geometry::Point& p) const + { + RangeType res; + getValue(res, p); + return res; + } + /** * @brief Gets the interpolated value at the point. */ inline - RangeType getValue(const Geometry::Point& p) const + void getValue(RangeType& res, const Geometry::Point& p) const { - RangeType out; - if constexpr (std::is_same_v) - out = NAN; - else if constexpr (std::is_same_v>) - out.setConstant(NAN); const auto& polytope = p.getPolytope(); const auto& polytopeMesh = polytope.getMesh(); const auto& fes = m_fes.get(); const auto& fesMesh = fes.getMesh(); if (polytopeMesh == fesMesh) { - interpolate(out, p); + interpolate(res, p); } else if (const auto inclusion = fesMesh.inclusion(p)) { - if constexpr (std::is_same_v) - out = interpolate(*inclusion); - else - interpolate(out, *inclusion); + interpolate(res, *inclusion); } else if (fesMesh.isSubMesh()) { @@ -1075,49 +1228,14 @@ namespace Rodin::Variational const auto restriction = submesh.restriction(p); if (restriction) { - if constexpr (std::is_same_v) - out = interpolate(*restriction); - else - interpolate(out, *restriction); + interpolate(res, *restriction); } else { - throw 1; + assert(false); } } else - { - assert(false); - throw 1; - } - return out; - } - - inline - constexpr - void getValue(Math::Vector& out, const Geometry::Point& p) const - { - static_assert(std::is_same_v>); - out.setConstant(NAN); - const auto& polytope = p.getPolytope(); - const auto& polytopeMesh = polytope.getMesh(); - const auto& fes = m_fes.get(); - const auto& fesMesh = fes.getMesh(); - if (polytopeMesh == fesMesh) - { - interpolate(out, p); - } - else if (const auto inclusion = fesMesh.inclusion(p)) - { - interpolate(out, *inclusion); - } - else if (fesMesh.isSubMesh()) - { - const auto& submesh = fesMesh.asSubMesh(); - const auto restriction = submesh.restriction(p); - interpolate(out, *restriction); - } - else { assert(false); } @@ -1129,10 +1247,11 @@ namespace Rodin::Variational */ inline constexpr - Real interpolate(const Geometry::Point& p) const + RangeType interpolate(const Geometry::Point& p) const { - static_assert(std::is_same_v); - return static_cast(*this).interpolate(p); + RangeType res; + res = static_cast(*this).interpolate(res, p); + return res; } /** @@ -1141,28 +1260,15 @@ namespace Rodin::Variational */ inline constexpr - void interpolate(Math::Vector& res, const Geometry::Point& p) const + void interpolate(RangeType& res, const Geometry::Point& p) const { - static_assert(std::is_same_v>); static_cast(*this).interpolate(res, p); } - /** - * @brief Deleted function. - */ - inline - constexpr - void getValue(Math::Matrix& res, const Geometry::Point& p) const = delete; - private: - void interpolate(Real& res, const Geometry::Point& p) const - { - res = interpolate(p); - } - - std::reference_wrapper m_fes; - Math::Matrix m_data; - std::optional> m_weights; + std::reference_wrapper m_fes; + DataType m_data; + std::optional m_weights; mutable Threads::Mutex m_mutex; }; } diff --git a/src/Rodin/Variational/IdentityMatrix.h b/src/Rodin/Variational/IdentityMatrix.h index 9f0c2ebc4..1940fe578 100644 --- a/src/Rodin/Variational/IdentityMatrix.h +++ b/src/Rodin/Variational/IdentityMatrix.h @@ -15,9 +15,13 @@ namespace Rodin::Variational * F(x) = I_n * @f$ */ - class IdentityMatrix : public MatrixFunctionBase + class IdentityMatrix : public MatrixFunctionBase { public: + using ScalarType = Real; + + using Parent = MatrixFunctionBase; + /** * @brief Constructs the identity matrix function. * @param[in] n Dimension of identity matrix diff --git a/src/Rodin/Variational/Jacobian.h b/src/Rodin/Variational/Jacobian.h index f152783dc..930f588b8 100644 --- a/src/Rodin/Variational/Jacobian.h +++ b/src/Rodin/Variational/Jacobian.h @@ -24,24 +24,27 @@ namespace Rodin::Variational /** * @brief Base class for Jacobian classes. */ - template + template class JacobianBase; /** * @ingroup JacobianSpecializations * @brief Jacobian of a P1 GridFunction */ - template - class JacobianBase> - : public MatrixFunctionBase>> + template + class JacobianBase, Derived> + : public MatrixFunctionBase< + typename FormLanguage::Traits::ScalarType, JacobianBase, Derived>> { public: using FESType = FES; - using OperandType = GridFunction; + using ScalarType = typename FormLanguage::Traits::ScalarType; - /// Parent class - using Parent = MatrixFunctionBase>; + using OperandType = GridFunction; + + using Parent = + MatrixFunctionBase>; /** * @brief Constructs the Jacobianient of a @f$ \mathbb{P}_1 @f$ function diff --git a/src/Rodin/Variational/LazyEvaluator.h b/src/Rodin/Variational/LazyEvaluator.h index fa6acd05f..6b844b83f 100644 --- a/src/Rodin/Variational/LazyEvaluator.h +++ b/src/Rodin/Variational/LazyEvaluator.h @@ -77,6 +77,14 @@ namespace Rodin::Variational return m_ref.get().getValue(p); } + template + inline + constexpr + void getValue(T& res, const Geometry::Point& p) const + { + m_ref.get().getValue(res, p); + } + inline LazyEvaluator* copy() const noexcept final override { return new LazyEvaluator(*this); diff --git a/src/Rodin/Variational/MatrixFunction.h b/src/Rodin/Variational/MatrixFunction.h index 262722649..b585426a4 100644 --- a/src/Rodin/Variational/MatrixFunction.h +++ b/src/Rodin/Variational/MatrixFunction.h @@ -16,6 +16,16 @@ #include "Function.h" #include "RangeShape.h" +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using ScalarType = Scalar; + using DerivedType = Derived; + }; +} + namespace Rodin::Variational { /** @@ -24,11 +34,15 @@ namespace Rodin::Variational * @see MatrixFunction */ - template - class MatrixFunctionBase : public FunctionBase> + template + class MatrixFunctionBase : public FunctionBase> { public: - using Parent = FunctionBase>; + using ScalarType = Scalar; + + using MatrixType = Math::Matrix; + + using Parent = FunctionBase>; MatrixFunctionBase() = default; @@ -55,6 +69,13 @@ namespace Rodin::Variational return static_cast(*this).getValue(p); } + inline + constexpr + void getValue(MatrixType& res, const Geometry::Point& p) const + { + return static_cast(*this).getValue(res, p); + } + inline constexpr RangeShape getRangeShape() const @@ -109,14 +130,18 @@ namespace Rodin::Variational /** * @ingroup MatrixFunctionSpecializations */ - template <> - class MatrixFunction> final - : public MatrixFunctionBase>> + template + class MatrixFunction> final + : public MatrixFunctionBase>> { public: - using Parent = MatrixFunctionBase>>; + using ScalarType = Scalar; + + using MatrixType = Math::Matrix; - MatrixFunction(std::reference_wrapper> matrix) + using Parent = MatrixFunctionBase>; + + MatrixFunction(const MatrixType& matrix) : m_matrix(matrix) {} @@ -131,11 +156,19 @@ namespace Rodin::Variational {} inline - const Math::Matrix& getValue(const Geometry::Point&) const + constexpr + const MatrixType& getValue(const Geometry::Point&) const { return m_matrix.get(); } + inline + constexpr + void getValue(MatrixType& res, const Geometry::Point&) const + { + res = m_matrix.get(); + } + inline constexpr MatrixFunction& traceOf(Geometry::Attribute) @@ -174,7 +207,7 @@ namespace Rodin::Variational } private: - std::reference_wrapper> m_matrix; + std::reference_wrapper m_matrix; }; MatrixFunction(std::reference_wrapper>) diff --git a/src/Rodin/Variational/Mult.h b/src/Rodin/Variational/Mult.h index 4d4bbf527..da20a3086 100644 --- a/src/Rodin/Variational/Mult.h +++ b/src/Rodin/Variational/Mult.h @@ -224,12 +224,12 @@ namespace Rodin::Variational { if constexpr (std::is_same_v && std::is_same_v>) { - getRHS().getValue(out, p); + getRHS().getDerived().getValue(out, p); out *= getLHS().getValue(p); } else if constexpr (std::is_same_v> && std::is_same_v) { - getLHS().getValue(out, p); + getLHS().getDerived().getValue(out, p); out *= getRHS().getValue(p); } else diff --git a/src/Rodin/Variational/Normal.h b/src/Rodin/Variational/Normal.h deleted file mode 100644 index aa2c14aae..000000000 --- a/src/Rodin/Variational/Normal.h +++ /dev/null @@ -1,63 +0,0 @@ -#ifndef RODIN_VARIATIONAL_NORMAL_H -#define RODIN_VARIATIONAL_NORMAL_H - -#include "Rodin/Geometry/Mesh.h" -#include "Rodin/Geometry/PolytopeTransformation.h" - -#include "ForwardDecls.h" -#include "VectorFunction.h" - -namespace Rodin::Variational -{ - /** - * @brief Outward unit normal. - */ - class Normal final : public VectorFunctionBase - { - public: - using Parent = VectorFunctionBase; - - /** - * @brief Constructs the outward unit normal. - */ - Normal(const Geometry::MeshBase& surface) - : m_dimension(surface.getSpaceDimension()) - { - assert(m_dimension > 0); - } - - Normal(const Normal& other) - : Parent(other), - m_dimension(other.m_dimension) - {} - - Normal(Normal&& other) - : Parent(std::move(other)), - m_dimension(std::move(other.m_dimension)) - {} - - inline - constexpr - size_t getDimension() const - { - return m_dimension; - } - - Math::SpatialVector getValue(const Geometry::Point& p) const - { - Math::SpatialVector res; - // this->getValue(res, p); - return res; - } - - inline Normal* copy() const noexcept override - { - return new Normal(*this); - } - - private: - const size_t m_dimension; - }; -} - -#endif diff --git a/src/Rodin/Variational/P0/GridFunction.h b/src/Rodin/Variational/P0/GridFunction.h index 612fb8b40..becf9989c 100644 --- a/src/Rodin/Variational/P0/GridFunction.h +++ b/src/Rodin/Variational/P0/GridFunction.h @@ -36,7 +36,7 @@ namespace Rodin::Variational */ template class GridFunction> final - : public GridFunctionBase>> + : public GridFunctionBase, GridFunction>> { public: using FESType = P0; @@ -49,8 +49,7 @@ namespace Rodin::Variational using ElementType = typename FormLanguage::Traits::ElementType; - /// Parent class - using Parent = GridFunctionBase>; + using Parent = GridFunctionBase>; using Parent::getValue; using Parent::operator=; @@ -97,7 +96,7 @@ namespace Rodin::Variational GridFunction& operator=(const GridFunction&) = delete; - Real interpolate(const Geometry::Point& p) const + void interpolate(Real& res, const Geometry::Point& p) const { static_assert(std::is_same_v); const auto& fes = this->getFiniteElementSpace(); @@ -108,7 +107,7 @@ namespace Rodin::Variational assert(d == polytope.getDimension()); const Index i = polytope.getIndex(); assert(fes.getFiniteElement(d, i).getCount() == 1); - return getValue({d, i}, 0); + res = getValue({ d, i }, 0); } // void interpolate(Math::Vector& res, const Geometry::Point& p) const diff --git a/src/Rodin/Variational/P1/Div.h b/src/Rodin/Variational/P1/Div.h index 52ffb93ae..916748362 100644 --- a/src/Rodin/Variational/P1/Div.h +++ b/src/Rodin/Variational/P1/Div.h @@ -14,22 +14,24 @@ namespace Rodin::FormLanguage { - template - struct Traits, Mesh>>>> + template + struct Traits, Mesh>>>> { - using FESType = Variational::P1, Mesh>; - using OperandType = Variational::GridFunction, Mesh>>; + using FESType = Variational::P1, Mesh>; + using ScalarType = Scalar; + using OperandType = Variational::GridFunction, Mesh>>; }; - template + template struct Traits< Variational::Div< - Variational::ShapeFunction, Mesh>, Space>>> + Variational::ShapeFunction, Mesh>, Space>>> { - using FESType = Variational::P1, Mesh>; + using FESType = Variational::P1, Mesh>; static constexpr Variational::ShapeFunctionSpaceType SpaceType = Space; + using ScalarType = Scalar; using OperandType = - Variational::ShapeFunction, Mesh>, Space>; + Variational::ShapeFunction, Mesh>, Space>; }; } @@ -39,12 +41,14 @@ namespace Rodin::Variational * @ingroup DivSpecializations * @brief Divient of a P1 GridFunction */ - template - class Div, Mesh>>> final - : public DivBase, Mesh>>>, GridFunction, Mesh>>> + template + class Div, Mesh>>> final + : public DivBase< + GridFunction, Mesh>>, + Div, Mesh>>>> { public: - using ScalarType = Number; + using ScalarType = Scalar; using FESType = Variational::P1, Mesh>; @@ -52,7 +56,7 @@ namespace Rodin::Variational using OperandType = GridFunction; /// Parent class - using Parent = DivBase, OperandType>; + using Parent = DivBase>; /** * @brief Constructs the Divient of an @f$ \mathbb{P}^1 @f$ function @@ -77,14 +81,7 @@ namespace Rodin::Variational : Parent(std::move(other)) {} - void interpolate(Math::Vector& out, const Geometry::Point& p) const - { - Math::SpatialVector tmp; - interpolate(tmp, p); - out = std::move(tmp); - } - - void interpolate(Real& out, const Geometry::Point& p) const + void interpolate(ScalarType& out, const Geometry::Point& p) const { const auto& polytope = p.getPolytope(); const auto& d = polytope.getDimension(); @@ -100,7 +97,7 @@ namespace Rodin::Variational if (inc.size() == 1) { const auto& tracePolytope = mesh.getPolytope(meshDim, *inc.begin()); - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -125,7 +122,7 @@ namespace Rodin::Variational const auto& tracePolytope = mesh.getPolytope(meshDim, idx); if (traceDomain.count(tracePolytope->getAttribute())) { - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -145,7 +142,7 @@ namespace Rodin::Variational const auto& vdim = fes.getVectorDimension(); const auto& fe = fes.getFiniteElement(d, i); const auto& rc = p.getReferenceCoordinates(); - Math::SpatialMatrix jacobian(vdim, d); + Math::SpatialMatrix jacobian(vdim, d); out = 0; for (size_t local = 0; local < fe.getCount(); local++) { diff --git a/src/Rodin/Variational/P1/Grad.h b/src/Rodin/Variational/P1/Grad.h index df5576c3c..eae629575 100644 --- a/src/Rodin/Variational/P1/Grad.h +++ b/src/Rodin/Variational/P1/Grad.h @@ -50,18 +50,20 @@ namespace Rodin::Variational * @ingroup GradSpecializations * @brief Gradient of a P1 GridFunction */ - template - class Grad>> final - : public GradBase>>, GridFunction>> + template + class Grad>> final + : public GradBase>, Grad>>> { public: - using FESType = P1; + using FESType = P1; + + using RangeType = Range; + + using ScalarType = typename FormLanguage::Traits::ScalarType; - /// Operand type using OperandType = GridFunction; - /// Parent class - using Parent = GradBase, OperandType>; + using Parent = GradBase>; /** * @brief Constructs the gradient of an @f$ \mathbb{P}^1 @f$ function diff --git a/src/Rodin/Variational/P1/GridFunction.h b/src/Rodin/Variational/P1/GridFunction.h index 64190d4d4..28da7556f 100644 --- a/src/Rodin/Variational/P1/GridFunction.h +++ b/src/Rodin/Variational/P1/GridFunction.h @@ -36,7 +36,7 @@ namespace Rodin::Variational */ template class GridFunction> final - : public GridFunctionBase>> + : public GridFunctionBase, GridFunction>> { public: /// Type of finite element space to which the GridFunction belongs to @@ -51,7 +51,7 @@ namespace Rodin::Variational using ElementType = typename FormLanguage::Traits::ElementType; /// Parent class - using Parent = GridFunctionBase>; + using Parent = GridFunctionBase>; using Parent::getValue; using Parent::operator=; @@ -98,7 +98,7 @@ namespace Rodin::Variational GridFunction& operator=(const GridFunction&) = delete; - Real interpolate(const Geometry::Point& p) const + void interpolate(Real& res, const Geometry::Point& p) const { static_assert(std::is_same_v); const auto& fes = this->getFiniteElementSpace(); @@ -109,13 +109,12 @@ namespace Rodin::Variational const Index i = polytope.getIndex(); const auto& fe = fes.getFiniteElement(d, i); const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); - Real res = 0; + res = 0; for (Index local = 0; local < fe.getCount(); local++) { const auto& basis = fe.getBasis(local); res += getValue({d, i}, local) * basis(r); } - return res; } void interpolate(Math::Vector& res, const Geometry::Point& p) const diff --git a/src/Rodin/Variational/P1/Jacobian.h b/src/Rodin/Variational/P1/Jacobian.h index bdd445f2a..6321f2e8d 100644 --- a/src/Rodin/Variational/P1/Jacobian.h +++ b/src/Rodin/Variational/P1/Jacobian.h @@ -45,14 +45,19 @@ namespace Rodin::Variational template class Jacobian, Mesh>>> final : public JacobianBase< - Jacobian, Mesh>>>, GridFunction, Mesh>>> + GridFunction, Mesh>>, + Jacobian, Mesh>>>> { public: using FESType = P1, Mesh>; + using ScalarType = typename FormLanguage::Traits; + + using MatrixType = Math::SpatialMatrix; + using OperandType = GridFunction; - using Parent = JacobianBase, OperandType>; + using Parent = JacobianBase>; /** * @brief Constructs the Jacobian matrix of an @f$ H^1 (\Omega)^d @f$ function diff --git a/src/Rodin/Variational/P1/P1Element.cpp b/src/Rodin/Variational/P1/P1Element.cpp index 351677031..38212d2e0 100644 --- a/src/Rodin/Variational/P1/P1Element.cpp +++ b/src/Rodin/Variational/P1/P1Element.cpp @@ -279,7 +279,7 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } @@ -302,7 +302,7 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } @@ -331,7 +331,7 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } @@ -358,13 +358,13 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } } assert(false); - return NAN; + return Math::nan(); } void @@ -396,7 +396,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -426,7 +426,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -463,7 +463,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -502,14 +502,14 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } } } assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); } Complex ComplexP1Element::BasisFunction::operator()(const Math::SpatialVector& r) const @@ -519,10 +519,22 @@ namespace Rodin::Variational { case Geometry::Polytope::Type::Point: { - if (m_i % 2 == 0) - return 1; - else - return 1i; + switch (m_i) + { + case 0: + { + return 1; + } + case 1: + { + return 1i; + } + default: + { + assert(false); + return Math::nan(); + } + } } case Geometry::Polytope::Type::Segment: { @@ -547,7 +559,7 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } @@ -582,7 +594,7 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } @@ -621,7 +633,7 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } @@ -664,13 +676,13 @@ namespace Rodin::Variational default: { assert(false); - return NAN; + return Math::nan(); } } } } assert(false); - return NAN; + return Math::nan(); } void @@ -713,7 +725,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -760,7 +772,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -821,7 +833,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -886,7 +898,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -953,7 +965,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -985,7 +997,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -1017,7 +1029,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -1056,7 +1068,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -1083,7 +1095,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -1115,7 +1127,7 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } @@ -1147,13 +1159,13 @@ namespace Rodin::Variational default: { assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); return; } } } } assert(false); - out.setConstant(NAN); + out.setConstant(Math::nan()); } } diff --git a/src/Rodin/Variational/P1/QuadratureRule.h b/src/Rodin/Variational/P1/QuadratureRule.h index 67a194320..b2e8b68d0 100644 --- a/src/Rodin/Variational/P1/QuadratureRule.h +++ b/src/Rodin/Variational/P1/QuadratureRule.h @@ -250,7 +250,7 @@ namespace Rodin::Variational } else if constexpr (std::is_same_v, LHSRangeType>) { - f.getValue(m_vv, p); + f.getDerived().getValue(m_vv, p); for (size_t i = 0; i < fe.getCount(); i++) { fe.getBasis(i)(m_vb, rc); diff --git a/src/Rodin/Variational/RealFunction.h b/src/Rodin/Variational/RealFunction.h index 53a1f922c..6ab41fcfd 100644 --- a/src/Rodin/Variational/RealFunction.h +++ b/src/Rodin/Variational/RealFunction.h @@ -4,8 +4,8 @@ * (See accompanying file LICENSE or copy at * https://www.boost.org/LICENSE_1_0.txt) */ -#ifndef RODIN_VARIATIONAL_SCALARFUNCTION_H -#define RODIN_VARIATIONAL_SCALARFUNCTION_H +#ifndef RODIN_VARIATIONAL_REALFUNCTION_H +#define RODIN_VARIATIONAL_REALFUNCTION_H #include #include @@ -21,6 +21,15 @@ #include "ScalarFunction.h" +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using ScalarType = Real; + using DerivedType = Derived; + }; +} namespace Rodin::Variational { @@ -72,17 +81,13 @@ namespace Rodin::Variational * @ingroup RealFunctionSpecializations */ template - class RealFunction> final - : public RealFunctionBase> + class RealFunction> final + : public RealFunctionBase> { public: - using Parent = RealFunctionBase>; - - using NestedRangeType = typename FormLanguage::Traits>::RangeType; - - static_assert(std::is_same_v); + using Parent = RealFunctionBase>; - RealFunction(const FunctionBase& nested) + RealFunction(const RealFunctionBase& nested) : m_nested(nested.copy()) {} @@ -117,14 +122,14 @@ namespace Rodin::Variational } private: - std::unique_ptr> m_nested; + std::unique_ptr> m_nested; }; /** * @brief CTAD for RealFunction. */ template - RealFunction(const FunctionBase&) -> RealFunction>; + RealFunction(const RealFunctionBase&) -> RealFunction>; /** * @ingroup RealFunctionSpecializations diff --git a/src/Rodin/Variational/ScalarFunction.h b/src/Rodin/Variational/ScalarFunction.h index 7b7cef892..471f4a5a5 100644 --- a/src/Rodin/Variational/ScalarFunction.h +++ b/src/Rodin/Variational/ScalarFunction.h @@ -4,8 +4,8 @@ * (See accompanying file LICENSE or copy at * https://www.boost.org/LICENSE_1_0.txt) */ -#ifndef RODIN_VARIATIONAL_NUMBERFUNCTION_H -#define RODIN_VARIATIONAL_NUMBERFUNCTION_H +#ifndef RODIN_VARIATIONAL_SCALARFUNCTION_H +#define RODIN_VARIATIONAL_SCALARFUNCTION_H #include #include @@ -17,24 +17,34 @@ #include "ForwardDecls.h" -#include "Function.h" #include "RangeShape.h" +#include "Function.h" + +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using ScalarType = Scalar; + using DerivedType = Derived; + }; +} + namespace Rodin::Variational { /** - * @defgroup ScalarFunctionSpecializations ScalarFunction Template Specializations - * @brief Template specializations of the ScalarFunction class. - * @see ScalarFunction + * @defgroup RealFunctionSpecializations RealFunction Template Specializations + * @brief Template specializations of the RealFunction class. + * @see RealFunction */ - template - class ScalarFunctionBase : public FunctionBase> + template + class ScalarFunctionBase + : public FunctionBase> { public: - using ScalarType = Number; - - using RangeType = ScalarType; + using ScalarType = Scalar; using Parent = FunctionBase>; @@ -67,24 +77,262 @@ namespace Rodin::Variational inline constexpr - void getValue(Math::Vector&, const Geometry::Point&) const = delete; + void getValue(ScalarType& res, const Geometry::Point& p) const + { + res = this->getValue(p); + } + + virtual ScalarFunctionBase* copy() const noexcept override = 0; + }; + + /** + * @ingroup ScalarFunctionSpecializations + */ + template + class ScalarFunction> final + : public ScalarFunctionBase>> + { + public: + using ScalarType = Scalar; + + using Parent = + ScalarFunctionBase>>; + + ScalarFunction(const ScalarFunctionBase& nested) + : m_nested(nested.copy()) + {} + + ScalarFunction(const ScalarFunction& other) + : Parent(other), + m_nested(other.m_nested->copy()) + {} + + ScalarFunction(ScalarFunction&& other) + : Parent(std::move(other)), + m_nested(std::move(other.m_nested)) + {} inline constexpr - void getValue(Math::Matrix&, const Geometry::Point&) const = delete; + auto getValue(const Geometry::Point& v) const + { + return m_nested->getValue(v); + } inline constexpr - RangeShape getRangeShape() const + ScalarFunction& traceOf(Geometry::Attribute attrs) { - return { 1, 1 }; + m_nested->traceOf(attrs); + return *this; } - virtual ScalarFunctionBase* copy() const noexcept override + inline ScalarFunction* copy() const noexcept override { - return static_cast(*this).copy(); + return new ScalarFunction(*this); } + + private: + std::unique_ptr> m_nested; }; + + /** + * @brief CTAD for ScalarFunction. + */ + template + ScalarFunction(const ScalarFunctionBase&) + -> ScalarFunction>; + + /** + * @ingroup ScalarFunctionSpecializations + * @brief Represents a constant scalar function with type Real. + */ + template <> + class ScalarFunction final + : public ScalarFunctionBase> + { + public: + using ScalarType = Real; + + using Parent = ScalarFunctionBase>; + + /** + * @brief Constructs a ScalarFunction from a constant scalar value. + * @param[in] x Constant scalar value + */ + ScalarFunction(const ScalarType& x) + : m_x(x) + {} + + ScalarFunction(const ScalarFunction& other) + : Parent(other), + m_x(other.m_x) + {} + + ScalarFunction(ScalarFunction&& other) + : Parent(std::move(other)), + m_x(other.m_x) + {} + + inline + constexpr + ScalarFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + ScalarType getValue() const + { + return m_x; + } + + inline + constexpr + ScalarType getValue(const Geometry::Point&) const + { + return m_x; + } + + inline ScalarFunction* copy() const noexcept override + { + return new ScalarFunction(*this); + } + + private: + const ScalarType m_x; + }; + + /** + * @brief CTAD for ScalarFunction. + */ + ScalarFunction(const Real&) -> ScalarFunction; + + /** + * @ingroup ScalarFunctionSpecializations + * @brief Represents a constant scalar function with type Complex. + */ + template <> + class ScalarFunction final + : public ScalarFunctionBase> + { + public: + using ScalarType = Complex; + + using Parent = ScalarFunctionBase>; + + /** + * @brief Constructs a ScalarFunction from a constant scalar value. + * @param[in] x Constant scalar value + */ + ScalarFunction(const ScalarType& x) + : m_x(x) + {} + + ScalarFunction(const ScalarFunction& other) + : Parent(other), + m_x(other.m_x) + {} + + ScalarFunction(ScalarFunction&& other) + : Parent(std::move(other)), + m_x(other.m_x) + {} + + inline + constexpr + ScalarFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + ScalarType getValue() const + { + return m_x; + } + + inline + constexpr + ScalarType getValue(const Geometry::Point&) const + { + return m_x; + } + + inline ScalarFunction* copy() const noexcept override + { + return new ScalarFunction(*this); + } + + private: + const ScalarType m_x; + }; + + /** + * @brief CTAD for ScalarFunction. + */ + ScalarFunction(const Complex&) -> ScalarFunction; + + /** + * @ingroup ScalarFunctionSpecializations + * @brief Represents a scalar function given by an arbitrary scalar function. + */ + template + class ScalarFunction final + : public ScalarFunctionBase, ScalarFunction> + { + + public: + using ScalarType = std::invoke_result_t; + + using Parent = ScalarFunctionBase>; + + static_assert(std::is_invocable_v); + + ScalarFunction(F f) + : m_f(f) + {} + + ScalarFunction(const ScalarFunction& other) + : Parent(other), + m_f(other.m_f) + {} + + ScalarFunction(ScalarFunction&& other) + : Parent(std::move(other)), + m_f(std::move(other.m_f)) + {} + + inline + constexpr + ScalarFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + ScalarType getValue(const Geometry::Point& v) const + { + return m_f(v); + } + + inline ScalarFunction* copy() const noexcept override + { + return new ScalarFunction(*this); + } + + private: + const F m_f; + }; + + /** + * @brief CTAD for ScalarFunction. + */ + template >> + ScalarFunction(F) -> ScalarFunction; } #endif diff --git a/src/Rodin/Variational/VectorFunction.h b/src/Rodin/Variational/VectorFunction.h index 1589c8496..743257291 100644 --- a/src/Rodin/Variational/VectorFunction.h +++ b/src/Rodin/Variational/VectorFunction.h @@ -19,6 +19,16 @@ #include "Function.h" #include "RealFunction.h" +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using ScalarType = Scalar; + using DerivedType = Derived; + }; +} + namespace Rodin::Variational { /** @@ -30,11 +40,15 @@ namespace Rodin::Variational /** * @brief Base class for vector-valued functions defined on a mesh. */ - template - class VectorFunctionBase : public FunctionBase> + template + class VectorFunctionBase : public FunctionBase> { public: - using Parent = FunctionBase>; + using ScalarType = Scalar; + + using VectorType = Math::Vector; + + using Parent = FunctionBase>; VectorFunctionBase() = default; @@ -139,9 +153,9 @@ namespace Rodin::Variational inline constexpr - void getValue(Math::Vector& res, const Geometry::Point& p) const + void getValue(VectorType& res, const Geometry::Point& p) const { - if constexpr (Internal::HasGetValueMethod&>::Value) + if constexpr (Internal::HasGetValueMethod::Value) { return static_cast(*this).getValue(res, p); } @@ -151,10 +165,6 @@ namespace Rodin::Variational } } - inline - constexpr - void getValue(Math::Matrix& res, const Geometry::Point& p) const = delete; - /** * @brief Gets the dimension of the vector object. * @returns Dimension of vector. @@ -172,12 +182,16 @@ namespace Rodin::Variational } }; - template <> - class VectorFunction> final - : public VectorFunctionBase>> + template + class VectorFunction> final + : public VectorFunctionBase>> { public: - using Parent = VectorFunctionBase>>; + using ScalarType = Scalar; + + using VectorType = Math::Vector; + + using Parent = VectorFunctionBase>; /** * @brief Constructs a vector with the given values. @@ -186,31 +200,38 @@ namespace Rodin::Variational * Each value passed must be convertible to any specialization of * RealFunction. */ - VectorFunction(const Math::Vector& v) - : m_value(v) + VectorFunction(const VectorType& v) + : m_vector(v) {} VectorFunction(const VectorFunction& other) : Parent(other), - m_value(other.m_value) + m_vector(other.m_vector) {} VectorFunction(VectorFunction&& other) : Parent(std::move(other)), - m_value(std::move(other.m_value)) + m_vector(std::move(other.m_vector)) {} inline - const Math::Vector& getValue(const Geometry::Point& p) const + const VectorType& getValue(const Geometry::Point& p) const { - return m_value.get(); + return m_vector.get(); + } + + inline + constexpr + void getValue(VectorType& res, const Geometry::Point&) const + { + res = m_vector.get(); } inline constexpr size_t getDimension() const { - return m_value.get().size(); + return m_vector.get().size(); } inline @@ -226,10 +247,11 @@ namespace Rodin::Variational } private: - std::reference_wrapper> m_value; + std::reference_wrapper m_vector; }; - VectorFunction(const Math::Vector&) -> VectorFunction>; + template + VectorFunction(const Math::Vector&) -> VectorFunction>; /** * @ingroup VectorFunctionSpecializations @@ -253,10 +275,16 @@ namespace Rodin::Variational */ template class VectorFunction final - : public VectorFunctionBase> + : public VectorFunctionBase> { public: - using Parent = VectorFunctionBase>; + using ScalarType = Real; + + using VectorType = Math::Vector; + + using FixedSizeVectorType = Math::FixedSizeVector; + + using Parent = VectorFunctionBase>; /** * @brief Constructs a vector with the given values. * @param[in] values Parameter pack of values @@ -281,20 +309,27 @@ namespace Rodin::Variational inline auto getValue(const Geometry::Point& p) const { - Math::FixedSizeVector res; + Math::FixedSizeVector res; Utility::ForIndex<1 + sizeof...(Values)>( [&](auto i){ res.coeffRef(static_cast(i)) = std::get(m_fs).getValue(p); }); return res; } inline - void getValue(Math::Vector& res, const Geometry::Point& p) const + void getValue(VectorType& res, const Geometry::Point& p) const { res.resize(1 + sizeof...(Values)); Utility::ForIndex<1 + sizeof...(Values)>( [&](auto i){ res.coeffRef(static_cast(i)) = std::get(m_fs).getValue(p); }); } + inline + void getValue(FixedSizeVectorType& res, const Geometry::Point& p) const + { + Utility::ForIndex<1 + sizeof...(Values)>( + [&](auto i){ res.coeffRef(static_cast(i)) = std::get(m_fs).getValue(p); }); + } + inline constexpr size_t getDimension() const @@ -323,14 +358,14 @@ namespace Rodin::Variational VectorFunction(const V&, const Values&...) -> VectorFunction; template - class VectorFunction final : public VectorFunctionBase> + class VectorFunction final : public VectorFunctionBase> { - static_assert( - std::is_invocable_r_v, F, const Geometry::Point&> - || std::is_invocable_v&, const Geometry::Point&>); - public: - using Parent = VectorFunctionBase>; + using ScalarType = Real; + + using VectorType = Math::Vector; + + using Parent = VectorFunctionBase>; VectorFunction(size_t vdim, F f) : m_vdim(vdim), m_f(f) @@ -356,41 +391,41 @@ namespace Rodin::Variational } inline - Math::Vector getValue(const Geometry::Point& p) const + VectorType getValue(const Geometry::Point& p) const { - if constexpr (std::is_invocable_r_v, F, const Geometry::Point&>) + if constexpr (std::is_invocable_r_v) { return m_f(p); } - else if constexpr (std::is_invocable_r_v&, const Geometry::Point&>) + else if constexpr (std::is_invocable_r_v) { - Math::Vector res; + VectorType res; m_f(res, p); return res; } else { assert(false); - Math::Vector res; + VectorType res; res.setConstant(NAN); } } inline - void getValue(Math::Vector& res, const Geometry::Point& p) const + void getValue(VectorType& res, const Geometry::Point& p) const { - if constexpr (std::is_invocable_r_v, F, const Geometry::Point&>) + if constexpr (std::is_invocable_r_v) { res = m_f(p); } - else if constexpr (std::is_invocable_v&, const Geometry::Point&>) + else if constexpr (std::is_invocable_v) { m_f(res, p); } else { assert(false); - Math::Vector res; + VectorType res; res.setConstant(NAN); } } @@ -405,9 +440,10 @@ namespace Rodin::Variational const F m_f; }; - template , F, const Geometry::Point&> - || std::is_invocable_v&, const Geometry::Point&>>> + template , F, const Geometry::Point&> || + std::is_invocable_v&, const Geometry::Point&>>> VectorFunction(size_t, F) -> VectorFunction; } diff --git a/src/Rodin/Variational/Zero.h b/src/Rodin/Variational/Zero.h index 4ef6e2673..9b58af965 100644 --- a/src/Rodin/Variational/Zero.h +++ b/src/Rodin/Variational/Zero.h @@ -26,13 +26,14 @@ namespace Rodin::Variational * @see Zero */ - template <> - class Zero final - : public RealFunctionBase> + template + class Zero final + : public ScalarFunctionBase> { public: - /// Parent class - using Parent = RealFunctionBase>; + using ScalarType = Scalar; + + using Parent = ScalarFunctionBase>; Zero() {} @@ -67,15 +68,16 @@ namespace Rodin::Variational /** * @brief CTAD for Zero. */ - Zero() -> Zero; + Zero() -> Zero; - template <> - class Zero final - : public VectorFunctionBase> + template + class Zero> final + : public VectorFunctionBase>> { public: - /// Parent class - using Parent = VectorFunctionBase>; + using VectorType = Math::Vector; + + using Parent = VectorFunctionBase>; Zero(size_t d) : m_d(d) @@ -101,11 +103,11 @@ namespace Rodin::Variational inline auto getValue(const Geometry::Point&) const { - return Math::Vector::Zero(m_d); + return VectorType::Zero(m_d); } inline - void getValue(Math::Vector& out, const Geometry::Point&) const + void getValue(VectorType& out, const Geometry::Point&) const { out.resize(m_d); out.setZero(); @@ -120,9 +122,10 @@ namespace Rodin::Variational const size_t m_d; }; - Zero(size_t) -> Zero; + Zero(size_t) -> Zero>; - using VectorZero = Zero; + template + using VectorZero = Zero>; } #endif From 078b94711ee002efaf47c80a28bddf68653ec868 Mon Sep 17 00:00:00 2001 From: Carlos Brito Date: Mon, 8 Jul 2024 01:13:05 +0200 Subject: [PATCH 2/4] CI --- src/Rodin/Variational/Function.h | 74 ++++++-------------------- src/Rodin/Variational/MatrixFunction.h | 9 +++- src/Rodin/Variational/Mult.h | 38 +++++++------ src/Rodin/Variational/ScalarFunction.h | 7 +++ src/Rodin/Variational/Sum.h | 43 +++++++++------ src/Rodin/Variational/VectorFunction.h | 2 +- 6 files changed, 76 insertions(+), 97 deletions(-) diff --git a/src/Rodin/Variational/Function.h b/src/Rodin/Variational/Function.h index f61c50957..bf59e89ef 100644 --- a/src/Rodin/Variational/Function.h +++ b/src/Rodin/Variational/Function.h @@ -158,38 +158,6 @@ namespace Rodin::Variational return static_cast(*this).getRangeShape(); } - // inline - // constexpr - // RangeType getRangeType() const - // { - // using R = typename FormLanguage::Traits>::RangeType; - // if constexpr (std::is_same_v) - // { - // return RangeType::Boolean; - // } - // else if constexpr (std::is_same_v) - // { - // return RangeType::Integer; - // } - // else if constexpr (std::is_same_v) - // { - // return RangeType::Real; - // } - // else if constexpr (Utility::IsSpecialization) - // { - // return RangeType::Vector; - // } - // else if constexpr (std::is_same_v>) - // { - // return RangeType::Matrix; - // } - // else - // { - // assert(false); - // static_assert(Utility::DependentFalse::Value); - // } - // } - /** * @brief Evaluates the function on a Point belonging to the mesh. * @note CRTP function to be overriden in Derived class. @@ -201,33 +169,21 @@ namespace Rodin::Variational return static_cast(*this).getValue(p); } - // inline - // constexpr - // void getValue(Math::Vector& res, const Geometry::Point& p) const - // { - // if constexpr (Internal::HasGetValueMethod&, const Geometry::Point&>::Value) - // { - // return static_cast(*this).getValue(res, p); - // } - // else - // { - // res = getValue(p); - // } - // } - - // inline - // constexpr - // void getValue(Math::Matrix& res, const Geometry::Point& p) const - // { - // if constexpr (Internal::HasGetValueMethod&, const Geometry::Point&>::Value) - // { - // return static_cast(*this).getValue(res, p); - // } - // else - // { - // res = getValue(p); - // } - // } + template + inline + constexpr + auto getValue(T& res, const Geometry::Point& p) const + { + if constexpr ( + Internal::HasGetValueMethod::Value) + { + return static_cast(*this).getValue(res, p); + } + else + { + res = getValue(p); + } + } inline auto coeff(size_t i, size_t j) const diff --git a/src/Rodin/Variational/MatrixFunction.h b/src/Rodin/Variational/MatrixFunction.h index b585426a4..62660a17b 100644 --- a/src/Rodin/Variational/MatrixFunction.h +++ b/src/Rodin/Variational/MatrixFunction.h @@ -73,7 +73,14 @@ namespace Rodin::Variational constexpr void getValue(MatrixType& res, const Geometry::Point& p) const { - return static_cast(*this).getValue(res, p); + if constexpr (Internal::HasGetValueMethod::Value) + { + return static_cast(*this).getValue(res, p); + } + else + { + res = getValue(p); + } } inline diff --git a/src/Rodin/Variational/Mult.h b/src/Rodin/Variational/Mult.h index da20a3086..fa384f521 100644 --- a/src/Rodin/Variational/Mult.h +++ b/src/Rodin/Variational/Mult.h @@ -34,6 +34,8 @@ namespace Rodin::FormLanguage using FESType = FES; static constexpr Variational::ShapeFunctionSpaceType SpaceType = Space; + using ScalarType = typename FormLanguage::Traits::ScalarType; + using LHSType = Variational::FunctionBase; @@ -49,32 +51,32 @@ namespace Rodin::FormLanguage using RangeType = std::conditional_t< // If - std::is_same_v, - // Then <----------------------------------------------- LHS is Real + std::is_same_v, + // Then <----------------------------------------------- LHS is Scalar RHSRangeType, // ------------------------------------------------------------------- // Else std::conditional_t< // If - std::is_same_v>, + std::is_same_v>, // Then <--------------------------------------------- LHS is Vector std::conditional_t< // If - std::is_same_v, + std::is_same_v, // Then - Math::Vector, + Math::Vector, // Else std::conditional_t< // If - std::is_same_v>, + std::is_same_v>, // Then void, // Else std::conditional_t< // If - std::is_same_v>, + std::is_same_v>, // Then - Math::Matrix, + Math::Matrix, // Else void > @@ -84,17 +86,17 @@ namespace Rodin::FormLanguage // Else std::conditional_t< // If - std::is_same_v>, + std::is_same_v>, // Then <------------------------------------------- LHS is Matrix std::conditional_t< - std::is_same_v, - Math::Matrix, + std::is_same_v, + Math::Matrix, std::conditional_t< - std::is_same_v>, - Math::Vector, + std::is_same_v>, + Math::Vector, std::conditional_t< - std::is_same_v>, - Math::Matrix, + std::is_same_v>, + Math::Matrix, void > > @@ -147,11 +149,7 @@ namespace Rodin::Variational Mult(const LHSType& lhs, const RHSType& rhs) : m_lhs(lhs.copy()), m_rhs(rhs.copy()) - { - assert(lhs.getRangeType() == RangeType::Real - || rhs.getRangeType() == RangeType::Real - || lhs.getRangeShape().width() == rhs.getRangeShape().height()); - } + {} Mult(const Mult& other) : Parent(other), diff --git a/src/Rodin/Variational/ScalarFunction.h b/src/Rodin/Variational/ScalarFunction.h index 471f4a5a5..41667c620 100644 --- a/src/Rodin/Variational/ScalarFunction.h +++ b/src/Rodin/Variational/ScalarFunction.h @@ -62,6 +62,13 @@ namespace Rodin::Variational virtual ~ScalarFunctionBase() = default; + inline + constexpr + RangeShape getRangeShape() const + { + return { 1, 1 }; + } + inline const Derived& getDerived() const { diff --git a/src/Rodin/Variational/Sum.h b/src/Rodin/Variational/Sum.h index 672360bde..4403df37f 100644 --- a/src/Rodin/Variational/Sum.h +++ b/src/Rodin/Variational/Sum.h @@ -14,6 +14,9 @@ #include "Function.h" #include "ShapeFunction.h" +#include "RealFunction.h" +#include "ComplexFunction.h" + #include "LinearFormIntegrator.h" #include "BilinearFormIntegrator.h" @@ -153,20 +156,11 @@ namespace Rodin::Variational return this->object(getLHS().getValue(p)) + this->object(getRHS().getValue(p)); } + template inline constexpr - void getValue(Math::Vector& res, const Geometry::Point& p) const - { - static_assert(FormLanguage::IsVectorRange::Value); - getLHS().getValue(res, p); - res += getRHS().getValue(p); - } - - inline - constexpr - void getValue(Math::Matrix& res, const Geometry::Point& p) const + void getValue(T& res, const Geometry::Point& p) const { - static_assert(FormLanguage::IsMatrixRange::Value); getLHS().getValue(res, p); res += getRHS().getValue(p); } @@ -194,24 +188,42 @@ namespace Rodin::Variational return Sum(lhs, rhs); } - template >> + template inline constexpr auto - operator+(const FunctionBase& lhs, Number rhs) + operator+(const FunctionBase& lhs, Real rhs) { return Sum(lhs, RealFunction(rhs)); } - template >> + template inline constexpr auto - operator+(Number lhs, const FunctionBase& rhs) + operator+(Real lhs, const FunctionBase& rhs) { return Sum(RealFunction(lhs), rhs); } + template + inline + constexpr + auto + operator+(const FunctionBase& lhs, Complex rhs) + { + return Sum(lhs, ComplexFunction(rhs)); + } + + template + inline + constexpr + auto + operator+(Complex lhs, const FunctionBase& rhs) + { + return Sum(ComplexFunction(lhs), rhs); + } + /** * @ingroup SumSpecializations */ @@ -241,7 +253,6 @@ namespace Rodin::Variational { assert(lhs.getRangeShape() == rhs.getRangeShape()); assert(lhs.getLeaf().getUUID() == rhs.getLeaf().getUUID()); - assert(lhs.getFiniteElementSpace() == rhs.getFiniteElementSpace()); } constexpr diff --git a/src/Rodin/Variational/VectorFunction.h b/src/Rodin/Variational/VectorFunction.h index 743257291..57e934b9c 100644 --- a/src/Rodin/Variational/VectorFunction.h +++ b/src/Rodin/Variational/VectorFunction.h @@ -155,7 +155,7 @@ namespace Rodin::Variational constexpr void getValue(VectorType& res, const Geometry::Point& p) const { - if constexpr (Internal::HasGetValueMethod::Value) + if constexpr (Internal::HasGetValueMethod::Value) { return static_cast(*this).getValue(res, p); } From f9e84abf9a0486d7cce444271987f600e5f0dd47 Mon Sep 17 00:00:00 2001 From: Carlos Brito Date: Tue, 9 Jul 2024 17:06:54 +0200 Subject: [PATCH 3/4] Helmholtz example --- examples/PDEs/Helmholtz.cpp | 7 +- src/Rodin/Assembly/Multithreaded.h | 203 ++++++++-- src/Rodin/Assembly/Sequential.h | 292 +++++++++------ src/Rodin/FormLanguage/Traits.h | 30 ++ src/Rodin/Math/Common.h | 9 + src/Rodin/Math/Kernels.h | 24 +- src/Rodin/Variational.h | 3 + src/Rodin/Variational/BilinearForm.h | 34 +- src/Rodin/Variational/BooleanFunction.h | 8 - src/Rodin/Variational/ComplexFunction.h | 294 ++++++++++++--- src/Rodin/Variational/Cosine.h | 30 +- src/Rodin/Variational/DenseProblem.h | 52 ++- src/Rodin/Variational/Div.h | 11 +- src/Rodin/Variational/ForwardDecls.h | 7 +- src/Rodin/Variational/Grad.h | 4 +- src/Rodin/Variational/GridFunction.h | 57 ++- src/Rodin/Variational/Im.h | 63 ++++ src/Rodin/Variational/Jacobian.h | 13 +- src/Rodin/Variational/LinearForm.h | 4 +- src/Rodin/Variational/MatrixFunction.h | 3 +- src/Rodin/Variational/Max.h | 6 +- src/Rodin/Variational/Min.h | 12 +- src/Rodin/Variational/P0/ForwardDecls.h | 5 +- src/Rodin/Variational/P0/GridFunction.h | 37 +- src/Rodin/Variational/P1/Div.h | 12 +- src/Rodin/Variational/P1/Grad.h | 19 +- src/Rodin/Variational/P1/GridFunction.h | 55 +-- src/Rodin/Variational/P1/Jacobian.h | 14 +- src/Rodin/Variational/P1/P1.cpp | 20 + src/Rodin/Variational/P1/P1.h | 259 +++++++++++-- src/Rodin/Variational/P1/P1Element.cpp | 207 +++++++++-- src/Rodin/Variational/P1/P1Element.h | 11 +- src/Rodin/Variational/Pow.h | 13 +- src/Rodin/Variational/Problem.h | 49 ++- src/Rodin/Variational/ProblemBody.h | 470 +++++++++++++----------- src/Rodin/Variational/QuadratureRule.h | 2 +- src/Rodin/Variational/RangeType.cpp | 7 +- src/Rodin/Variational/RangeType.h | 1 + src/Rodin/Variational/Re.h | 62 ++++ src/Rodin/Variational/RealFunction.h | 32 +- src/Rodin/Variational/ScalarFunction.h | 117 +++--- src/Rodin/Variational/Sine.h | 30 +- src/Rodin/Variational/Sqrt.h | 4 +- src/Rodin/Variational/Tangent.h | 33 +- src/Rodin/Variational/Traits.h | 52 ++- src/Rodin/Variational/VectorFunction.h | 3 +- 46 files changed, 1871 insertions(+), 809 deletions(-) create mode 100644 src/Rodin/Variational/Im.h create mode 100644 src/Rodin/Variational/Re.h diff --git a/examples/PDEs/Helmholtz.cpp b/examples/PDEs/Helmholtz.cpp index e6d94a64b..a1840c3ff 100644 --- a/examples/PDEs/Helmholtz.cpp +++ b/examples/PDEs/Helmholtz.cpp @@ -23,9 +23,10 @@ int main(int, char**) P1 vh(mesh); GridFunction gf(vh); - Math::Vector m(2); - // m << Complex(3, 1), Complex(2, 9); - m.setOnes(); + TrialFunction u(vh); + TestFunction v(vh); + + Problem helmholtz(u, v); // TrialFunction uRe(vh), uIm(vh); // TestFunction vRe(vh), vIm(vh); diff --git a/src/Rodin/Assembly/Multithreaded.h b/src/Rodin/Assembly/Multithreaded.h index 0cfe0340f..cfcb7a38b 100644 --- a/src/Rodin/Assembly/Multithreaded.h +++ b/src/Rodin/Assembly/Multithreaded.h @@ -53,14 +53,31 @@ namespace Rodin::Assembly template class Multithreaded< - std::vector>, - Variational::BilinearForm>>> + std::vector::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>> final : public AssemblyBase< - std::vector>, - Variational::BilinearForm>>> + std::vector::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>> { public: - using ScalarType = Real; + using ScalarType = + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>()); using OperatorType = std::vector>; @@ -271,39 +288,96 @@ namespace Rodin::Assembly template thread_local - std::vector> + std::vector< + Eigen::Triplet< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>> Multithreaded< - std::vector>, - Variational::BilinearForm>>>::tl_triplets; + std::vector< + Eigen::Triplet< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>>::tl_triplets; template thread_local - std::unique_ptr> + std::unique_ptr< + Variational::LocalBilinearFormIntegratorBase< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>> Multithreaded< - std::vector>, - Variational::BilinearForm>>>::tl_lbfi; + std::vector< + Eigen::Triplet< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>>::tl_lbfi; template thread_local - std::unique_ptr> + std::unique_ptr< + Variational::GlobalBilinearFormIntegratorBase< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>> Multithreaded< - std::vector>, - Variational::BilinearForm>>>::tl_gbfi; + std::vector< + Eigen::Triplet< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>>::tl_gbfi; /** - * @brief Multithreaded assembly of the Math::SparseMatrix associated to a + * @brief Multithreaded assembly of the Math::SparseMatrix associated to a * BilinearFormBase object. */ template class Multithreaded< - Math::SparseMatrix, - Variational::BilinearForm>> final + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm< + TrialFES, TestFES, + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>> final : public AssemblyBase< - Math::SparseMatrix, - Variational::BilinearForm>> + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm< + TrialFES, TestFES, + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>> { public: - using ScalarType = Real; + using ScalarType = + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>()); using OperatorType = Math::SparseMatrix; @@ -370,14 +444,32 @@ namespace Rodin::Assembly template class Multithreaded< - Math::Matrix, - Variational::BilinearForm>> + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm< + TrialFES, TestFES, + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>> final : public AssemblyBase< - Math::Matrix, - Variational::BilinearForm>> + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>> { public: - using ScalarType = Real; + using ScalarType = + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>()); using OperatorType = Math::Matrix; @@ -564,20 +656,57 @@ namespace Rodin::Assembly mutable std::variant> m_pool; }; - template - thread_local - Math::Matrix - Multithreaded, Variational::BilinearForm>>::tl_res; + template + thread_local + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())> + Multithreaded< + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>::tl_res; - template - thread_local - std::unique_ptr> - Multithreaded, Variational::BilinearForm>>::tl_lbfi; + template + thread_local + std::unique_ptr< + Variational::LocalBilinearFormIntegratorBase< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>> + Multithreaded< + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>::tl_lbfi; - template - thread_local - std::unique_ptr> - Multithreaded, Variational::BilinearForm>>::tl_gbfi; + template + thread_local + std::unique_ptr< + Variational::GlobalBilinearFormIntegratorBase< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>> + Multithreaded< + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>::tl_gbfi; /** * @brief %Multithreaded assembly of the Math::Vector associated to a diff --git a/src/Rodin/Assembly/Sequential.h b/src/Rodin/Assembly/Sequential.h index 836c1502d..f39ee07da 100644 --- a/src/Rodin/Assembly/Sequential.h +++ b/src/Rodin/Assembly/Sequential.h @@ -34,25 +34,118 @@ namespace Rodin::Assembly::Internal namespace Rodin::Assembly { - template + /** + * @brief %Sequential assembly of the Math::Vector associated to a LinearFormBase + * object. + */ + template class Sequential< - std::vector>, - Variational::BilinearForm>>> final + Math::Vector::ScalarType>, + Variational::LinearForm::ScalarType>>> final : public AssemblyBase< - std::vector>, - Variational::BilinearForm>>> + Math::Vector::ScalarType>, + Variational::LinearForm::ScalarType>>> { public: - using ScalarType = Real; + using FESType = FES; - using OperatorType = std::vector>; + using ScalarType = typename FormLanguage::Traits::ScalarType; - using BilinearFormType = Variational::BilinearForm; + using VectorType = Math::Vector; + + using LinearFormType = Variational::LinearForm; + + using Parent = AssemblyBase; + + using InputType = typename Parent::InputType; + + Sequential() = default; + + Sequential(const Sequential& other) + : Parent(other) + {} + + Sequential(Sequential&& other) + : Parent(std::move(other)) + {} + + /** + * @brief Executes the assembly and returns the vector associated to the + * linear form. + */ + VectorType execute(const InputType& input) const override + { + VectorType res(input.getFES().getSize()); + res.setZero(); + const auto& mesh = input.getFES().getMesh(); + for (auto& lfi : input.getLFIs()) + { + const auto& attrs = lfi.getAttributes(); + Internal::SequentialIteration seq(mesh, lfi.getRegion()); + for (auto it = seq.getIterator(); it; ++it) + { + if (attrs.size() == 0 || attrs.count(it->getAttribute())) + { + lfi.setPolytope(*it); + const size_t d = it.getDimension(); + const size_t i = it->getIndex(); + const auto& dofs = input.getFES().getDOFs(d, i); + for (size_t l = 0; l < static_cast(dofs.size()); l++) + res(dofs(l)) += lfi.integrate(l); + } + } + } + return res; + } + + inline + Sequential* copy() const noexcept override + { + return new Sequential(*this); + } + }; + + /** + * @brief Sequential assembly of the Math::SparseMatrix associated to a + * BilinearFormBase object. + */ + template + class Sequential< + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm< + TrialFES, TestFES, + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>> final + : public AssemblyBase< + Math::Matrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>> + { + public: + using ScalarType = + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>()); + + using OperatorType = Math::Matrix; using LocalBilinearFormIntegratorBaseType = Variational::LocalBilinearFormIntegratorBase; using GlobalBilinearFormIntegratorBaseType = Variational::GlobalBilinearFormIntegratorBase; + using BilinearFormType = Variational::BilinearForm; + using Parent = AssemblyBase; using InputType = typename Parent::InputType; @@ -73,9 +166,9 @@ namespace Rodin::Assembly */ OperatorType execute(const InputType& input) const override { - OperatorType res; + OperatorType res(input.getTestFES().getSize(), input.getTrialFES().getSize()); + res.setZero(); const auto& mesh = input.getTrialFES().getMesh(); - res.reserve(input.getTestFES().getSize() * std::log(input.getTrialFES().getSize())); for (auto& bfi : input.getLocalBFIs()) { const auto& attrs = bfi.getAttributes(); @@ -88,14 +181,8 @@ namespace Rodin::Assembly const auto& rows = input.getTestFES().getDOFs(it.getDimension(), it->getIndex()); const auto& cols = input.getTrialFES().getDOFs(it.getDimension(), it->getIndex()); for (size_t l = 0; l < static_cast(rows.size()); l++) - { for (size_t m = 0; m < static_cast(cols.size()); m++) - { - const ScalarType s = bfi.integrate(m, l); - if (s != ScalarType(0)) - res.emplace_back(rows(l), cols(m), s); - } - } + res(rows(l), cols(m)) += bfi.integrate(m, l); } } } @@ -103,12 +190,12 @@ namespace Rodin::Assembly { const auto& trialAttrs = bfi.getTrialAttributes(); const auto& testAttrs = bfi.getTestAttributes(); + Internal::SequentialIteration trialseq(mesh, bfi.getTrialRegion()); Internal::SequentialIteration testseq(mesh, bfi.getTestRegion()); for (auto teIt = testseq.getIterator(); teIt; ++teIt) { if (testAttrs.size() == 0 || testAttrs.count(teIt->getAttribute())) { - Internal::SequentialIteration trialseq(mesh, bfi.getTrialRegion()); for (auto trIt = trialseq.getIterator(); trIt; ++trIt) { if (trialAttrs.size() == 0 || trialAttrs.count(trIt->getAttribute())) @@ -117,14 +204,8 @@ namespace Rodin::Assembly const auto& rows = input.getTestFES().getDOFs(teIt.getDimension(), teIt->getIndex()); const auto& cols = input.getTrialFES().getDOFs(trIt.getDimension(), trIt->getIndex()); for (size_t l = 0; l < static_cast(rows.size()); l++) - { for (size_t m = 0; m < static_cast(cols.size()); m++) - { - const ScalarType s = bfi.integrate(m, l); - if (s != ScalarType(0)) - res.emplace_back(rows(l), cols(m), s); - } - } + res(rows(l), cols(m)) += bfi.integrate(m, l); } } } @@ -141,19 +222,38 @@ namespace Rodin::Assembly }; /** - * @brief Sequential assembly of the Math::SparseMatrix associated to a + * @brief Sequential assembly of the Math::SparseMatrix associated to a * BilinearFormBase object. */ template class Sequential< - Math::SparseMatrix, - Variational::BilinearForm>> final + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm< + TrialFES, TestFES, + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>> final : public AssemblyBase< - Math::SparseMatrix, - Variational::BilinearForm>> + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>, + Variational::BilinearForm< + TrialFES, TestFES, + Math::SparseMatrix< + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>())>>> { public: - using ScalarType = Real; + using ScalarType = + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>()); using OperatorType = Math::SparseMatrix; @@ -180,8 +280,8 @@ namespace Rodin::Assembly OperatorType execute(const InputType& input) const override { Sequential< - std::vector>, - Variational::BilinearForm>>> assembly; + std::vector>, + Variational::BilinearForm>>> assembly; const auto triplets = assembly.execute({ input.getTrialFES(), input.getTestFES(), @@ -198,29 +298,42 @@ namespace Rodin::Assembly } }; - /** - * @brief Sequential assembly of the Math::SparseMatrix associated to a - * BilinearFormBase object. - */ template class Sequential< - Math::Matrix, - Variational::BilinearForm>> final + std::vector::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>> final : public AssemblyBase< - Math::Matrix, - Variational::BilinearForm>> + std::vector::ScalarType>() * + std::declval::ScalarType>())>>, + Variational::BilinearForm::ScalarType>() * + std::declval::ScalarType>())>>>> { public: - using ScalarType = Real; + using ScalarType = + decltype( + std::declval::ScalarType>() * + std::declval::ScalarType>()); - using OperatorType = Math::Matrix; + using OperatorType = std::vector>; + + using BilinearFormType = Variational::BilinearForm; using LocalBilinearFormIntegratorBaseType = Variational::LocalBilinearFormIntegratorBase; using GlobalBilinearFormIntegratorBaseType = Variational::GlobalBilinearFormIntegratorBase; - using BilinearFormType = Variational::BilinearForm; - using Parent = AssemblyBase; using InputType = typename Parent::InputType; @@ -241,9 +354,9 @@ namespace Rodin::Assembly */ OperatorType execute(const InputType& input) const override { - OperatorType res(input.getTestFES().getSize(), input.getTrialFES().getSize()); - res.setZero(); + OperatorType res; const auto& mesh = input.getTrialFES().getMesh(); + res.reserve(input.getTestFES().getSize() * std::log(input.getTrialFES().getSize())); for (auto& bfi : input.getLocalBFIs()) { const auto& attrs = bfi.getAttributes(); @@ -256,8 +369,14 @@ namespace Rodin::Assembly const auto& rows = input.getTestFES().getDOFs(it.getDimension(), it->getIndex()); const auto& cols = input.getTrialFES().getDOFs(it.getDimension(), it->getIndex()); for (size_t l = 0; l < static_cast(rows.size()); l++) + { for (size_t m = 0; m < static_cast(cols.size()); m++) - res(rows(l), cols(m)) += bfi.integrate(m, l); + { + const ScalarType s = bfi.integrate(m, l); + if (s != ScalarType(0)) + res.emplace_back(rows(l), cols(m), s); + } + } } } } @@ -265,12 +384,12 @@ namespace Rodin::Assembly { const auto& trialAttrs = bfi.getTrialAttributes(); const auto& testAttrs = bfi.getTestAttributes(); - Internal::SequentialIteration trialseq(mesh, bfi.getTrialRegion()); Internal::SequentialIteration testseq(mesh, bfi.getTestRegion()); for (auto teIt = testseq.getIterator(); teIt; ++teIt) { if (testAttrs.size() == 0 || testAttrs.count(teIt->getAttribute())) { + Internal::SequentialIteration trialseq(mesh, bfi.getTrialRegion()); for (auto trIt = trialseq.getIterator(); trIt; ++trIt) { if (trialAttrs.size() == 0 || trialAttrs.count(trIt->getAttribute())) @@ -279,8 +398,14 @@ namespace Rodin::Assembly const auto& rows = input.getTestFES().getDOFs(teIt.getDimension(), teIt->getIndex()); const auto& cols = input.getTrialFES().getDOFs(trIt.getDimension(), trIt->getIndex()); for (size_t l = 0; l < static_cast(rows.size()); l++) + { for (size_t m = 0; m < static_cast(cols.size()); m++) - res(rows(l), cols(m)) += bfi.integrate(m, l); + { + const ScalarType s = bfi.integrate(m, l); + if (s != ScalarType(0)) + res.emplace_back(rows(l), cols(m), s); + } + } } } } @@ -479,73 +604,6 @@ namespace Rodin::Assembly return new Sequential(*this); } }; - - /** - * @brief %Sequential assembly of the Math::Vector associated to a LinearFormBase - * object. - */ - template - class Sequential, Variational::LinearForm>> final - : public AssemblyBase, Variational::LinearForm>> - { - public: - using FESType = FES; - - using ScalarType = typename FormLanguage::Traits::ScalarType; - - using VectorType = Math::Vector; - - using LinearFormType = Variational::LinearForm; - - using Parent = AssemblyBase; - - using InputType = typename Parent::InputType; - - Sequential() = default; - - Sequential(const Sequential& other) - : Parent(other) - {} - - Sequential(Sequential&& other) - : Parent(std::move(other)) - {} - - /** - * @brief Executes the assembly and returns the vector associated to the - * linear form. - */ - VectorType execute(const InputType& input) const override - { - VectorType res(input.getFES().getSize()); - res.setZero(); - const auto& mesh = input.getFES().getMesh(); - for (auto& lfi : input.getLFIs()) - { - const auto& attrs = lfi.getAttributes(); - Internal::SequentialIteration seq(mesh, lfi.getRegion()); - for (auto it = seq.getIterator(); it; ++it) - { - if (attrs.size() == 0 || attrs.count(it->getAttribute())) - { - lfi.setPolytope(*it); - const size_t d = it.getDimension(); - const size_t i = it->getIndex(); - const auto& dofs = input.getFES().getDOFs(d, i); - for (size_t l = 0; l < static_cast(dofs.size()); l++) - res(dofs(l)) += lfi.integrate(l); - } - } - } - return res; - } - - inline - Sequential* copy() const noexcept override - { - return new Sequential(*this); - } - }; } #endif diff --git a/src/Rodin/FormLanguage/Traits.h b/src/Rodin/FormLanguage/Traits.h index aed16d73b..2cd305bbb 100644 --- a/src/Rodin/FormLanguage/Traits.h +++ b/src/Rodin/FormLanguage/Traits.h @@ -45,6 +45,36 @@ namespace Rodin::FormLanguage { using ScalarType = Complex; }; + + template + struct Sum + { + using Type = decltype(std::declval() + std::declval()); + }; + + template + struct Minus + { + using Type = decltype(std::declval() - std::declval()); + }; + + template + struct UnaryMinus + { + using Type = decltype(-std::declval()); + }; + + template + struct Division + { + using Type = decltype(std::declval() / std::declval()); + }; + + template + struct Mult + { + using Type = decltype(std::declval() * std::declval()); + }; } #endif diff --git a/src/Rodin/Math/Common.h b/src/Rodin/Math/Common.h index d9f94696a..d80cc2f60 100644 --- a/src/Rodin/Math/Common.h +++ b/src/Rodin/Math/Common.h @@ -99,6 +99,15 @@ namespace Rodin::Math return std::cos(x); } + template + constexpr + inline + std::enable_if_t()))>, bool> + sin(const T& x) + { + return std::sin(x); + } + template constexpr inline diff --git a/src/Rodin/Math/Kernels.h b/src/Rodin/Math/Kernels.h index 6e08a0d29..48a1bdb40 100644 --- a/src/Rodin/Math/Kernels.h +++ b/src/Rodin/Math/Kernels.h @@ -12,20 +12,20 @@ namespace Rodin::Math::Kernels { - inline + template static void eliminate( - SparseMatrix& stiffness, Vector& mass, - const IndexMap& dofs, size_t offset = 0) + SparseMatrix& stiffness, Vector& mass, + const IndexMap& dofs, size_t offset = 0) { - Real* const valuePtr = stiffness.valuePtr(); - SparseMatrix::StorageIndex* const outerPtr = stiffness.outerIndexPtr(); - SparseMatrix::StorageIndex* const innerPtr = stiffness.innerIndexPtr(); + auto* const valuePtr = stiffness.valuePtr(); + auto* const outerPtr = stiffness.outerIndexPtr(); + auto* const innerPtr = stiffness.innerIndexPtr(); // Move essential degrees of freedom in the LHS to the RHS for (const auto& kv : dofs) { const Index& global = kv.first + offset; const auto& dof = kv.second; - for (SparseMatrix::InnerIterator it(stiffness, global); it; ++it) + for (typename SparseMatrix::InnerIterator it(stiffness, global); it; ++it) mass.coeffRef(it.row()) -= it.value() * dof; } for (const auto& [global, dof] : dofs) @@ -39,7 +39,7 @@ namespace Rodin::Math::Kernels assert(innerPtr[i] >= 0); // Assumes CCS format const Index row = innerPtr[i]; - valuePtr[i] = Real(row == global + offset); + valuePtr[i] = (row == global + offset); if (row != global + offset) { for (auto k = outerPtr[row]; 1; k++) @@ -55,9 +55,11 @@ namespace Rodin::Math::Kernels } } - inline - static void replace(const Vector& row, SparseMatrix& stiffness, Vector& mass, - const IndexMap& dofs, size_t offset = 0) + template + static void replace( + const Vector& row, + SparseMatrix& stiffness, Vector& mass, + const IndexMap& dofs, size_t offset = 0) { for (const auto& kv : dofs) { diff --git a/src/Rodin/Variational.h b/src/Rodin/Variational.h index abe0336bf..9d543bcb6 100644 --- a/src/Rodin/Variational.h +++ b/src/Rodin/Variational.h @@ -62,6 +62,9 @@ #include "Variational/Max.h" #include "Variational/Min.h" +#include "Variational/Re.h" +#include "Variational/Im.h" + #include "Variational/Sine.h" #include "Variational/Cosine.h" #include "Variational/Tangent.h" diff --git a/src/Rodin/Variational/BilinearForm.h b/src/Rodin/Variational/BilinearForm.h index d4664fe66..4a60593aa 100644 --- a/src/Rodin/Variational/BilinearForm.h +++ b/src/Rodin/Variational/BilinearForm.h @@ -23,6 +23,15 @@ #include "TestFunction.h" #include "BilinearFormIntegrator.h" +namespace Rodin::FormLanguage +{ + template + struct Traits> + { + using ScalarType = typename FormLanguage::Traits::ScalarType; + }; +} + namespace Rodin::Variational { /** @@ -31,10 +40,12 @@ namespace Rodin::Variational * @see BilinearForm */ - template + template class BilinearFormBase : public FormLanguage::Base { public: + using OperatorType = Operator; + BilinearFormBase() {} @@ -100,8 +111,7 @@ namespace Rodin::Variational using TestFESScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = decltype( - std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; /// Type of operator associated to the bilinear form using OperatorType = Operator; @@ -112,11 +122,11 @@ namespace Rodin::Variational using GlobalBilinearFormIntegratorBaseType = GlobalBilinearFormIntegratorBase; - using LocalBilinearFormIntegratorBaseListType - = FormLanguage::List; + using LocalBilinearFormIntegratorBaseListType = + FormLanguage::List; - using GlobalBilinearFormIntegratorBaseListType - = FormLanguage::List; + using GlobalBilinearFormIntegratorBaseListType = + FormLanguage::List; /// Parent class using Parent = BilinearFormBase; @@ -408,10 +418,12 @@ namespace Rodin::Variational template BilinearForm(TrialFunction&, TestFunction&) - -> BilinearForm::ScalarType>() * - std::declval::ScalarType>())>>; + -> BilinearForm::ScalarType, + typename FormLanguage::Traits::ScalarType> + ::Type>>; } #include "BilinearForm.hpp" diff --git a/src/Rodin/Variational/BooleanFunction.h b/src/Rodin/Variational/BooleanFunction.h index 49608d6fa..dc307ed31 100644 --- a/src/Rodin/Variational/BooleanFunction.h +++ b/src/Rodin/Variational/BooleanFunction.h @@ -73,14 +73,6 @@ namespace Rodin::Variational return static_cast(*this).getValue(p); } - inline - constexpr - void getValue(Math::Vector&, const Geometry::Point&) const = delete; - - inline - constexpr - void getValue(Math::Matrix&, const Geometry::Point&) const = delete; - inline constexpr RangeShape getRangeShape() const diff --git a/src/Rodin/Variational/ComplexFunction.h b/src/Rodin/Variational/ComplexFunction.h index e6e036349..2d4c479c9 100644 --- a/src/Rodin/Variational/ComplexFunction.h +++ b/src/Rodin/Variational/ComplexFunction.h @@ -25,9 +25,9 @@ namespace Rodin::Variational { /** - * @defgroup RealFunctionSpecializations RealFunction Template Specializations - * @brief Template specializations of the RealFunction class. - * @see RealFunction + * @defgroup ComplexFunctionSpecializations ComplexFunction Template Specializations + * @brief Template specializations of the ComplexFunction class. + * @see ComplexFunction */ template @@ -68,47 +68,107 @@ namespace Rodin::Variational virtual ComplexFunctionBase* copy() const noexcept override = 0; }; - /** - * @ingroup ComplexFunctionSpecializations - */ - template - class ComplexFunction> final - : public ComplexFunctionBase> + template <> + class ComplexFunction final + : public ComplexFunctionBase> { public: - using Parent = ComplexFunctionBase>; + using Parent = ComplexFunctionBase>; - using NestedRangeType = typename FormLanguage::Traits>::RangeType; + /** + * @brief Constructs a ComplexFunction from an integer value. + * @param[in] x Constant integer value + */ + ComplexFunction(const Integer& x) + : m_x(x) + {} - static_assert(std::is_same_v); + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_x(other.m_x) + {} - ComplexFunction(const FunctionBase& nested) - : m_nested(nested.copy()) + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_x(other.m_x) + {} + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute) + { + return *this; + } + + inline + constexpr + const Integer& getValue() const + { + return m_x; + } + + inline + constexpr + Complex getValue(const Geometry::Point&) const + { + return static_cast(m_x); + } + + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + const Integer m_x; + }; + + ComplexFunction(Integer) -> ComplexFunction; + + template <> + class ComplexFunction final + : public ComplexFunctionBase> + { + public: + using Parent = ComplexFunctionBase>; + + /** + * @brief Constructs a ComplexFunction from an Real value. + * @param[in] x Constant Real value + */ + ComplexFunction(const Real& x) + : m_x(x) {} ComplexFunction(const ComplexFunction& other) : Parent(other), - m_nested(other.m_nested->copy()) + m_x(other.m_x) {} ComplexFunction(ComplexFunction&& other) : Parent(std::move(other)), - m_nested(std::move(other.m_nested)) + m_x(other.m_x) {} inline constexpr - auto getValue(const Geometry::Point& v) const + ComplexFunction& traceOf(Geometry::Attribute) { - return m_nested->getValue(v); + return *this; } inline constexpr - ComplexFunction& traceOf(Geometry::Attribute attrs) + const Real& getValue() const { - m_nested->traceOf(attrs); - return *this; + return m_x; + } + + inline + constexpr + Complex getValue(const Geometry::Point&) const + { + return static_cast(m_x); } inline ComplexFunction* copy() const noexcept override @@ -117,14 +177,10 @@ namespace Rodin::Variational } private: - std::unique_ptr> m_nested; + const Real m_x; }; - /** - * @brief CTAD for ComplexFunction. - */ - template - ComplexFunction(const FunctionBase&) -> ComplexFunction>; + ComplexFunction(Real) -> ComplexFunction; /** * @ingroup ComplexFunctionSpecializations @@ -141,7 +197,7 @@ namespace Rodin::Variational * @brief Constructs a ComplexFunction from a constant scalar value. * @param[in] x Constant scalar value */ - ComplexFunction(Complex x) + ComplexFunction(const Complex& x) : m_x(x) {} @@ -164,7 +220,7 @@ namespace Rodin::Variational inline constexpr - Complex getValue() const + const Complex& getValue() const { return m_x; } @@ -188,52 +244,115 @@ namespace Rodin::Variational /** * @brief CTAD for ComplexFunction. */ - ComplexFunction(Complex) -> ComplexFunction; + ComplexFunction(const Complex&) -> ComplexFunction; - template <> - class ComplexFunction final - : public ComplexFunctionBase> + /** + * @ingroup ComplexFunctionSpecializations + */ + template + class ComplexFunction> final + : public ComplexFunctionBase> { public: - using Parent = ComplexFunctionBase>; + using Parent = ComplexFunctionBase>; - /** - * @brief Constructs a ComplexFunction from an integer value. - * @param[in] x Constant integer value - */ - ComplexFunction(Integer x) - : m_x(x) + using ScalarType = Complex; + + ComplexFunction(const FunctionBase& nested) + : m_nested(nested.copy()) {} ComplexFunction(const ComplexFunction& other) : Parent(other), - m_x(other.m_x) + m_nested(other.m_nested->copy()) {} ComplexFunction(ComplexFunction&& other) : Parent(std::move(other)), - m_x(other.m_x) + m_nested(std::move(other.m_nested)) {} inline constexpr - ComplexFunction& traceOf(Geometry::Attribute) + ScalarType getValue(const Geometry::Point& v) const { + return m_nested->getValue(v); + } + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute attrs) + { + m_nested->traceOf(attrs); return *this; } + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + std::unique_ptr> m_nested; + }; + + /** + * @brief CTAD for ComplexFunction. + */ + template + ComplexFunction(const FunctionBase&) -> ComplexFunction>; + + /** + * @ingroup ComplexFunctionSpecializations + */ + template + class ComplexFunction, FunctionBase> final + : public ComplexFunctionBase, FunctionBase>> + { + public: + using RealFunctionType = FunctionBase; + + using ImagFunctionType = FunctionBase; + + using Parent = + ComplexFunctionBase>; + + using RealFunctionRangeType = typename FormLanguage::Traits::RangeType; + + using ImagFunctionRangeType = typename FormLanguage::Traits::RangeType; + + static_assert(std::is_same_v); + + static_assert(std::is_same_v); + + ComplexFunction(const RealFunctionType& re, const ImagFunctionType& imag) + : m_re(re.copy()), m_imag(imag.copy()) + {} + + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_re(m_re->copy()), m_imag(m_imag->copy()) + {} + + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_re(std::move(m_re)), m_imag(std::move(m_imag)) + {} + inline constexpr - Complex getValue() const + Complex getValue(const Geometry::Point& p) const { - return m_x; + return { m_re->getValue(p), m_imag->getValue(p) }; } inline constexpr - Complex getValue(const Geometry::Point&) const + ComplexFunction& traceOf(Geometry::Attribute attrs) { - return static_cast(m_x); + m_re->traceOf(attrs); + m_imag->traceOf(attrs); + return *this; } inline ComplexFunction* copy() const noexcept override @@ -242,10 +361,16 @@ namespace Rodin::Variational } private: - const Integer m_x; + std::unique_ptr m_re; + std::unique_ptr m_imag; }; - ComplexFunction(Integer) -> ComplexFunction; + /** + * @brief CTAD for ComplexFunction. + */ + template + ComplexFunction(const FunctionBase&, const FunctionBase&) + -> ComplexFunction, FunctionBase>; /** * @ingroup ComplexFunctionSpecializations @@ -259,7 +384,7 @@ namespace Rodin::Variational public: using Parent = ComplexFunctionBase>; - ComplexFunction(F f) + ComplexFunction(const F& f) : m_f(f) {} @@ -301,7 +426,74 @@ namespace Rodin::Variational */ template >> - ComplexFunction(F) -> ComplexFunction; + ComplexFunction(const F&) -> ComplexFunction; + + template + class ComplexFunction final + : public ComplexFunctionBase> + { + static_assert(std::is_invocable_v); + static_assert(std::is_invocable_v); + + public: + using Parent = ComplexFunctionBase>; + + ComplexFunction(const FReal& re, const FImag& imag) + : m_re(re), m_imag(imag) + {} + + ComplexFunction(const ComplexFunction& other) + : Parent(other), + m_re(other.m_re), m_imag(other.m_imag) + {} + + ComplexFunction(ComplexFunction&& other) + : Parent(std::move(other)), + m_re(std::move(other.m_re)), m_imag(std::move(other.m_imag)) + {} + + inline + constexpr + ComplexFunction& traceOf(Geometry::Attribute attr) + { + m_re->traceOf(attr); + m_imag->traceOf(attr); + return *this; + } + + inline + constexpr + ComplexFunction& traceOf(const FlatSet& attrs) + { + m_re->traceOf(attrs); + m_imag->traceOf(attrs); + return *this; + } + + inline + constexpr + Complex getValue(const Geometry::Point& p) const + { + return { m_re->getValue(p), m_imag->getValue(p) }; + } + + inline ComplexFunction* copy() const noexcept override + { + return new ComplexFunction(*this); + } + + private: + const FReal m_re; + const FImag m_imag; + }; + + /** + * @brief CTAD for ComplexFunction. + */ + template && std::is_invocable_v>> + ComplexFunction(const FReal&, const FImag&) -> ComplexFunction; } #endif diff --git a/src/Rodin/Variational/Cosine.h b/src/Rodin/Variational/Cosine.h index d659af3a8..239a033d6 100644 --- a/src/Rodin/Variational/Cosine.h +++ b/src/Rodin/Variational/Cosine.h @@ -30,41 +30,51 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; + using ScalarType = Real; + using Parent = RealFunctionBase>>; Cos(const OperandType& v) - : m_v(v.copy()) + : m_operand(v.copy()) {} Cos(const Cos& other) : Parent(other), - m_v(other.m_v->copy()) + m_operand(other.m_operand->copy()) {} Cos(Cos&& other) : Parent(std::move(other)), - m_v(std::move(other.m_v)) + m_operand(std::move(other.m_operand)) {} inline constexpr - Cos& traceOf(Geometry::Attribute attrs) + Cos& traceOf(Geometry::Attribute attr) + { + m_operand->traceOf(attr); + return *this; + } + + inline + constexpr + Cos& traceOf(const FlatSet& attrs) { - m_v.traceOf(attrs); + m_operand->traceOf(attrs); return *this; } inline - auto getValue(const Geometry::Point& p) const + ScalarType getValue(const Geometry::Point& p) const { - return std::cos(static_cast(getOperand().getValue(p))); + return Math::cos(getOperand().getValue(p)); } inline const OperandType& getOperand() const { - assert(m_v); - return *m_v; + assert(m_operand); + return *m_operand; } inline Cos* copy() const noexcept override @@ -73,7 +83,7 @@ namespace Rodin::Variational } private: - std::unique_ptr m_v; + std::unique_ptr m_operand; }; template diff --git a/src/Rodin/Variational/DenseProblem.h b/src/Rodin/Variational/DenseProblem.h index 6e75100d5..7738a5693 100644 --- a/src/Rodin/Variational/DenseProblem.h +++ b/src/Rodin/Variational/DenseProblem.h @@ -31,7 +31,12 @@ namespace Rodin::Variational { template DenseProblem(TrialFunction&, TestFunction&) - -> DenseProblem, Math::Vector>; + -> DenseProblem::ScalarType, + typename FormLanguage::Traits::ScalarType>::Type>, + Math::Vector::ScalarType>>; /** * @defgroup DenseProblemSpecializations DenseProblem Template Specializations @@ -41,17 +46,26 @@ namespace Rodin::Variational /** * @ingroup DenseProblemSpecializations - * @brief General class to assemble linear systems with `Math::Matrix` - * and `Math::Vector` types in a serial context. + * @brief General class to assemble linear systems with `Math::Matrix` + * and `Math::Vector` types in a serial context. */ template class DenseProblem< TrialFES, TestFES, - Math::Matrix::ScalarType>, + Math::Matrix< + typename FormLanguage::Mult< + typename FormLanguage::Traits::ScalarType, + typename FormLanguage::Traits::ScalarType>::Type>, Math::Vector::ScalarType>> : public ProblemBase< - Math::Matrix::ScalarType>, - Math::Vector::ScalarType>> + Math::Matrix< + typename FormLanguage::Mult< + typename FormLanguage::Traits::ScalarType, + typename FormLanguage::Traits::ScalarType>::Type>, + Math::Vector::ScalarType>, + typename FormLanguage::Mult< + typename FormLanguage::Traits::ScalarType, + typename FormLanguage::Traits::ScalarType>::Type> { public: using TrialFESScalarType = @@ -61,19 +75,21 @@ namespace Rodin::Variational typename FormLanguage::Traits::ScalarType; using OperatorScalarType = - decltype(std::declval() * std::declval()); + typename FormLanguage::Mult::Type; using VectorScalarType = TestFESScalarType; + using ScalarType = OperatorScalarType; + using ContextType = Context::Sequential; - using OperatorType = Math::Matrix; + using OperatorType = Math::Matrix; - using VectorType = Math::Vector; + using VectorType = Math::Vector; using LinearFormIntegratorBaseType = LinearFormIntegratorBase; - using Parent = ProblemBase, Math::Vector>; + using Parent = ProblemBase; /** * @brief Constructs an empty DenseProblem involving the trial function @f$ u @f$ @@ -191,7 +207,7 @@ namespace Rodin::Variational // Eliminate the parent column, adding it to the child columns for (size_t i = 0; i < count; i++) { - const Real coeff = coeffs.coeff(i); + const ScalarType coeff = coeffs.coeff(i); const Index child = children.coeff(i); m_stiffness.col(child) += coeff * m_stiffness.col(parent); } @@ -199,13 +215,13 @@ namespace Rodin::Variational m_stiffness.col(parent).setZero(); // Eliminate the parent row, adding it to the child rows - IndexMap parentLookup; - std::vector> childrenLookup(children.size()); + IndexMap parentLookup; + std::vector> childrenLookup(children.size()); for (size_t col = 0; col < static_cast(m_stiffness.cols()); col++) { bool parentFound = false; size_t childrenFound = 0; - for (Math::Matrix::InnerIterator it(m_stiffness, col); it; ++it) + for (typename OperatorType::InnerIterator it(m_stiffness, col); it; ++it) { if (parentFound && childrenFound == count) { @@ -240,7 +256,7 @@ namespace Rodin::Variational { for (size_t i = 0; i < count; i++) { - const Real coeff = coeffs.coeff(i); + const ScalarType coeff = coeffs.coeff(i); childrenLookup[i][col] += coeff * value; } } @@ -255,7 +271,7 @@ namespace Rodin::Variational // Eliminate the parent entry, adding it to the child entries for (size_t i = 0; i < count; i++) { - const Real coeff = coeffs.coeff(i); + const ScalarType coeff = coeffs.coeff(i); const Index child = children.coeff(i); m_mass.coeffRef(child) += coeff * m_mass.coeff(parent); } @@ -269,7 +285,7 @@ namespace Rodin::Variational assert(children.size() >= 0); for (size_t i = 0; i < static_cast(children.size()); i++) { - const Real coeff = coeffs.coeff(i); + const ScalarType coeff = coeffs.coeff(i); const Index child = children.coeff(i); m_stiffness.coeffRef(parent, child) = -coeff; } @@ -359,7 +375,7 @@ namespace Rodin::Variational getTrialFunction().getSolution().setWeights(std::move(m_guess)); } - DenseProblem& operator=(const ProblemBody& rhs) override + DenseProblem& operator=(const ProblemBody& rhs) override { for (auto& bfi : rhs.getLocalBFIs()) m_bilinearForm.add(bfi); diff --git a/src/Rodin/Variational/Div.h b/src/Rodin/Variational/Div.h index ced60d2f6..23731fed8 100644 --- a/src/Rodin/Variational/Div.h +++ b/src/Rodin/Variational/Div.h @@ -34,6 +34,8 @@ namespace Rodin::Variational public: using FESType = FES; + using ScalarType = typename FormLanguage::Traits::ScalarType; + using OperandType = GridFunction; /// Parent class @@ -65,9 +67,9 @@ namespace Rodin::Variational {} inline - Real getValue(const Geometry::Point& p) const + ScalarType getValue(const Geometry::Point& p) const { - Real out; + ScalarType out; const auto& polytope = p.getPolytope(); const auto& polytopeMesh = polytope.getMesh(); const auto& gf = getOperand(); @@ -90,7 +92,6 @@ namespace Rodin::Variational else { assert(false); - out = NAN; } return out; } @@ -107,9 +108,9 @@ namespace Rodin::Variational */ inline constexpr - auto interpolate(Real& out, const Geometry::Point& p) const + void interpolate(Real& out, const Geometry::Point& p) const { - return static_cast(*this).interpolate(out, p); + static_cast(*this).interpolate(out, p); } /** diff --git a/src/Rodin/Variational/ForwardDecls.h b/src/Rodin/Variational/ForwardDecls.h index 8b0b7ca7e..352c6efcc 100644 --- a/src/Rodin/Variational/ForwardDecls.h +++ b/src/Rodin/Variational/ForwardDecls.h @@ -1050,13 +1050,16 @@ namespace Rodin::Variational template class PeriodicBC; - template + template + class ProblemBodyBase; + + template class ProblemBody; /** * @brief Base class for variational problem objects. */ - template + template class ProblemBase; /** diff --git a/src/Rodin/Variational/Grad.h b/src/Rodin/Variational/Grad.h index 75b3c562d..d3218f329 100644 --- a/src/Rodin/Variational/Grad.h +++ b/src/Rodin/Variational/Grad.h @@ -127,9 +127,9 @@ namespace Rodin::Variational */ inline constexpr - auto interpolate(Math::SpatialVector& out, const Geometry::Point& p) const + void interpolate(Math::SpatialVector& out, const Geometry::Point& p) const { - return static_cast(*this).interpolate(out, p); + static_cast(*this).interpolate(out, p); } inline diff --git a/src/Rodin/Variational/GridFunction.h b/src/Rodin/Variational/GridFunction.h index 1701acfc6..564d85812 100644 --- a/src/Rodin/Variational/GridFunction.h +++ b/src/Rodin/Variational/GridFunction.h @@ -524,13 +524,13 @@ namespace Rodin::Variational template Derived& project(const FunctionBase& fn, const FlatSet& attrs) { - using FunctionType = FunctionBase; - using FunctionRangeType = typename FormLanguage::Traits::RangeType; - static_assert(std::is_same_v); + // using FunctionType = FunctionBase; + // using FunctionRangeType = typename FormLanguage::Traits::RangeType; + // static_assert(std::is_same_v); const auto& fes = getFiniteElementSpace(); const auto& mesh = fes.getMesh(); const size_t d = mesh.getDimension(); - std::vector ns(fes.getSize(), 0); + std::vector ns(fes.getSize(), 0); if constexpr (std::is_same_v) { #ifdef RODIN_MULTITHREADED @@ -568,7 +568,7 @@ namespace Rodin::Variational { const Index global = is[i]; m_data(global) = - (vs[i] + ns[global] * m_data(global)) / (ns[global] + 1.0); + (vs[i] + ns[global] * m_data(global)) / (ns[global] + 1); ns[global] += 1; } m_mutex.unlock(); @@ -590,7 +590,7 @@ namespace Rodin::Variational const Index global = fes.getGlobalIndex({ d, i }, local); assert(m_data.rows() == 1); m_data(global) = - (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1.0); + (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1); ns[global] += 1; } } @@ -614,7 +614,7 @@ namespace Rodin::Variational const Index global = fes.getGlobalIndex({ d, i }, local); fn.getDerived().getValue(value, p); m_data.col(global) = - (value + ns[global] * m_data.col(global)) / (ns[global] + 1.0); + (value + ns[global] * m_data.col(global)) / (ns[global] + 1); ns[global] += 1; } } @@ -700,13 +700,13 @@ namespace Rodin::Variational template Derived& projectOnBoundary(const FunctionBase& fn, const FlatSet& attrs) { - using FunctionType = FunctionBase; - using FunctionRangeType = typename FormLanguage::Traits::RangeType; - static_assert(std::is_same_v); + // using FunctionType = FunctionBase; + // using FunctionRangeType = typename FormLanguage::Traits::RangeType; + // static_assert(std::is_same_v); const auto& fes = getFiniteElementSpace(); const auto& mesh = fes.getMesh(); const size_t d = mesh.getDimension() - 1; - std::vector ns(fes.getSize(), 0); + std::vector ns(fes.getSize(), 0); for (auto it = mesh.getBoundary(); !it.end(); ++it) { const auto& polytope = *it; @@ -723,12 +723,12 @@ namespace Rodin::Variational { assert(m_data.rows() == 1); m_data(global) = - (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1.0); + (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1); } else if constexpr (std::is_same_v>) { m_data.col(global) = - (fn.getValue(p) + ns[global] * m_data.col(global)) / (ns[global] + 1.0); + (fn.getValue(p) + ns[global] * m_data.col(global)) / (ns[global] + 1); } else { @@ -812,13 +812,13 @@ namespace Rodin::Variational template Derived& projectOnFaces(const FunctionBase& fn, const FlatSet& attrs) { - using FunctionType = FunctionBase; - using FunctionRangeType = typename FormLanguage::Traits::RangeType; - static_assert(std::is_same_v); + // using FunctionType = FunctionBase; + // using FunctionRangeType = typename FormLanguage::Traits::RangeType; + // static_assert(std::is_same_v); const auto& fes = getFiniteElementSpace(); const auto& mesh = fes.getMesh(); const size_t d = mesh.getDimension() - 1; - std::vector ns(fes.getSize(), 0); + std::vector ns(fes.getSize(), 0); for (auto it = mesh.getFace(); !it.end(); ++it) { const auto& polytope = *it; @@ -835,12 +835,12 @@ namespace Rodin::Variational { assert(m_data.rows() == 1); m_data(global) = - (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1.0); + (fn.getValue(p) + ns[global] * m_data(global)) / (ns[global] + 1); } else if constexpr (std::is_same_v>) { m_data.col(global) = - (fn.getValue(p) + ns[global] * m_data.col(global)) / (ns[global] + 1.0); + (fn.getValue(p) + ns[global] * m_data.col(global)) / (ns[global] + 1); } else { @@ -925,9 +925,9 @@ namespace Rodin::Variational Derived& projectOnInterfaces( const FunctionBase& fn, const FlatSet& attrs) { - using FunctionType = FunctionBase; - using FunctionRangeType = typename FormLanguage::Traits::RangeType; - static_assert(std::is_same_v); + // using FunctionType = FunctionBase; + // using FunctionRangeType = typename FormLanguage::Traits::RangeType; + // static_assert(std::is_same_v); const auto& fes = getFiniteElementSpace(); const auto& mesh = fes.getMesh(); const size_t d = mesh.getDimension() - 1; @@ -1241,19 +1241,6 @@ namespace Rodin::Variational } } - /** - * @brief Interpolates the GridFunction at the given point. - * @note CRTP function to be overriden in Derived class. - */ - inline - constexpr - RangeType interpolate(const Geometry::Point& p) const - { - RangeType res; - res = static_cast(*this).interpolate(res, p); - return res; - } - /** * @brief Interpolates the GridFunction at the given point. * @note CRTP function to be overriden in Derived class. diff --git a/src/Rodin/Variational/Im.h b/src/Rodin/Variational/Im.h new file mode 100644 index 000000000..a5b02a962 --- /dev/null +++ b/src/Rodin/Variational/Im.h @@ -0,0 +1,63 @@ +#ifndef RODDIN_VARIATIONAL_IM_H +#define RODDIN_VARIATIONAL_IM_H + +#include "ForwardDecls.h" +#include "RealFunction.h" + +namespace Rodin::Variational +{ + template + class Im; + + template + class Im> : public RealFunctionBase>> + { + public: + using OperandType = FunctionBase; + + using ScalarType = Real; + + using Parent = RealFunctionBase>>; + + Im(const OperandType& f) + : m_operand(f.copy()) + {} + + Im(const Im& other) + : Parent(other), + m_operand(other.m_operand->copy()) + {} + + Im(Im&& other) + : Parent(std::move(other)), + m_operand(std::move(other.m_operand)) + {} + + const OperandType& getOperand() const + { + assert(m_operand); + return *m_operand; + } + + inline + constexpr + ScalarType getValue(const Geometry::Point& p) const + { + return getOperand().getValue(p).imag(); + } + + inline Im* copy() const noexcept override + { + return new Im(*this); + } + + private: + std::unique_ptr m_operand; + }; + + template + Im(const FunctionBase&) -> Im>; + +} + +#endif diff --git a/src/Rodin/Variational/Jacobian.h b/src/Rodin/Variational/Jacobian.h index 930f588b8..9aa5c4dff 100644 --- a/src/Rodin/Variational/Jacobian.h +++ b/src/Rodin/Variational/Jacobian.h @@ -41,6 +41,8 @@ namespace Rodin::Variational using ScalarType = typename FormLanguage::Traits::ScalarType; + using SpatialMatrixType = Math::SpatialMatrix; + using OperandType = GridFunction; using Parent = @@ -86,17 +88,16 @@ namespace Rodin::Variational } inline - Math::SpatialMatrix getValue(const Geometry::Point& p) const + SpatialMatrixType getValue(const Geometry::Point& p) const { - Math::SpatialMatrix out; + SpatialMatrixType out; getValue(out, p); return out; } inline - void getValue(Math::SpatialMatrix& out, const Geometry::Point& p) const + void getValue(SpatialMatrixType& out, const Geometry::Point& p) const { - out.setConstant(NAN); const auto& polytope = p.getPolytope(); const auto& polytopeMesh = polytope.getMesh(); const auto& gf = getOperand(); @@ -134,9 +135,9 @@ namespace Rodin::Variational */ inline constexpr - auto interpolate(Math::SpatialMatrix& out, const Geometry::Point& p) const + void interpolate(SpatialMatrixType& out, const Geometry::Point& p) const { - return static_cast(*this).interpolate(out, p); + static_cast(*this).interpolate(out, p); } /** diff --git a/src/Rodin/Variational/LinearForm.h b/src/Rodin/Variational/LinearForm.h index 7678c7dbb..8925c5996 100644 --- a/src/Rodin/Variational/LinearForm.h +++ b/src/Rodin/Variational/LinearForm.h @@ -23,10 +23,12 @@ namespace Rodin::Variational { - template + template class LinearFormBase : public FormLanguage::Base { public: + using VectorType = Vector; + LinearFormBase() = default; LinearFormBase(const LinearFormBase& other) diff --git a/src/Rodin/Variational/MatrixFunction.h b/src/Rodin/Variational/MatrixFunction.h index 62660a17b..e46c21e71 100644 --- a/src/Rodin/Variational/MatrixFunction.h +++ b/src/Rodin/Variational/MatrixFunction.h @@ -40,8 +40,6 @@ namespace Rodin::Variational public: using ScalarType = Scalar; - using MatrixType = Math::Matrix; - using Parent = FunctionBase>; MatrixFunctionBase() = default; @@ -69,6 +67,7 @@ namespace Rodin::Variational return static_cast(*this).getValue(p); } + template inline constexpr void getValue(MatrixType& res, const Geometry::Point& p) const diff --git a/src/Rodin/Variational/Max.h b/src/Rodin/Variational/Max.h index 99c82cd87..5cffe0d44 100644 --- a/src/Rodin/Variational/Max.h +++ b/src/Rodin/Variational/Max.h @@ -26,14 +26,14 @@ namespace Rodin::Variational */ template class Max, FunctionBase> final - : public RealFunctionBase, FunctionBase>> + : public FunctionBase, FunctionBase>> { public: using LHSType = FunctionBase; using RHSType = FunctionBase; - using Parent = RealFunctionBase, FunctionBase>>; + using Parent = FunctionBase, FunctionBase>>; Max(const LHSType& a, const RHSType& b) : m_lhs(a.copy()), m_rhs(b.copy()) @@ -60,7 +60,7 @@ namespace Rodin::Variational inline constexpr - Real getValue(const Geometry::Point& p) const + auto getValue(const Geometry::Point& p) const { const auto lhs = getLHS().getValue(p); const auto rhs = getRHS().getValue(p); diff --git a/src/Rodin/Variational/Min.h b/src/Rodin/Variational/Min.h index 6347a88b5..7176f7579 100644 --- a/src/Rodin/Variational/Min.h +++ b/src/Rodin/Variational/Min.h @@ -27,14 +27,14 @@ namespace Rodin::Variational */ template class Min, FunctionBase> final - : public RealFunctionBase, FunctionBase>> + : public FunctionBase, FunctionBase>> { public: using LHSType = FunctionBase; using RHSType = FunctionBase; - using Parent = RealFunctionBase, FunctionBase>>; + using Parent = FunctionBase, FunctionBase>>; Min(const LHSType& a, const RHSType& b) : m_lhs(a.copy()), m_rhs(b.copy()) @@ -61,7 +61,7 @@ namespace Rodin::Variational inline constexpr - Real getValue(const Geometry::Point& p) const + auto getValue(const Geometry::Point& p) const { const auto lhs = getLHS().getValue(p); const auto rhs = getRHS().getValue(p); @@ -106,7 +106,9 @@ namespace Rodin::Variational public: using LHSType = FunctionBase; - using RHSType = Real; + using ScalarType = Real; + + using RHSType = ScalarType; using Parent = RealFunctionBase, RHSType>>; @@ -137,7 +139,7 @@ namespace Rodin::Variational inline constexpr - Real getValue(const Geometry::Point& p) const + ScalarType getValue(const Geometry::Point& p) const { const auto lhs = getLHS().getValue(p); const auto& rhs = getRHS(); diff --git a/src/Rodin/Variational/P0/ForwardDecls.h b/src/Rodin/Variational/P0/ForwardDecls.h index b76e72b3b..ae7b4ed89 100644 --- a/src/Rodin/Variational/P0/ForwardDecls.h +++ b/src/Rodin/Variational/P0/ForwardDecls.h @@ -28,7 +28,7 @@ namespace Rodin::Variational /** * @brief Degree 0 Lagrange finite element space - * @tparam Range Range value type. Either Rodin::Real or Math::Vector. + * @tparam Range Range value type * @tparam Context Context type * @tparam Args Additional arguments * @@ -59,6 +59,9 @@ namespace Rodin::Variational */ using RealP0Element = P0Element; + /** + * @brief Alias for P0Element + */ using ComplexP0Element = P0Element; /** diff --git a/src/Rodin/Variational/P0/GridFunction.h b/src/Rodin/Variational/P0/GridFunction.h index becf9989c..3a32f8d5a 100644 --- a/src/Rodin/Variational/P0/GridFunction.h +++ b/src/Rodin/Variational/P0/GridFunction.h @@ -41,6 +41,8 @@ namespace Rodin::Variational public: using FESType = P0; + using ScalarType = typename FormLanguage::Traits::ScalarType; + using MeshType = typename FormLanguage::Traits::MeshType; using RangeType = typename FormLanguage::Traits::RangeType; @@ -96,7 +98,7 @@ namespace Rodin::Variational GridFunction& operator=(const GridFunction&) = delete; - void interpolate(Real& res, const Geometry::Point& p) const + void interpolate(RangeType& res, const Geometry::Point& p) const { static_assert(std::is_same_v); const auto& fes = this->getFiniteElementSpace(); @@ -107,29 +109,20 @@ namespace Rodin::Variational assert(d == polytope.getDimension()); const Index i = polytope.getIndex(); assert(fes.getFiniteElement(d, i).getCount() == 1); - res = getValue({ d, i }, 0); + if constexpr (std::is_same_v) + { + res = getValue({ d, i }, 0); + } + else if constexpr (std::is_same_v>) + { + assert(false); + } + else + { + assert(false); + } } - // void interpolate(Math::Vector& res, const Geometry::Point& p) const - // { - // static_assert(std::is_same_v>); - // const auto& fes = this->getFiniteElementSpace(); - // const auto& polytope = p.getPolytope(); - // const size_t d = polytope.getDimension(); - // const Index i = polytope.getIndex(); - // const auto& fe = fes.getFiniteElement(d, i); - // const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); - // const size_t vdim = fes.getVectorDimension(); - // const size_t dofs = fe.getCount(); - // res.resize(vdim); - // res.setZero(); - // for (Index local = 0; local < dofs; local++) - // { - // const auto& basis = fe.getBasis(local); - // res += getValue({d, i}, local).coeff(local % vdim) * basis(r); - // } - // } - GridFunction& setWeights() { auto& data = this->getData(); diff --git a/src/Rodin/Variational/P1/Div.h b/src/Rodin/Variational/P1/Div.h index 916748362..ed0eb41f7 100644 --- a/src/Rodin/Variational/P1/Div.h +++ b/src/Rodin/Variational/P1/Div.h @@ -48,9 +48,11 @@ namespace Rodin::Variational Div, Mesh>>>> { public: - using ScalarType = Scalar; + using FESType = Variational::P1, Mesh>; - using FESType = Variational::P1, Mesh>; + using ScalarType = typename FormLanguage::Traits::ScalarType; + + using SpatialMatrixType = Math::SpatialMatrix; /// Operand type using OperandType = GridFunction; @@ -97,7 +99,7 @@ namespace Rodin::Variational if (inc.size() == 1) { const auto& tracePolytope = mesh.getPolytope(meshDim, *inc.begin()); - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const auto rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -122,7 +124,7 @@ namespace Rodin::Variational const auto& tracePolytope = mesh.getPolytope(meshDim, idx); if (traceDomain.count(tracePolytope->getAttribute())) { - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const auto rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -142,7 +144,7 @@ namespace Rodin::Variational const auto& vdim = fes.getVectorDimension(); const auto& fe = fes.getFiniteElement(d, i); const auto& rc = p.getReferenceCoordinates(); - Math::SpatialMatrix jacobian(vdim, d); + SpatialMatrixType jacobian(vdim, d); out = 0; for (size_t local = 0; local < fe.getCount(); local++) { diff --git a/src/Rodin/Variational/P1/Grad.h b/src/Rodin/Variational/P1/Grad.h index eae629575..35171475d 100644 --- a/src/Rodin/Variational/P1/Grad.h +++ b/src/Rodin/Variational/P1/Grad.h @@ -61,6 +61,8 @@ namespace Rodin::Variational using ScalarType = typename FormLanguage::Traits::ScalarType; + using SpatialVectorType = Math::SpatialVector; + using OperandType = GridFunction; using Parent = GradBase>; @@ -88,14 +90,7 @@ namespace Rodin::Variational : Parent(std::move(other)) {} - void interpolate(Math::Vector& out, const Geometry::Point& p) const - { - Math::SpatialVector tmp; - interpolate(tmp, p); - out = std::move(tmp); - } - - void interpolate(Math::SpatialVector& out, const Geometry::Point& p) const + void interpolate(SpatialVectorType& out, const Geometry::Point& p) const { const auto& polytope = p.getPolytope(); const auto& d = polytope.getDimension(); @@ -111,7 +106,7 @@ namespace Rodin::Variational if (inc.size() == 1) { const auto& tracePolytope = mesh.getPolytope(meshDim, *inc.begin()); - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const auto rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -136,7 +131,7 @@ namespace Rodin::Variational const auto& tracePolytope = mesh.getPolytope(meshDim, idx); if (traceDomain.count(tracePolytope->getAttribute())) { - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const auto rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -155,8 +150,8 @@ namespace Rodin::Variational const auto& fes = gf.getFiniteElementSpace(); const auto& fe = fes.getFiniteElement(d, i); const auto& rc = p.getReferenceCoordinates(); - Math::SpatialVector grad(d); - Math::SpatialVector res(d); + SpatialVectorType grad(d); + SpatialVectorType res(d); res.setZero(); for (size_t local = 0; local < fe.getCount(); local++) { diff --git a/src/Rodin/Variational/P1/GridFunction.h b/src/Rodin/Variational/P1/GridFunction.h index 28da7556f..b8f138cc7 100644 --- a/src/Rodin/Variational/P1/GridFunction.h +++ b/src/Rodin/Variational/P1/GridFunction.h @@ -42,6 +42,8 @@ namespace Rodin::Variational /// Type of finite element space to which the GridFunction belongs to using FESType = P1; + using ScalarType = typename FormLanguage::Traits::ScalarType; + using MeshType = typename FormLanguage::Traits::MeshType; using RangeType = typename FormLanguage::Traits::RangeType; @@ -98,9 +100,8 @@ namespace Rodin::Variational GridFunction& operator=(const GridFunction&) = delete; - void interpolate(Real& res, const Geometry::Point& p) const + void interpolate(RangeType& res, const Geometry::Point& p) const { - static_assert(std::is_same_v); const auto& fes = this->getFiniteElementSpace(); const auto& fesMesh = fes.getMesh(); const auto& polytope = p.getPolytope(); @@ -109,31 +110,37 @@ namespace Rodin::Variational const Index i = polytope.getIndex(); const auto& fe = fes.getFiniteElement(d, i); const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); - res = 0; - for (Index local = 0; local < fe.getCount(); local++) + if constexpr (std::is_same_v) { - const auto& basis = fe.getBasis(local); - res += getValue({d, i}, local) * basis(r); + res = Real(0); + for (Index local = 0; local < fe.getCount(); local++) + res += getValue({d, i}, local) * fe.getBasis(local)(r); } - } - - void interpolate(Math::Vector& res, const Geometry::Point& p) const - { - static_assert(std::is_same_v>); - const auto& fes = this->getFiniteElementSpace(); - const auto& polytope = p.getPolytope(); - const size_t d = polytope.getDimension(); - const Index i = polytope.getIndex(); - const auto& fe = fes.getFiniteElement(d, i); - const auto& r = p.getCoordinates(Geometry::Point::Coordinates::Reference); - const size_t vdim = fes.getVectorDimension(); - const size_t dofs = fe.getCount(); - res.resize(vdim); - res.setZero(); - for (Index local = 0; local < dofs; local++) + else if constexpr (std::is_same_v) { - const auto& basis = fe.getBasis(local); - res += getValue({d, i}, local).coeff(local % vdim) * basis(r); + res = Complex(0, 0); + for (Index local = 0; local < fe.getCount(); local++) + { + if (local % 2 == 0) + res += getValue({d, i}, local).real() * fe.getBasis(local)(r); + else + res += getValue({d, i}, local).imag() * fe.getBasis(local)(r); + } + } + else if constexpr (std::is_same_v>) + { + const size_t vdim = fes.getVectorDimension(); + res.resize(vdim); + res.setZero(); + for (Index local = 0; local < fe.getCount(); local++) + { + const auto& basis = fe.getBasis(local); + res += getValue({d, i}, local).coeff(local % vdim) * basis(r); + } + } + else + { + assert(false); } } diff --git a/src/Rodin/Variational/P1/Jacobian.h b/src/Rodin/Variational/P1/Jacobian.h index 6321f2e8d..de71dec17 100644 --- a/src/Rodin/Variational/P1/Jacobian.h +++ b/src/Rodin/Variational/P1/Jacobian.h @@ -51,9 +51,9 @@ namespace Rodin::Variational public: using FESType = P1, Mesh>; - using ScalarType = typename FormLanguage::Traits; + using ScalarType = typename FormLanguage::Traits::ScalarType; - using MatrixType = Math::SpatialMatrix; + using SpatialMatrixType = Math::SpatialMatrix; using OperandType = GridFunction; @@ -76,7 +76,7 @@ namespace Rodin::Variational : Parent(std::move(other)) {} - void interpolate(Math::SpatialMatrix& out, const Geometry::Point& p) const + void interpolate(SpatialMatrixType& out, const Geometry::Point& p) const { const auto& polytope = p.getPolytope(); const auto& d = polytope.getDimension(); @@ -92,7 +92,7 @@ namespace Rodin::Variational if (inc.size() == 1) { const auto& tracePolytope = mesh.getPolytope(meshDim, *inc.begin()); - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const auto rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -117,7 +117,7 @@ namespace Rodin::Variational const auto& tracePolytope = mesh.getPolytope(meshDim, idx); if (traceDomain.count(tracePolytope->getAttribute())) { - const Math::SpatialVector rc = tracePolytope->getTransformation().inverse(pc); + const auto rc = tracePolytope->getTransformation().inverse(pc); const Geometry::Point np(*tracePolytope, std::cref(rc), pc); interpolate(out, np); return; @@ -137,8 +137,8 @@ namespace Rodin::Variational const auto& vdim = fes.getVectorDimension(); const auto& fe = fes.getFiniteElement(d, i); const auto& rc = p.getReferenceCoordinates(); - Math::SpatialMatrix jacobian(vdim, d); - Math::SpatialMatrix res(vdim, d); + SpatialMatrixType jacobian(vdim, d); + SpatialMatrixType res(vdim, d); res.setZero(); for (size_t local = 0; local < fe.getCount(); local++) { diff --git a/src/Rodin/Variational/P1/P1.cpp b/src/Rodin/Variational/P1/P1.cpp index 7f171d40c..fa9648340 100644 --- a/src/Rodin/Variational/P1/P1.cpp +++ b/src/Rodin/Variational/P1/P1.cpp @@ -8,6 +8,26 @@ namespace Rodin::Variational { + const Geometry::GeometryIndexed> + P1>::s_elements = + { + { Geometry::Polytope::Type::Point, P1Element(Geometry::Polytope::Type::Point) }, + { Geometry::Polytope::Type::Segment, P1Element(Geometry::Polytope::Type::Segment) }, + { Geometry::Polytope::Type::Triangle, P1Element(Geometry::Polytope::Type::Triangle) }, + { Geometry::Polytope::Type::Quadrilateral, P1Element(Geometry::Polytope::Type::Quadrilateral) }, + { Geometry::Polytope::Type::Tetrahedron, P1Element(Geometry::Polytope::Type::Tetrahedron) } + }; + + const Geometry::GeometryIndexed> + P1>::s_elements = + { + { Geometry::Polytope::Type::Point, P1Element(Geometry::Polytope::Type::Point) }, + { Geometry::Polytope::Type::Segment, P1Element(Geometry::Polytope::Type::Segment) }, + { Geometry::Polytope::Type::Triangle, P1Element(Geometry::Polytope::Type::Triangle) }, + { Geometry::Polytope::Type::Quadrilateral, P1Element(Geometry::Polytope::Type::Quadrilateral) }, + { Geometry::Polytope::Type::Tetrahedron, P1Element(Geometry::Polytope::Type::Tetrahedron) } + }; + P1, Geometry::Mesh> ::P1(const MeshType& mesh, size_t vdim) : m_mesh(mesh), m_vdim(vdim) diff --git a/src/Rodin/Variational/P1/P1.h b/src/Rodin/Variational/P1/P1.h index 335344d08..0b63b6aa5 100644 --- a/src/Rodin/Variational/P1/P1.h +++ b/src/Rodin/Variational/P1/P1.h @@ -22,21 +22,21 @@ namespace Rodin::FormLanguage { - template - struct Traits> + template + struct Traits> { using MeshType = Mesh; - using ScalarType = Number; + using ScalarType = Scalar; using RangeType = ScalarType; using ContextType = typename FormLanguage::Traits::ContextType; using ElementType = Variational::P1Element; }; - template - struct Traits, Mesh>> + template + struct Traits, Mesh>> { using MeshType = Mesh; - using ScalarType = Number; + using ScalarType = Scalar; using RangeType = Math::Vector; using ContextType = typename FormLanguage::Traits::ContextType; using ElementType = Variational::P1Element; @@ -67,16 +67,16 @@ namespace Rodin::Variational * This class is scalar valued, i.e. evaluations of the function are of * Rodin::Real type. */ - template - class P1> final - : public FiniteElementSpace>> + template <> + class P1> final + : public FiniteElementSpace>> { using KeyLeft = std::tuple; using KeyRight = Index; using IndexMap = FlatMap; public: - using ScalarType = Number; + using ScalarType = Real; /// Range type of value using RangeType = ScalarType; @@ -184,7 +184,7 @@ namespace Rodin::Variational P1& operator=(P1&& other) = default; inline - const P1Element& getFiniteElement(size_t d, Index i) const + const ElementType& getFiniteElement(size_t d, Index i) const { return s_elements[getMesh().getGeometry(d, i)]; } @@ -285,6 +285,232 @@ namespace Rodin::Variational template using RealP1 = P1; + /** + * @ingroup P1Specializations + * @brief Real valued Lagrange finite element space + * + * Represents the finite element space composed of scalar valued continuous, + * piecewise linear functions: + * @f[ + * \mathbb{P}_1 (\mathcal{T}_h) = \{ v \in C^0(\mathcal{T}_h) : v|_{\tau} \in \mathbb{P}_1(\tau), \ \tau \in \mathcal{T}_h \} \ . + * @f] + * + * This class is scalar valued, i.e. evaluations of the function are of + * Rodin::Complex type. + */ + template <> + class P1> final + : public FiniteElementSpace>> + { + using KeyLeft = std::tuple; + using KeyRight = Index; + using IndexMap = FlatMap; + + public: + using ScalarType = Complex; + + /// Range type of value + using RangeType = ScalarType; + + /// Represents the Context of the P1 space + using ContextType = Context::Sequential; + + /// Type of mesh on which the finite element space is built + using MeshType = Geometry::Mesh; + + /// Type of finite element + using ElementType = P1Element; + + /// Parent class + using Parent = FiniteElementSpace>; + + /** + * @brief Mapping for the scalar/complex P1 space. + */ + template + class Mapping : public FiniteElementSpaceMappingBase> + { + public: + using FunctionType = FunctionBase; + + Mapping(const Geometry::Polytope& polytope, const FunctionType& v) + : m_polytope(polytope), m_trans(m_polytope.getTransformation()), m_v(v.copy()) + {} + + Mapping(const Mapping&) = default; + + inline + auto operator()(const Math::SpatialVector& r) const + { + const Geometry::Point p(m_polytope, m_trans.get(), r); + return getFunction()(p); + } + + inline + constexpr + const FunctionType& getFunction() const + { + assert(m_v); + return *m_v; + } + + private: + Geometry::Polytope m_polytope; + std::reference_wrapper m_trans; + std::unique_ptr m_v; + }; + + /** + * @brief Inverse mapping for the scalar/complex P1 space. + */ + template + class InverseMapping : public FiniteElementSpaceInverseMappingBase> + { + public: + using FunctionType = CallableType; + + /** + * @param[in] polytope Reference to polytope on the mesh. + * @param[in] v Reference to the function defined on the reference + * space. + */ + InverseMapping(const FunctionType& v) + : m_v(v) + {} + + InverseMapping(const InverseMapping&) = default; + + inline + constexpr + auto operator()(const Geometry::Point& p) const + { + return getFunction()(p.getReferenceCoordinates()); + } + + inline + constexpr + const FunctionType& getFunction() const + { + return m_v.get(); + } + + private: + std::reference_wrapper m_v; + }; + + P1(const MeshType& mesh) + : m_mesh(mesh) + {} + + P1(const P1& other) + : Parent(other), + m_mesh(other.m_mesh) + {} + + P1(P1&& other) + : Parent(std::move(other)), + m_mesh(other.m_mesh) + {} + + P1& operator=(P1&& other) = default; + + inline + const ElementType& getFiniteElement(size_t d, Index i) const + { + return s_elements[getMesh().getGeometry(d, i)]; + } + + inline + size_t getSize() const override + { + return 2 * m_mesh.get().getVertexCount(); + } + + inline + size_t getVectorDimension() const override + { + return 1; + } + + inline + const MeshType& getMesh() const override + { + return m_mesh.get(); + } + + inline + const IndexArray& getDOFs(size_t d, Index i) const override + { + return getMesh().getConnectivity().getPolytope(d, i); + } + + inline + Index getGlobalIndex(const std::pair& idx, Index local) const override + { + const auto [d, i] = idx; + const auto& p = getMesh().getConnectivity().getPolytope(d, i); + const size_t q = local / 2; + const size_t r = local % 2; + assert(q < static_cast(p.size())); + return p(q) + r * m_mesh.get().getVertexCount(); + } + + /** + * @brief Returns the mapping of the function from the physical element + * to the reference element. + * @param[in] idx Index of the element in the mesh + * @param[in] v Function defined on an element of the mesh + * + * For all @f$ \tau \in \mathcal{T}_h @f$ in the mesh, the finite + * element space is generated by the bijective mapping: + * @f[ + * \psi_\tau : \mathbb{P}_1(\tau) \rightarrow \mathbb{P}_1(R) + * @f] + * taking a function @f$ v \in V(\tau) @f$ from the global element @f$ + * \tau @f$ element @f$ R @f$. + */ + template + inline + auto getMapping(const std::pair& idx, const FunctionBase& v) const + { + const auto [d, i] = idx; + const auto& mesh = getMesh(); + return Mapping(*mesh.getPolytope(d, i), v); + } + + template + inline + auto getMapping(const Geometry::Polytope& polytope, const FunctionBase& v) const + { + return Mapping(polytope, v); + } + + /** + * @brief Returns the inverse mapping of the function from the physical + * element to the reference element. + * @param[in] idx Index of the element in the mesh. + * @param[in] v Callable type + */ + template + inline + auto getInverseMapping(const std::pair& idx, const CallableType& v) const + { + return InverseMapping(v); + } + + template + inline + auto getInverseMapping(const Geometry::Polytope& polytope, const CallableType& v) const + { + return InverseMapping(v); + } + + private: + static const Geometry::GeometryIndexed> s_elements; + + std::reference_wrapper m_mesh; + }; + template using ComplexP1 = P1; @@ -488,17 +714,6 @@ namespace Rodin::Variational std::array, RODIN_P1_MAX_VECTOR_DIMENSION> initVectorP1Elements(); } - - template - const Geometry::GeometryIndexed> - P1>::s_elements = - { - { Geometry::Polytope::Type::Point, P1Element(Geometry::Polytope::Type::Point) }, - { Geometry::Polytope::Type::Segment, P1Element(Geometry::Polytope::Type::Segment) }, - { Geometry::Polytope::Type::Triangle, P1Element(Geometry::Polytope::Type::Triangle) }, - { Geometry::Polytope::Type::Quadrilateral, P1Element(Geometry::Polytope::Type::Quadrilateral) }, - { Geometry::Polytope::Type::Tetrahedron, P1Element(Geometry::Polytope::Type::Tetrahedron) } - }; } #endif diff --git a/src/Rodin/Variational/P1/P1Element.cpp b/src/Rodin/Variational/P1/P1Element.cpp index 38212d2e0..46d78efa0 100644 --- a/src/Rodin/Variational/P1/P1Element.cpp +++ b/src/Rodin/Variational/P1/P1Element.cpp @@ -2,8 +2,26 @@ namespace Rodin::Variational { - const Geometry::GeometryIndexed> - RealP1Element::s_ls = + const Geometry::GeometryIndexed RealP1Element::s_nodes = + { + { Geometry::Polytope::Type::Point, + Math::PointMatrix{{0}} }, + { Geometry::Polytope::Type::Segment, + Math::PointMatrix{{0, 1}} }, + { Geometry::Polytope::Type::Triangle, + Math::PointMatrix{{0, 1, 0}, + {0, 0, 1}} }, + { Geometry::Polytope::Type::Quadrilateral, + Math::PointMatrix{{0, 1, 0, 1}, + {0, 0, 1, 1}} }, + { Geometry::Polytope::Type::Tetrahedron, + Math::PointMatrix{{0, 1, 0, 0}, + {0, 0, 1, 0}, + {0, 0, 0, 1}} }, + }; + + const Geometry::GeometryIndexed> + RealP1Element::s_basis = { { Geometry::Polytope::Type::Point, { @@ -41,8 +59,8 @@ namespace Rodin::Variational } }; - const Geometry::GeometryIndexed> - RealP1Element::s_basis = + const Geometry::GeometryIndexed> + RealP1Element::s_gradient = { { Geometry::Polytope::Type::Point, { @@ -80,8 +98,8 @@ namespace Rodin::Variational } }; - const Geometry::GeometryIndexed> - RealP1Element::s_gradient = + const Geometry::GeometryIndexed> + RealP1Element::s_ls = { { Geometry::Polytope::Type::Point, { @@ -119,24 +137,6 @@ namespace Rodin::Variational } }; - const Geometry::GeometryIndexed RealP1Element::s_nodes = - { - { Geometry::Polytope::Type::Point, - Math::PointMatrix{{0}} }, - { Geometry::Polytope::Type::Segment, - Math::PointMatrix{{0, 1}} }, - { Geometry::Polytope::Type::Triangle, - Math::PointMatrix{{0, 1, 0}, - {0, 0, 1}} }, - { Geometry::Polytope::Type::Quadrilateral, - Math::PointMatrix{{0, 1, 0, 1}, - {0, 0, 1, 1}} }, - { Geometry::Polytope::Type::Tetrahedron, - Math::PointMatrix{{0, 1, 0, 0}, - {0, 0, 1, 0}, - {0, 0, 0, 1}} }, - }; - const Geometry::GeometryIndexed ComplexP1Element::s_nodes = { { Geometry::Polytope::Type::Point, @@ -155,6 +155,165 @@ namespace Rodin::Variational {0, 0, 0, 0, 0, 0, 1, 1}} }, }; + const Geometry::GeometryIndexed> + ComplexP1Element::s_basis = + { + { Geometry::Polytope::Type::Point, + { + { 0, Geometry::Polytope::Type::Point }, + { 1, Geometry::Polytope::Type::Point } + } + }, + { Geometry::Polytope::Type::Segment, + { + { 0, Geometry::Polytope::Type::Segment }, + { 1, Geometry::Polytope::Type::Segment }, + { 2, Geometry::Polytope::Type::Segment }, + { 3, Geometry::Polytope::Type::Segment } + } + }, + { Geometry::Polytope::Type::Triangle, + { + { 0, Geometry::Polytope::Type::Triangle }, + { 1, Geometry::Polytope::Type::Triangle }, + { 2, Geometry::Polytope::Type::Triangle }, + { 3, Geometry::Polytope::Type::Triangle }, + { 4, Geometry::Polytope::Type::Triangle }, + { 5, Geometry::Polytope::Type::Triangle } + } + }, + { Geometry::Polytope::Type::Quadrilateral, + { + { 0, Geometry::Polytope::Type::Quadrilateral }, + { 1, Geometry::Polytope::Type::Quadrilateral }, + { 2, Geometry::Polytope::Type::Quadrilateral }, + { 3, Geometry::Polytope::Type::Quadrilateral }, + { 4, Geometry::Polytope::Type::Quadrilateral }, + { 5, Geometry::Polytope::Type::Quadrilateral }, + { 6, Geometry::Polytope::Type::Quadrilateral }, + { 7, Geometry::Polytope::Type::Quadrilateral } + } + }, + { Geometry::Polytope::Type::Tetrahedron, + { + { 0, Geometry::Polytope::Type::Tetrahedron }, + { 1, Geometry::Polytope::Type::Tetrahedron }, + { 2, Geometry::Polytope::Type::Tetrahedron }, + { 3, Geometry::Polytope::Type::Tetrahedron }, + { 4, Geometry::Polytope::Type::Tetrahedron }, + { 5, Geometry::Polytope::Type::Tetrahedron }, + { 6, Geometry::Polytope::Type::Tetrahedron }, + { 7, Geometry::Polytope::Type::Tetrahedron } + } + } + }; + + const Geometry::GeometryIndexed> + ComplexP1Element::s_gradient = + { + { Geometry::Polytope::Type::Point, + { + { 0, Geometry::Polytope::Type::Point }, + { 1, Geometry::Polytope::Type::Point } + } + }, + { Geometry::Polytope::Type::Segment, + { + { 0, Geometry::Polytope::Type::Segment }, + { 1, Geometry::Polytope::Type::Segment }, + { 2, Geometry::Polytope::Type::Segment }, + { 3, Geometry::Polytope::Type::Segment } + } + }, + { Geometry::Polytope::Type::Triangle, + { + { 0, Geometry::Polytope::Type::Triangle }, + { 1, Geometry::Polytope::Type::Triangle }, + { 2, Geometry::Polytope::Type::Triangle }, + { 3, Geometry::Polytope::Type::Triangle }, + { 4, Geometry::Polytope::Type::Triangle }, + { 5, Geometry::Polytope::Type::Triangle } + } + }, + { Geometry::Polytope::Type::Quadrilateral, + { + { 0, Geometry::Polytope::Type::Quadrilateral }, + { 1, Geometry::Polytope::Type::Quadrilateral }, + { 2, Geometry::Polytope::Type::Quadrilateral }, + { 3, Geometry::Polytope::Type::Quadrilateral }, + { 4, Geometry::Polytope::Type::Quadrilateral }, + { 5, Geometry::Polytope::Type::Quadrilateral }, + { 6, Geometry::Polytope::Type::Quadrilateral }, + { 7, Geometry::Polytope::Type::Quadrilateral } + } + }, + { Geometry::Polytope::Type::Tetrahedron, + { + { 0, Geometry::Polytope::Type::Tetrahedron }, + { 1, Geometry::Polytope::Type::Tetrahedron }, + { 2, Geometry::Polytope::Type::Tetrahedron }, + { 3, Geometry::Polytope::Type::Tetrahedron }, + { 4, Geometry::Polytope::Type::Tetrahedron }, + { 5, Geometry::Polytope::Type::Tetrahedron }, + { 6, Geometry::Polytope::Type::Tetrahedron }, + { 7, Geometry::Polytope::Type::Tetrahedron } + } + } + }; + + const Geometry::GeometryIndexed> + ComplexP1Element::s_ls = + { + { Geometry::Polytope::Type::Point, + { + { 0, Geometry::Polytope::Type::Point }, + { 1, Geometry::Polytope::Type::Point } + } + }, + { Geometry::Polytope::Type::Segment, + { + { 0, Geometry::Polytope::Type::Segment }, + { 1, Geometry::Polytope::Type::Segment }, + { 2, Geometry::Polytope::Type::Segment }, + { 3, Geometry::Polytope::Type::Segment } + } + }, + { Geometry::Polytope::Type::Triangle, + { + { 0, Geometry::Polytope::Type::Triangle }, + { 1, Geometry::Polytope::Type::Triangle }, + { 2, Geometry::Polytope::Type::Triangle }, + { 3, Geometry::Polytope::Type::Triangle }, + { 4, Geometry::Polytope::Type::Triangle }, + { 5, Geometry::Polytope::Type::Triangle } + } + }, + { Geometry::Polytope::Type::Quadrilateral, + { + { 0, Geometry::Polytope::Type::Quadrilateral }, + { 1, Geometry::Polytope::Type::Quadrilateral }, + { 2, Geometry::Polytope::Type::Quadrilateral }, + { 3, Geometry::Polytope::Type::Quadrilateral }, + { 4, Geometry::Polytope::Type::Quadrilateral }, + { 5, Geometry::Polytope::Type::Quadrilateral }, + { 6, Geometry::Polytope::Type::Quadrilateral }, + { 7, Geometry::Polytope::Type::Quadrilateral } + } + }, + { Geometry::Polytope::Type::Tetrahedron, + { + { 0, Geometry::Polytope::Type::Tetrahedron }, + { 1, Geometry::Polytope::Type::Tetrahedron }, + { 2, Geometry::Polytope::Type::Tetrahedron }, + { 3, Geometry::Polytope::Type::Tetrahedron }, + { 4, Geometry::Polytope::Type::Tetrahedron }, + { 5, Geometry::Polytope::Type::Tetrahedron }, + { 6, Geometry::Polytope::Type::Tetrahedron }, + { 7, Geometry::Polytope::Type::Tetrahedron } + } + } + }; + const std::array, RODIN_P1_MAX_VECTOR_DIMENSION> VectorP1Element::s_nodes = Internal::initVectorP1Nodes(); diff --git a/src/Rodin/Variational/P1/P1Element.h b/src/Rodin/Variational/P1/P1Element.h index 84a5700bf..4d1e525a5 100644 --- a/src/Rodin/Variational/P1/P1Element.h +++ b/src/Rodin/Variational/P1/P1Element.h @@ -292,7 +292,7 @@ namespace Rodin::Variational LinearForm(size_t i, Geometry::Polytope::Type g) : m_i(i), m_g(g) { - assert(i < Geometry::Polytope::getVertexCount(g)); + assert(i < 2 * Geometry::Polytope::getVertexCount(g)); } constexpr @@ -312,7 +312,10 @@ namespace Rodin::Variational constexpr auto operator()(const T& v) const { - return v(s_nodes[m_g].col(m_i)); + if (m_i % 2 == 0) + return v(s_nodes[m_g].col(m_i)).real(); + else + return v(s_nodes[m_g].col(m_i)).imag(); } private: @@ -332,7 +335,7 @@ namespace Rodin::Variational BasisFunction(size_t i, Geometry::Polytope::Type g) : m_i(i), m_g(g) { - assert(i < Geometry::Polytope::getVertexCount(g)); + assert(i < 2 * Geometry::Polytope::getVertexCount(g)); } constexpr @@ -363,7 +366,7 @@ namespace Rodin::Variational GradientFunction(size_t i, Geometry::Polytope::Type g) : m_i(i), m_g(g) { - assert(i < Geometry::Polytope::getVertexCount(g)); + assert(i < 2 * Geometry::Polytope::getVertexCount(g)); } constexpr diff --git a/src/Rodin/Variational/Pow.h b/src/Rodin/Variational/Pow.h index f7308a9fa..68d2c52a9 100644 --- a/src/Rodin/Variational/Pow.h +++ b/src/Rodin/Variational/Pow.h @@ -96,8 +96,17 @@ namespace Rodin::Variational constexpr auto getValue(const Geometry::Point& p) const { - assert(m_s); - return Math::pow(m_s->getValue(p), m_p); + return Math::pow(getBase().getValue(p), getExponent()); + } + + const BaseType& getBase() const + { + return *m_s; + } + + const ExponentType& getExponent() const + { + return m_p; } inline Pow* copy() const noexcept override diff --git a/src/Rodin/Variational/Problem.h b/src/Rodin/Variational/Problem.h index 0e3169d75..0e2f4cdea 100644 --- a/src/Rodin/Variational/Problem.h +++ b/src/Rodin/Variational/Problem.h @@ -45,7 +45,7 @@ namespace Rodin::Variational /** * @brief Abstract base class for variational problems. */ - template + template class ProblemBase : public FormLanguage::Base { public: @@ -57,13 +57,15 @@ namespace Rodin::Variational using OperatorScalarType = typename FormLanguage::Traits::ScalarType; + using ScalarType = Scalar; + ProblemBase() = default; ProblemBase(ProblemBase&& other) = default; ProblemBase(const ProblemBase& other) = default; - virtual ProblemBase& operator=(const ProblemBody& rhs) = 0; + virtual ProblemBase& operator=(const ProblemBody& rhs) = 0; virtual void solve(Solver::SolverBase& solver) = 0; @@ -112,16 +114,22 @@ namespace Rodin::Variational class Problem< TrialFES, TestFES, Math::SparseMatrix< - decltype( - std::declval::ScalarType>() * - std::declval::ScalarType>())>, + typename FormLanguage::Mult< + typename FormLanguage::Traits::ScalarType, + typename FormLanguage::Traits::ScalarType> + ::Type>, Math::Vector::ScalarType>> : public ProblemBase< Math::SparseMatrix< - decltype( - std::declval::ScalarType>() * - std::declval::ScalarType>())>, - Math::Vector::ScalarType>> + typename FormLanguage::Mult< + typename FormLanguage::Traits::ScalarType, + typename FormLanguage::Traits::ScalarType> + ::Type>, + Math::Vector::ScalarType>, + typename FormLanguage::Mult< + typename FormLanguage::Traits::ScalarType, + typename FormLanguage::Traits::ScalarType> + ::Type> { public: using TrialFESScalarType = @@ -131,10 +139,12 @@ namespace Rodin::Variational typename FormLanguage::Traits::ScalarType; using OperatorScalarType = - decltype(std::declval() * std::declval()); + typename FormLanguage::Mult::Type; using VectorScalarType = TestFESScalarType; + using ScalarType = OperatorScalarType; + using ContextType = Context::Sequential; using OperatorType = Math::SparseMatrix; @@ -143,7 +153,7 @@ namespace Rodin::Variational using LinearFormIntegratorBaseType = LinearFormIntegratorBase; - using Parent = ProblemBase; + using Parent = ProblemBase; /** * @brief Constructs an empty problem involving the trial function @f$ u @f$ @@ -420,7 +430,7 @@ namespace Rodin::Variational getTrialFunction().getSolution().setWeights(std::move(m_guess)); } - Problem& operator=(const ProblemBody& rhs) override + Problem& operator=(const ProblemBody& rhs) override { for (auto& bfi : rhs.getLocalBFIs()) m_bilinearForm.add(bfi); @@ -483,15 +493,20 @@ namespace Rodin::Variational OperatorType m_stiffness; }; - template Problem(TrialFunction&, TestFunction&) - -> Problem, Math::Vector>; + -> Problem::ScalarType, + typename FormLanguage::Traits::ScalarType>::Type>, + Math::Vector< + typename FormLanguage::Traits::ScalarType>>; template class Problem< Tuple, Math::SparseMatrix, Math::Vector> - : public ProblemBase, Math::Vector> + : public ProblemBase, Math::Vector, Real> { template @@ -511,7 +526,7 @@ namespace Rodin::Variational using VectorType = Math::Vector; - using Parent = ProblemBase; + using Parent = ProblemBase; private: template @@ -750,7 +765,7 @@ namespace Rodin::Variational }); } - Problem& operator=(const ProblemBody& rhs) override + Problem& operator=(const ProblemBody& rhs) override { for (auto& bfi : rhs.getLocalBFIs()) { diff --git a/src/Rodin/Variational/ProblemBody.h b/src/Rodin/Variational/ProblemBody.h index 66a5a130d..cf991e351 100644 --- a/src/Rodin/Variational/ProblemBody.h +++ b/src/Rodin/Variational/ProblemBody.h @@ -28,10 +28,11 @@ namespace Rodin::Variational /** * @brief Represents the body of a variational problem. */ + template class ProblemBodyBase : public FormLanguage::Base { public: - using ScalarType = Real; + using ScalarType = Scalar; using LinearFormIntegratorBaseType = LinearFormIntegratorBase; @@ -67,7 +68,6 @@ namespace Rodin::Variational m_periodicBdr(std::move(other.m_periodicBdr)) {} - inline ProblemBodyBase& operator=(ProblemBodyBase&& other) { m_lfis = std::move(other.m_lfis); @@ -142,12 +142,12 @@ namespace Rodin::Variational PeriodicBoundary m_periodicBdr; }; - template <> - class ProblemBody - : public ProblemBodyBase + template + class ProblemBody + : public ProblemBodyBase { public: - using Parent = ProblemBodyBase; + using Parent = ProblemBodyBase; ProblemBody() = default; @@ -165,8 +165,9 @@ namespace Rodin::Variational } }; - template - class ProblemBody : public ProblemBodyBase + template + class ProblemBody + : public ProblemBodyBase { public: using OperatorType = Operator; @@ -177,7 +178,7 @@ namespace Rodin::Variational using BilinearFormBaseListType = FormLanguage::List; - using Parent = ProblemBodyBase; + using Parent = ProblemBodyBase; ProblemBody() = default; @@ -210,8 +211,8 @@ namespace Rodin::Variational BilinearFormBaseListType m_bfs; }; - template - class ProblemBody : public ProblemBodyBase + template + class ProblemBody : public ProblemBodyBase { public: using VectorType = Vector; @@ -220,7 +221,7 @@ namespace Rodin::Variational using LinearFormBaseListType = FormLanguage::List; - using Parent = ProblemBodyBase; + using Parent = ProblemBodyBase; ProblemBody() = default; @@ -234,13 +235,11 @@ namespace Rodin::Variational m_lfs(std::move(other.m_lfs)) {} - inline LinearFormBaseListType& getLFs() { return m_lfs; } - inline const LinearFormBaseListType& getLFs() const { return m_lfs; @@ -255,8 +254,8 @@ namespace Rodin::Variational LinearFormBaseListType m_lfs; }; - template - class ProblemBody : public ProblemBodyBase + template + class ProblemBody : public ProblemBodyBase { public: using VectorType = Vector; @@ -287,7 +286,7 @@ namespace Rodin::Variational using GlobalBilinearFormIntegratorBaseListType = FormLanguage::List; - using Parent = ProblemBodyBase; + using Parent = ProblemBodyBase; ProblemBody(const LocalBilinearFormIntegratorBaseType& bfi) { @@ -309,19 +308,19 @@ namespace Rodin::Variational this->getGlobalBFIs().add(bfis); } - ProblemBody(const ProblemBody& pbo) + ProblemBody(const ProblemBody& pbo) : Parent(pbo) { m_bfs.add(pbo.getBFs()); } - ProblemBody(const ProblemBody& pbv) + ProblemBody(const ProblemBody& pbv) : Parent(pbv) { m_lfs.add(pbv.getLFs()); } - ProblemBody(const ProblemBody& parent) + ProblemBody(const ProblemBody& parent) : Parent(parent) {} @@ -337,25 +336,21 @@ namespace Rodin::Variational m_bfs(std::move(other.m_bfs)) {} - inline LinearFormBaseListType& getLFs() { return m_lfs; } - inline BilinearFormBaseListType& getBFs() { return m_bfs; } - inline const LinearFormBaseListType& getLFs() const { return m_lfs; } - inline const BilinearFormBaseListType& getBFs() const { return m_bfs; @@ -371,403 +366,430 @@ namespace Rodin::Variational BilinearFormBaseListType m_bfs; }; - ProblemBody(const LocalBilinearFormIntegratorBase&) - -> ProblemBody; + template + ProblemBody(const LocalBilinearFormIntegratorBase&) + -> ProblemBody; - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, const LinearFormIntegratorBase& lfi) + template + auto + operator+(const LocalBilinearFormIntegratorBase& bfi, const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getLFIs().add(lfi); return res; } - template - inline - ProblemBody operator+( - const LinearFormIntegratorBase& lfi, const LocalBilinearFormIntegratorBase& bfi) + template + auto + operator+(const LinearFormIntegratorBase& lfi, const LocalBilinearFormIntegratorBase& bfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getLFIs().add(lfi); return res; } - template - inline - ProblemBody operator-( - const LocalBilinearFormIntegratorBase& bfi, const LinearFormIntegratorBase& lfi) + template + auto + operator-(const LocalBilinearFormIntegratorBase& bfi, const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getLFIs().add(UnaryMinus(lfi)); return res; } - template - inline - ProblemBody operator-( - const GlobalBilinearFormIntegratorBase& bfi, const LinearFormIntegratorBase& lfi) + template + auto + operator-(const GlobalBilinearFormIntegratorBase& bfi, const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getGlobalBFIs().add(bfi); res.getLFIs().add(UnaryMinus(lfi)); return res; } - template - inline - ProblemBody operator-( - const GlobalBilinearFormIntegratorBase& bfi, - const FormLanguage::List>& lfis) + template + auto + operator-( + const GlobalBilinearFormIntegratorBase& bfi, + const FormLanguage::List>& lfis) { - ProblemBody res; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getGlobalBFIs().add(bfi); res.getLFIs().add(UnaryMinus(lfis)); return res; } - template - inline - ProblemBody operator-( - const LinearFormIntegratorBase& lfi, const LocalBilinearFormIntegratorBase& bfi) + template + auto operator-( + const LinearFormIntegratorBase& lfi, const LocalBilinearFormIntegratorBase& bfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getLocalBFIs().add(UnaryMinus(bfi)); res.getLFIs().add(lfi); return res; } - template - inline - ProblemBody operator+( - const FormLanguage::List>& bfis, - const LinearFormIntegratorBase& lfi) + template + auto operator+( + const FormLanguage::List>& bfis, + const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(bfis); res.getLFIs().add(lfi); return res; } - template - inline - ProblemBody operator+( - const FormLanguage::List>& lbfis, - const GlobalBilinearFormIntegratorBase& gbfi) + template + auto operator+( + const FormLanguage::List>& lbfis, + const GlobalBilinearFormIntegratorBase& gbfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(lbfis); res.getGlobalBFIs().add(gbfi); return res; } - template - inline - ProblemBody operator+( - const FormLanguage::List>& lbfis, - const FormLanguage::List>& gbfis) + template + auto operator+( + const FormLanguage::List>& lbfis, + const FormLanguage::List>& gbfis) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(lbfis); res.getGlobalBFIs().add(gbfis); return res; } - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& lbfi, - const FormLanguage::List>& gbfis) + template + auto operator+( + const LocalBilinearFormIntegratorBase& lbfi, + const FormLanguage::List>& gbfis) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(lbfi); res.getGlobalBFIs().add(gbfis); return res; } - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& lbfi, - const GlobalBilinearFormIntegratorBase& gbfi) + template + auto operator+( + const LocalBilinearFormIntegratorBase& lbfi, + const GlobalBilinearFormIntegratorBase& gbfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(lbfi); res.getGlobalBFIs().add(gbfi); return res; } - template - inline - ProblemBody operator-( - const FormLanguage::List>& bfis, - const LinearFormIntegratorBase& lfi) + template + auto operator-( + const FormLanguage::List>& bfis, + const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getLocalBFIs().add(bfis); res.getLFIs().add(UnaryMinus(lfi)); return res; } - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, const DirichletBCBase& dbc) + template + auto + operator+(const LocalBilinearFormIntegratorBase& bfi, const DirichletBCBase& dbc) { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getDBCs().add(dbc); return res; } - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, const FormLanguage::List& dbcs) + template + auto + operator+(const LocalBilinearFormIntegratorBase& bfi, const FormLanguage::List& dbcs) { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getDBCs().add(dbcs); return res; } - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, const PeriodicBCBase& pbc) + template + auto operator+( + const LocalBilinearFormIntegratorBase& bfi, const PeriodicBCBase& pbc) { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getPBCs().add(pbc); return res; } - template - inline - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, const FormLanguage::List& pbcs) + template + auto + operator+(const LocalBilinearFormIntegratorBase& bfi, const FormLanguage::List& pbcs) { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getPBCs().add(pbcs); return res; } - template - inline - ProblemBody operator+( - const FormLanguage::List>& bfis, const DirichletBCBase& dbc) + template + auto operator+( + const FormLanguage::List>& bfis, const DirichletBCBase& dbc) { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfis); res.getDBCs().add(dbc); return res; } - template - inline - ProblemBody operator+( - const FormLanguage::List>& bfis, const FormLanguage::List& dbcs) + template + auto + operator+( + const FormLanguage::List>& bfis, + const FormLanguage::List& dbcs) { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfis); res.getDBCs().add(dbcs); return res; } - template - inline - ProblemBody operator+( - const FormLanguage::List>& bfis, const PeriodicBCBase& pbc) + template + auto operator+( + const FormLanguage::List>& bfis, + const FormLanguage::List& pbcs) { - ProblemBody res; - res.getLocalBFIs().add(bfis); - res.getPBCs().add(pbc); - return res; - } - - template - inline - ProblemBody operator+( - const FormLanguage::List>& bfis, const FormLanguage::List& pbcs) - { - ProblemBody res; + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res; res.getLocalBFIs().add(bfis); res.getPBCs().add(pbcs); return res; } - template - inline - ProblemBody operator+( - const ProblemBody& pb, - const LinearFormIntegratorBase& lfi) + template + auto operator+( + const ProblemBody& pb, + const LinearFormIntegratorBase& lfi) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res(pb); res.getLFIs().add(lfi); return res; } - template - inline - ProblemBody operator+( - const ProblemBody& pb, - const LocalBilinearFormIntegratorBase& lbfi) + template + auto operator+( + const ProblemBody& pb, + const LocalBilinearFormIntegratorBase& lbfi) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res(pb); res.getLocalBFIs().add(lbfi); return res; } - template - inline - ProblemBody operator-( - const ProblemBody& pb, - const LocalBilinearFormIntegratorBase& lbfi) + template + auto operator-( + const ProblemBody& pb, + const LocalBilinearFormIntegratorBase& lbfi) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res(pb); res.getLocalBFIs().add(UnaryMinus(lbfi)); return res; } - template - inline - ProblemBody operator+( - const ProblemBody& pb, - const GlobalBilinearFormIntegratorBase& gbfi) + template + auto + operator+( + const ProblemBody& pb, + const GlobalBilinearFormIntegratorBase& gbfi) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res(pb); res.getGlobalBFIs().add(gbfi); return res; } - template - ProblemBody operator+( - const ProblemBody& pb, - const LinearFormIntegratorBase& lfi) - { - ProblemBody res(pb); - res.getLFIs().add(lfi); - return res; - } - - template - ProblemBody operator+( - const ProblemBody& pb, - const FormLanguage::List>& lfis) + template + auto + operator+( + const ProblemBody& pb, + const FormLanguage::List>& lfis) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res(pb); res.getLFIs().add(lfis); return res; } - template - ProblemBody operator-( - const ProblemBody& pb, - const LinearFormIntegratorBase& lfi) + template + auto + operator-( + const ProblemBody& pb, + const LinearFormIntegratorBase& lfi) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res(pb); res.getLFIs().add(UnaryMinus(lfi)); return res; } - template - ProblemBody operator-( - const ProblemBody& pb, - const FormLanguage::List>& lfis) + template + auto + operator-( + const ProblemBody& pb, + const FormLanguage::List>& lfis) { - ProblemBody res(pb); + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res(pb); res.getLFIs().add(UnaryMinus(lfis)); return res; } - template - ProblemBody operator+( - const ProblemBody& pb, const DirichletBCBase& dbc) + template + auto + operator+( + const ProblemBody& pb, const DirichletBCBase& dbc) { - ProblemBody res(pb); + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res(pb); res.getDBCs().add(dbc); return res; } - template - ProblemBody operator+( - const ProblemBody& pb, const FormLanguage::List& dbcs) + template + auto + operator+( + const ProblemBody& pb, const FormLanguage::List& dbcs) { - ProblemBody res(pb); + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res(pb); res.getEssentialBoundary().add(dbcs); return res; } - template - ProblemBody operator+( - const ProblemBody& pb, const PeriodicBCBase& pbc) + template + auto + operator+( + const ProblemBody& pb, + const PeriodicBCBase& pbc) { - ProblemBody res(pb); + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res(pb); res.getPBCs().add(pbc); return res; } - template - ProblemBody operator+( - const ProblemBody& pb, + template + auto + operator+( + const ProblemBody& pb, const BilinearFormBase& bf) { - ProblemBody res(pb); + using RHSScalar = typename FormLanguage::Traits>::ScalarType; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res(pb); res.getBFs().add(bf); return res; } - template - ProblemBody operator+( - const ProblemBody& pb, const FormLanguage::List& pbcs) + template + auto + operator+( + const ProblemBody& pb, + const FormLanguage::List& pbcs) { - ProblemBody res(pb); + using RHSScalar = Real; + using ScalarType = RHSScalar; + ProblemBody res(pb); res.getPBCs().add(pbcs); return res; } - template - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, const BilinearFormBase& bf) + template + auto + operator+( + const LocalBilinearFormIntegratorBase& bfi, const BilinearFormBase& bf) { - ProblemBody res; + using RHSScalar = typename FormLanguage::Traits>::ScalarType; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getBFs().add(bf); return res; } - template - ProblemBody operator-( - const BilinearFormBase& bf, const LinearFormIntegratorBase& lfi) + template + auto + operator-( + const BilinearFormBase& bf, const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using LHSScalar = typename FormLanguage::Traits>::ScalarType; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getBFs().add(bf); res.getLFIs().add(UnaryMinus(lfi)); return res; } - template - ProblemBody operator-( + template + auto + operator-( const FormLanguage::List>& bfs, - const LinearFormIntegratorBase& lfi) + const LinearFormIntegratorBase& lfi) { - ProblemBody res; + using LHSScalar = typename FormLanguage::Traits>::ScalarType; + using ScalarType = typename FormLanguage::Minus::Type; + ProblemBody res; res.getBFs().add(bfs); res.getLFIs().add(UnaryMinus(lfi)); return res; } - template - ProblemBody operator+( - const LocalBilinearFormIntegratorBase& bfi, + template + auto + operator+( + const LocalBilinearFormIntegratorBase& bfi, const FormLanguage::List>& bfs) { - ProblemBody res; + using RHSScalar = typename FormLanguage::Traits>::ScalarType; + using ScalarType = typename FormLanguage::Sum::Type; + ProblemBody res; res.getLocalBFIs().add(bfi); res.getBFs().add(bfs); return res; diff --git a/src/Rodin/Variational/QuadratureRule.h b/src/Rodin/Variational/QuadratureRule.h index 5a64c1ef3..fae788e7a 100644 --- a/src/Rodin/Variational/QuadratureRule.h +++ b/src/Rodin/Variational/QuadratureRule.h @@ -271,7 +271,7 @@ namespace Rodin::Variational std::reference_wrapper> m_u; TestFunction m_v; - FlatSet m_attrs; + FlatSet m_attrs; LinearForm> m_lf; std::optional m_value; diff --git a/src/Rodin/Variational/RangeType.cpp b/src/Rodin/Variational/RangeType.cpp index 71a79e20d..fc251c2c1 100644 --- a/src/Rodin/Variational/RangeType.cpp +++ b/src/Rodin/Variational/RangeType.cpp @@ -19,7 +19,7 @@ namespace Rodin::Variational } case RangeType::Integer: { - os << "Boolean"; + os << "Integer"; break; } case RangeType::Real: @@ -27,6 +27,11 @@ namespace Rodin::Variational os << "Real"; break; } + case RangeType::Complex: + { + os << "Complex"; + break; + } case RangeType::Vector: { os << "Vector"; diff --git a/src/Rodin/Variational/RangeType.h b/src/Rodin/Variational/RangeType.h index d5be0faf4..75926539b 100644 --- a/src/Rodin/Variational/RangeType.h +++ b/src/Rodin/Variational/RangeType.h @@ -17,6 +17,7 @@ namespace Rodin::Variational Boolean, Integer, Real, + Complex, Vector, Matrix }; diff --git a/src/Rodin/Variational/Re.h b/src/Rodin/Variational/Re.h new file mode 100644 index 000000000..fe382649e --- /dev/null +++ b/src/Rodin/Variational/Re.h @@ -0,0 +1,62 @@ +#ifndef RODDIN_VARIATIONAL_RE_H +#define RODDIN_VARIATIONAL_RE_H + +#include "ForwardDecls.h" +#include "RealFunction.h" + +namespace Rodin::Variational +{ + template + class Re; + + template + class Re> : public RealFunctionBase>> + { + public: + using OperandType = FunctionBase; + + using ScalarType = Real; + + using Parent = RealFunctionBase>>; + + Re(const OperandType& f) + : m_operand(f.copy()) + {} + + Re(const Re& other) + : Parent(other), + m_operand(other.m_operand->copy()) + {} + + Re(Re&& other) + : Parent(std::move(other)), + m_operand(std::move(other.m_operand)) + {} + + const OperandType& getOperand() const + { + assert(m_operand); + return *m_operand; + } + + inline + constexpr + ScalarType getValue(const Geometry::Point& p) const + { + return getOperand().getValue(p).real(); + } + + inline Re* copy() const noexcept override + { + return new Re(*this); + } + + private: + std::unique_ptr m_operand; + }; + + template + Re(const FunctionBase&) -> Re>; +} + +#endif diff --git a/src/Rodin/Variational/RealFunction.h b/src/Rodin/Variational/RealFunction.h index 6ab41fcfd..315459239 100644 --- a/src/Rodin/Variational/RealFunction.h +++ b/src/Rodin/Variational/RealFunction.h @@ -81,11 +81,13 @@ namespace Rodin::Variational * @ingroup RealFunctionSpecializations */ template - class RealFunction> final - : public RealFunctionBase> + class RealFunction> final + : public RealFunctionBase> { public: - using Parent = RealFunctionBase>; + using ScalarType = Real; + + using Parent = RealFunctionBase>; RealFunction(const RealFunctionBase& nested) : m_nested(nested.copy()) @@ -103,7 +105,7 @@ namespace Rodin::Variational inline constexpr - auto getValue(const Geometry::Point& v) const + ScalarType getValue(const Geometry::Point& v) const { return m_nested->getValue(v); } @@ -122,7 +124,7 @@ namespace Rodin::Variational } private: - std::unique_ptr> m_nested; + std::unique_ptr> m_nested; }; /** @@ -140,13 +142,15 @@ namespace Rodin::Variational : public RealFunctionBase> { public: + using ScalarType = Real; + using Parent = RealFunctionBase>; /** * @brief Constructs a RealFunction from a constant scalar value. * @param[in] x Constant scalar value */ - RealFunction(Real x) + RealFunction(const Real& x) : m_x(x) {} @@ -169,14 +173,14 @@ namespace Rodin::Variational inline constexpr - Real getValue() const + const Real& getValue() const { return m_x; } inline constexpr - Real getValue(const Geometry::Point&) const + ScalarType getValue(const Geometry::Point&) const { return m_x; } @@ -200,6 +204,8 @@ namespace Rodin::Variational : public RealFunctionBase> { public: + using ScalarType = Real; + using Parent = RealFunctionBase>; /** @@ -229,16 +235,16 @@ namespace Rodin::Variational inline constexpr - Real getValue() const + const Integer& getValue() const { return m_x; } inline constexpr - Real getValue(const Geometry::Point&) const + ScalarType getValue(const Geometry::Point&) const { - return static_cast(m_x); + return m_x; } inline RealFunction* copy() const noexcept override @@ -262,6 +268,8 @@ namespace Rodin::Variational static_assert(std::is_invocable_r_v); public: + using ScalarType = Real; + using Parent = RealFunctionBase>; RealFunction(F f) @@ -287,7 +295,7 @@ namespace Rodin::Variational inline constexpr - Real getValue(const Geometry::Point& v) const + ScalarType getValue(const Geometry::Point& v) const { return m_f(v); } diff --git a/src/Rodin/Variational/ScalarFunction.h b/src/Rodin/Variational/ScalarFunction.h index 41667c620..570ae392d 100644 --- a/src/Rodin/Variational/ScalarFunction.h +++ b/src/Rodin/Variational/ScalarFunction.h @@ -86,7 +86,14 @@ namespace Rodin::Variational constexpr void getValue(ScalarType& res, const Geometry::Point& p) const { - res = this->getValue(p); + if constexpr (Internal::HasGetValueMethod::Value) + { + return static_cast(*this).getValue(res, p); + } + else + { + res = getValue(p); + } } virtual ScalarFunctionBase* copy() const noexcept override = 0; @@ -94,44 +101,54 @@ namespace Rodin::Variational /** * @ingroup ScalarFunctionSpecializations + * @brief Represents a constant scalar function with type Real. */ - template - class ScalarFunction> final - : public ScalarFunctionBase>> + template <> + class ScalarFunction final + : public ScalarFunctionBase> { public: - using ScalarType = Scalar; + using ScalarType = Real; - using Parent = - ScalarFunctionBase>>; + using Parent = ScalarFunctionBase>; - ScalarFunction(const ScalarFunctionBase& nested) - : m_nested(nested.copy()) + /** + * @brief Constructs a ScalarFunction from a constant scalar value. + * @param[in] x Constant scalar value + */ + ScalarFunction(const ScalarType& x) + : m_x(x) {} ScalarFunction(const ScalarFunction& other) : Parent(other), - m_nested(other.m_nested->copy()) + m_x(other.m_x) {} ScalarFunction(ScalarFunction&& other) : Parent(std::move(other)), - m_nested(std::move(other.m_nested)) + m_x(other.m_x) {} inline constexpr - auto getValue(const Geometry::Point& v) const + ScalarFunction& traceOf(Geometry::Attribute) { - return m_nested->getValue(v); + return *this; } inline constexpr - ScalarFunction& traceOf(Geometry::Attribute attrs) + const Real& getValue() const { - m_nested->traceOf(attrs); - return *this; + return m_x; + } + + inline + constexpr + ScalarType getValue(const Geometry::Point&) const + { + return m_x; } inline ScalarFunction* copy() const noexcept override @@ -140,28 +157,26 @@ namespace Rodin::Variational } private: - std::unique_ptr> m_nested; + const Real m_x; }; /** * @brief CTAD for ScalarFunction. */ - template - ScalarFunction(const ScalarFunctionBase&) - -> ScalarFunction>; + ScalarFunction(const Real&) -> ScalarFunction; /** * @ingroup ScalarFunctionSpecializations - * @brief Represents a constant scalar function with type Real. + * @brief Represents a constant scalar function with type Complex. */ template <> - class ScalarFunction final - : public ScalarFunctionBase> + class ScalarFunction final + : public ScalarFunctionBase> { public: - using ScalarType = Real; + using ScalarType = Complex; - using Parent = ScalarFunctionBase>; + using Parent = ScalarFunctionBase>; /** * @brief Constructs a ScalarFunction from a constant scalar value. @@ -190,7 +205,7 @@ namespace Rodin::Variational inline constexpr - ScalarType getValue() const + const Complex& getValue() const { return m_x; } @@ -208,64 +223,54 @@ namespace Rodin::Variational } private: - const ScalarType m_x; + const Complex m_x; }; /** * @brief CTAD for ScalarFunction. */ - ScalarFunction(const Real&) -> ScalarFunction; + ScalarFunction(const Complex&) -> ScalarFunction; /** * @ingroup ScalarFunctionSpecializations - * @brief Represents a constant scalar function with type Complex. */ - template <> - class ScalarFunction final - : public ScalarFunctionBase> + template + class ScalarFunction> final + : public ScalarFunctionBase>> { public: - using ScalarType = Complex; + using ScalarType = Scalar; - using Parent = ScalarFunctionBase>; + using Parent = + ScalarFunctionBase>>; - /** - * @brief Constructs a ScalarFunction from a constant scalar value. - * @param[in] x Constant scalar value - */ - ScalarFunction(const ScalarType& x) - : m_x(x) + ScalarFunction(const ScalarFunctionBase& nested) + : m_nested(nested.copy()) {} ScalarFunction(const ScalarFunction& other) : Parent(other), - m_x(other.m_x) + m_nested(other.m_nested->copy()) {} ScalarFunction(ScalarFunction&& other) : Parent(std::move(other)), - m_x(other.m_x) + m_nested(std::move(other.m_nested)) {} inline constexpr - ScalarFunction& traceOf(Geometry::Attribute) - { - return *this; - } - - inline - constexpr - ScalarType getValue() const + auto getValue(const Geometry::Point& v) const { - return m_x; + return m_nested->getValue(v); } inline constexpr - ScalarType getValue(const Geometry::Point&) const + ScalarFunction& traceOf(Geometry::Attribute attrs) { - return m_x; + m_nested->traceOf(attrs); + return *this; } inline ScalarFunction* copy() const noexcept override @@ -274,13 +279,15 @@ namespace Rodin::Variational } private: - const ScalarType m_x; + std::unique_ptr> m_nested; }; /** * @brief CTAD for ScalarFunction. */ - ScalarFunction(const Complex&) -> ScalarFunction; + template + ScalarFunction(const ScalarFunctionBase&) + -> ScalarFunction>; /** * @ingroup ScalarFunctionSpecializations diff --git a/src/Rodin/Variational/Sine.h b/src/Rodin/Variational/Sine.h index b840d48f1..25a1e418d 100644 --- a/src/Rodin/Variational/Sine.h +++ b/src/Rodin/Variational/Sine.h @@ -30,41 +30,51 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; + using ScalarType = Real; + using Parent = RealFunctionBase>>; Sin(const OperandType& v) - : m_v(v.copy()) + : m_operand(v.copy()) {} Sin(const Sin& other) : Parent(other), - m_v(other.m_v->copy()) + m_operand(other.m_operand->copy()) {} Sin(Sin&& other) : Parent(std::move(other)), - m_v(std::move(other.m_v)) + m_operand(std::move(other.m_operand)) {} inline constexpr - Sin& traceOf(Geometry::Attribute attrs) + Sin& traceOf(Geometry::Attribute attr) + { + m_operand->traceOf(attr); + return *this; + } + + inline + constexpr + Sin& traceOf(const FlatSet& attrs) { - m_v.traceOf(attrs); + m_operand->traceOf(attrs); return *this; } inline - auto getValue(const Geometry::Point& p) const + ScalarType getValue(const Geometry::Point& p) const { - return std::sin(static_cast(getOperand().getValue(p))); + return Math::sin(getOperand().getValue(p)); } inline const OperandType& getOperand() const { - assert(m_v); - return *m_v; + assert(m_operand); + return *m_operand; } inline Sin* copy() const noexcept override @@ -73,7 +83,7 @@ namespace Rodin::Variational } private: - std::unique_ptr m_v; + std::unique_ptr m_operand; }; template diff --git a/src/Rodin/Variational/Sqrt.h b/src/Rodin/Variational/Sqrt.h index 49f1ab3ad..5902ff074 100644 --- a/src/Rodin/Variational/Sqrt.h +++ b/src/Rodin/Variational/Sqrt.h @@ -55,9 +55,9 @@ namespace Rodin::Variational } inline - auto getValue(const Geometry::Point& p) const + Real getValue(const Geometry::Point& p) const { - return std::sqrt(static_cast(getOperand().getValue(p))); + return std::sqrt(getOperand().getValue(p)); } inline diff --git a/src/Rodin/Variational/Tangent.h b/src/Rodin/Variational/Tangent.h index d9cfcc23b..cc9a100e1 100644 --- a/src/Rodin/Variational/Tangent.h +++ b/src/Rodin/Variational/Tangent.h @@ -30,37 +30,54 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; + using ScalarType = Real; + using Parent = RealFunctionBase>>; constexpr Tan(const OperandType& v) - : m_v(v) + : m_operand(v.copy()) {} constexpr Tan(const Tan& other) : Parent(other), - m_v(other.m_v) + m_operand(other.m_v->copy()) {} constexpr Tan(Tan&& other) : Parent(std::move(other)), - m_v(std::move(other.m_v)) + m_operand(std::move(other.m_v)) {} inline constexpr - Tan& traceOf(Geometry::Attribute attrs) + Tan& traceOf(Geometry::Attribute attr) { - m_v.traceOf(attrs); + m_operand->traceOf(attr); return *this; } inline - Real getValue(const Geometry::Point& p) const + constexpr + Tan& traceOf(const FlatSet& attrs) + { + m_operand->traceOf(attrs); + return *this; + } + + inline + ScalarType getValue(const Geometry::Point& p) const + { + return Math::tan(Real(getOperand().getValue(p))); + } + + inline + const OperandType& getOperand() const { - return std::tan(Real(m_v.getValue(p))); + assert(m_operand); + return *m_operand; } inline @@ -71,7 +88,7 @@ namespace Rodin::Variational } private: - OperandType m_v; + std::unique_ptr m_operand; }; template diff --git a/src/Rodin/Variational/Traits.h b/src/Rodin/Variational/Traits.h index a4cbfc7e2..efedd7409 100644 --- a/src/Rodin/Variational/Traits.h +++ b/src/Rodin/Variational/Traits.h @@ -27,54 +27,47 @@ namespace Rodin::FormLanguage { namespace Internal { - template + template struct RangeOfSAT {}; template <> - struct RangeOfSAT + struct RangeOfSAT { using Type = Boolean; Variational::RangeType Value = Variational::RangeType::Boolean; }; template <> - struct RangeOfSAT - { - using Type = Real; - Variational::RangeType Value = Variational::RangeType::Real; - }; - - template <> - struct RangeOfSAT + struct RangeOfSAT { using Type = Integer; Variational::RangeType Value = Variational::RangeType::Integer; }; template <> - struct RangeOfSAT + struct RangeOfSAT { using Type = Real; Variational::RangeType Value = Variational::RangeType::Real; }; template <> - struct RangeOfSAT + struct RangeOfSAT { - using Type = Real; - Variational::RangeType Value = Variational::RangeType::Real; + using Type = Complex; + Variational::RangeType Value = Variational::RangeType::Complex; }; template <> - struct RangeOfSAT + struct RangeOfSAT { using Type = Math::Vector; Variational::RangeType Value = Variational::RangeType::Vector; }; template <> - struct RangeOfSAT + struct RangeOfSAT { using Type = Math::Matrix; Variational::RangeType Value = Variational::RangeType::Matrix; @@ -105,18 +98,16 @@ namespace Rodin::FormLanguage static constexpr const bool Value = std::is_same_v; }; - template - struct IsMatrixRange + template + struct IsRealRange { - static constexpr const bool Value = false; + static constexpr const bool Value = std::is_same_v; }; template - struct IsMatrixRange, std::void_t> + struct IsComplexRange { - static constexpr const bool Value = - (T::ColsAtCompileTime == Eigen::Dynamic) || - (T::ColsAtCompileTime > 1); + static constexpr const bool Value = std::is_same_v; }; template @@ -131,13 +122,18 @@ namespace Rodin::FormLanguage static constexpr const bool Value = T::ColsAtCompileTime == 1; }; + template + struct IsMatrixRange + { + static constexpr const bool Value = false; + }; + template - struct IsRealRange + struct IsMatrixRange, std::void_t> { static constexpr const bool Value = - std::is_convertible_v && - !IsVectorRange::Value && - !IsMatrixRange::Value; + (T::ColsAtCompileTime == Eigen::Dynamic) || + (T::ColsAtCompileTime > 1); }; template @@ -168,6 +164,7 @@ namespace Rodin::FormLanguage IsBooleanRange::Value, IsIntegerRange::Value, IsRealRange::Value, + IsComplexRange::Value, IsVectorRange::Value, IsMatrixRange::Value> ::Type; @@ -182,6 +179,7 @@ namespace Rodin::FormLanguage IsBooleanRange::Value, IsIntegerRange::Value, IsRealRange::Value, + IsComplexRange::Value, IsVectorRange::Value, IsMatrixRange::Value> ::Type; diff --git a/src/Rodin/Variational/VectorFunction.h b/src/Rodin/Variational/VectorFunction.h index 57e934b9c..236995bfe 100644 --- a/src/Rodin/Variational/VectorFunction.h +++ b/src/Rodin/Variational/VectorFunction.h @@ -46,8 +46,6 @@ namespace Rodin::Variational public: using ScalarType = Scalar; - using VectorType = Math::Vector; - using Parent = FunctionBase>; VectorFunctionBase() = default; @@ -151,6 +149,7 @@ namespace Rodin::Variational return static_cast(*this).getValue(p); } + template inline constexpr void getValue(VectorType& res, const Geometry::Point& p) const From 721fc1e96dedfddcdce11e96a5f7c5060b369e7e Mon Sep 17 00:00:00 2001 From: Carlos Brito Date: Wed, 10 Jul 2024 11:24:48 +0200 Subject: [PATCH 4/4] Work towards Complex number support --- src/Rodin/Variational/Cosine.h | 4 +- src/Rodin/Variational/Div.h | 6 +- src/Rodin/Variational/Division.h | 21 +--- src/Rodin/Variational/Dot.h | 18 +-- src/Rodin/Variational/Function.h | 60 +++------ src/Rodin/Variational/Grad.h | 11 +- src/Rodin/Variational/Im.h | 4 +- src/Rodin/Variational/MatrixFunction.h | 23 +--- src/Rodin/Variational/Max.h | 8 -- src/Rodin/Variational/Min.h | 4 - src/Rodin/Variational/Mult.h | 22 ++-- src/Rodin/Variational/P1/Div.h | 4 +- src/Rodin/Variational/P1/GridFunction.h | 13 +- src/Rodin/Variational/Re.h | 4 +- src/Rodin/Variational/ShapeFunction.h | 143 +++++++++++----------- src/Rodin/Variational/Sine.h | 4 +- src/Rodin/Variational/Sqrt.h | 2 +- src/Rodin/Variational/Sum.h | 35 +++--- src/Rodin/Variational/Tangent.h | 11 +- src/Rodin/Variational/Trace.h | 21 +--- src/Rodin/Variational/Traits.h | 156 ++++++------------------ src/Rodin/Variational/Transpose.h | 23 +--- src/Rodin/Variational/UnaryMinus.h | 22 +--- src/Rodin/Variational/Zero.h | 2 +- 24 files changed, 209 insertions(+), 412 deletions(-) diff --git a/src/Rodin/Variational/Cosine.h b/src/Rodin/Variational/Cosine.h index 239a033d6..c49b3ab25 100644 --- a/src/Rodin/Variational/Cosine.h +++ b/src/Rodin/Variational/Cosine.h @@ -30,8 +30,6 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; - using ScalarType = Real; - using Parent = RealFunctionBase>>; Cos(const OperandType& v) @@ -65,7 +63,7 @@ namespace Rodin::Variational } inline - ScalarType getValue(const Geometry::Point& p) const + Real getValue(const Geometry::Point& p) const { return Math::cos(getOperand().getValue(p)); } diff --git a/src/Rodin/Variational/Div.h b/src/Rodin/Variational/Div.h index 23731fed8..f3f40ac17 100644 --- a/src/Rodin/Variational/Div.h +++ b/src/Rodin/Variational/Div.h @@ -29,7 +29,7 @@ namespace Rodin::Variational */ template class DivBase, Derived> - : public RealFunctionBase, Derived>> + : public ScalarFunctionBase::ScalarType, DivBase, Derived>> { public: using FESType = FES; @@ -39,7 +39,7 @@ namespace Rodin::Variational using OperandType = GridFunction; /// Parent class - using Parent = RealFunctionBase>; + using Parent = ScalarFunctionBase>; /** * @brief Constructs the Div of a @f$ \mathbb{P}_1 @f$ function @f$ u @@ -108,7 +108,7 @@ namespace Rodin::Variational */ inline constexpr - void interpolate(Real& out, const Geometry::Point& p) const + void interpolate(ScalarType& out, const Geometry::Point& p) const { static_cast(*this).interpolate(out, p); } diff --git a/src/Rodin/Variational/Division.h b/src/Rodin/Variational/Division.h index 41dd6d73d..bf22071f3 100644 --- a/src/Rodin/Variational/Division.h +++ b/src/Rodin/Variational/Division.h @@ -38,12 +38,11 @@ namespace Rodin::Variational using Parent = FunctionBase, FunctionBase>>; - static_assert(std::is_same_v); - Division(const FunctionBase& lhs, const FunctionBase& rhs) : m_lhs(lhs.copy()), m_rhs(rhs.copy()) {} + Division(const Division& other) : Parent(other), m_lhs(other.m_lhs->copy()), m_rhs(other.m_rhs->copy()) @@ -54,14 +53,12 @@ namespace Rodin::Variational m_lhs(std::move(other.m_lhs)), m_rhs(std::move(other.m_rhs)) {} - inline constexpr RangeShape getRangeShape() const { return getLHS().getRangeShape(); } - inline constexpr Division& traceOf(Geometry::Attribute attr) { @@ -70,7 +67,6 @@ namespace Rodin::Variational return *this; } - inline constexpr const LHSType& getLHS() const { @@ -78,7 +74,6 @@ namespace Rodin::Variational return *m_lhs; } - inline constexpr const RHSType& getRHS() const { @@ -86,28 +81,20 @@ namespace Rodin::Variational return *m_rhs; } - inline + constexpr auto getValue(const Geometry::Point& p) const { return this->object(getLHS().getValue(p)) / this->object(getRHS().getValue(p)); } - inline + template constexpr - void getValue(Math::Vector& res, const Geometry::Point& p) const + void getValue(T& res, const Geometry::Point& p) const { getLHS().getDerived().getValue(res, p); res /= getRHS().getValue(p); } - inline - constexpr - void getValue(Math::Matrix& res, const Geometry::Point& p) const - { - getLHS().getValue(res, p); - res /= getRHS().getValue(p); - } - inline Division* copy() const noexcept final override { diff --git a/src/Rodin/Variational/Dot.h b/src/Rodin/Variational/Dot.h index a1db812c6..affaf07ea 100644 --- a/src/Rodin/Variational/Dot.h +++ b/src/Rodin/Variational/Dot.h @@ -39,8 +39,7 @@ namespace Rodin::FormLanguage using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = - decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using RangeType = ScalarType; }; @@ -66,8 +65,7 @@ namespace Rodin::FormLanguage using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = - decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using RangeType = ScalarType; }; @@ -90,8 +88,7 @@ namespace Rodin::FormLanguage using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = - decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using RangeType = ScalarType; }; @@ -125,8 +122,7 @@ namespace Rodin::Variational using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = - decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using RangeType = ScalarType; @@ -245,8 +241,7 @@ namespace Rodin::Variational using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = - decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using RangeType = ScalarType; @@ -397,8 +392,7 @@ namespace Rodin::Variational using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = - decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using RangeType = ScalarType; diff --git a/src/Rodin/Variational/Function.h b/src/Rodin/Variational/Function.h index bf59e89ef..a1ddd8f13 100644 --- a/src/Rodin/Variational/Function.h +++ b/src/Rodin/Variational/Function.h @@ -109,20 +109,17 @@ namespace Rodin::Variational virtual ~FunctionBase() = default; - inline FunctionBase& operator=(FunctionBase&& other) { m_traceDomain = std::move(other.m_traceDomain); return *this; } - inline Derived& getDerived() { return static_cast(*this); } - inline const Derived& getDerived() const { return static_cast(*this); @@ -137,21 +134,18 @@ namespace Rodin::Variational * boundaries. If the trace domain is empty, then this has the * semantic value that it has not been specified yet. */ - inline constexpr const TraceDomain& getTraceDomain() const { return m_traceDomain; } - inline constexpr Transpose T() const { return Transpose(*this); } - inline constexpr RangeShape getRangeShape() const { @@ -162,7 +156,6 @@ namespace Rodin::Variational * @brief Evaluates the function on a Point belonging to the mesh. * @note CRTP function to be overriden in Derived class. */ - inline constexpr auto getValue(const Geometry::Point& p) const { @@ -170,14 +163,12 @@ namespace Rodin::Variational } template - inline constexpr - auto getValue(T& res, const Geometry::Point& p) const + void getValue(T& res, const Geometry::Point& p) const { - if constexpr ( - Internal::HasGetValueMethod::Value) + if constexpr (Internal::HasGetValueMethod::Value) { - return static_cast(*this).getValue(res, p); + static_cast(*this).getValue(res, p); } else { @@ -185,42 +176,34 @@ namespace Rodin::Variational } } - inline - auto coeff(size_t i, size_t j) const - { - return Component(*this, i, j); - } - - inline - auto coeff(size_t i) const - { - return Component(*this, i); - } - - // inline - // void operator()(Math::Vector& res, const Geometry::Point& p) const - // { - // return getValue(res, p); - // } - - // inline - // void operator()(Math::Matrix& res, const Geometry::Point& p) const - // { - // return getValue(res, p); - // } - /** * @brief Evaluates the function on a Point belonging to the mesh. * * This calls the function get getValue(const Geometry::Point&). */ - inline constexpr auto operator()(const Geometry::Point& p) const { return getValue(p); } + template + constexpr + void operator()(T& res, const Geometry::Point& p) const + { + getValue(res, p); + } + + auto coeff(size_t i, size_t j) const + { + return Component(*this, i, j); + } + + auto coeff(size_t i) const + { + return Component(*this, i); + } + /** * @brief Sets an attribute which will be interpreted as the domain to * trace. @@ -230,7 +213,6 @@ namespace Rodin::Variational * * @returns Reference to self (for method chaining) */ - inline constexpr Derived& traceOf(Geometry::Attribute attr) { @@ -238,14 +220,12 @@ namespace Rodin::Variational } template - inline constexpr Derived& traceOf(A1 a1, A2 a2, As ... as) { return traceOf(FlatSet{ a1, a2, as... }); } - inline constexpr Derived& traceOf(const FlatSet& attr) { diff --git a/src/Rodin/Variational/Grad.h b/src/Rodin/Variational/Grad.h index d3218f329..4aff07cb8 100644 --- a/src/Rodin/Variational/Grad.h +++ b/src/Rodin/Variational/Grad.h @@ -45,7 +45,7 @@ namespace Rodin::Variational using ScalarType = typename FormLanguage::Traits::ScalarType; - using VectorType = Math::SpatialVector; + using SpatialVectorType = Math::SpatialVector; using OperandType = GridFunction; @@ -86,17 +86,16 @@ namespace Rodin::Variational } inline - Math::SpatialVector getValue(const Geometry::Point& p) const + SpatialVectorType getValue(const Geometry::Point& p) const { - Math::SpatialVector out; + SpatialVectorType out; getValue(out, p); return out; } inline - void getValue(Math::SpatialVector& out, const Geometry::Point& p) const + void getValue(SpatialVectorType& out, const Geometry::Point& p) const { - out.setConstant(NAN); const auto& polytope = p.getPolytope(); const auto& polytopeMesh = polytope.getMesh(); const auto& gf = getOperand(); @@ -127,7 +126,7 @@ namespace Rodin::Variational */ inline constexpr - void interpolate(Math::SpatialVector& out, const Geometry::Point& p) const + void interpolate(SpatialVectorType& out, const Geometry::Point& p) const { static_cast(*this).interpolate(out, p); } diff --git a/src/Rodin/Variational/Im.h b/src/Rodin/Variational/Im.h index a5b02a962..039419942 100644 --- a/src/Rodin/Variational/Im.h +++ b/src/Rodin/Variational/Im.h @@ -15,8 +15,6 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; - using ScalarType = Real; - using Parent = RealFunctionBase>>; Im(const OperandType& f) @@ -41,7 +39,7 @@ namespace Rodin::Variational inline constexpr - ScalarType getValue(const Geometry::Point& p) const + Real getValue(const Geometry::Point& p) const { return getOperand().getValue(p).imag(); } diff --git a/src/Rodin/Variational/MatrixFunction.h b/src/Rodin/Variational/MatrixFunction.h index e46c21e71..4047c2d0e 100644 --- a/src/Rodin/Variational/MatrixFunction.h +++ b/src/Rodin/Variational/MatrixFunction.h @@ -54,13 +54,11 @@ namespace Rodin::Variational virtual ~MatrixFunctionBase() = default; - inline const Derived& getDerived() const { return static_cast(*this); } - inline constexpr auto getValue(const Geometry::Point& p) const { @@ -68,7 +66,6 @@ namespace Rodin::Variational } template - inline constexpr void getValue(MatrixType& res, const Geometry::Point& p) const { @@ -82,7 +79,6 @@ namespace Rodin::Variational } } - inline constexpr RangeShape getRangeShape() const { @@ -93,7 +89,6 @@ namespace Rodin::Variational * @brief Gets the number of rows in the matrix * @returns Number of rows */ - inline constexpr size_t getRows() const { @@ -104,14 +99,12 @@ namespace Rodin::Variational * @brief Gets the number of columns in the matrix * @returns Number of columns */ - inline constexpr size_t getColumns() const { return static_cast(*this).getColumns(); } - inline constexpr MatrixFunctionBase& traceOf(Geometry::Attribute attr) { @@ -119,7 +112,6 @@ namespace Rodin::Variational return *this; } - inline constexpr MatrixFunctionBase& traceOf(const FlatSet& attrs) { @@ -127,7 +119,7 @@ namespace Rodin::Variational return *this; } - virtual inline MatrixFunctionBase* copy() const noexcept override + virtual MatrixFunctionBase* copy() const noexcept override { return static_cast(*this).copy(); } @@ -161,35 +153,30 @@ namespace Rodin::Variational m_matrix(std::move(other.m_matrix)) {} - inline constexpr const MatrixType& getValue(const Geometry::Point&) const { return m_matrix.get(); } - inline constexpr void getValue(MatrixType& res, const Geometry::Point&) const { res = m_matrix.get(); } - inline constexpr MatrixFunction& traceOf(Geometry::Attribute) { return *this; } - inline constexpr MatrixFunction& traceOf(const FlatSet& attr) { return *this; } - inline constexpr size_t getRows() const { @@ -200,14 +187,13 @@ namespace Rodin::Variational * @brief Gets the number of columns in the matrix * @returns Number of columns */ - inline constexpr size_t getColumns() const { return m_matrix.get().cols(); } - inline MatrixFunction* copy() const noexcept override + MatrixFunction* copy() const noexcept override { return new MatrixFunction(*this); } @@ -216,8 +202,9 @@ namespace Rodin::Variational std::reference_wrapper m_matrix; }; - MatrixFunction(std::reference_wrapper>) - -> MatrixFunction>; + template + MatrixFunction(std::reference_wrapper>) + -> MatrixFunction>; } #endif diff --git a/src/Rodin/Variational/Max.h b/src/Rodin/Variational/Max.h index 5cffe0d44..787256141 100644 --- a/src/Rodin/Variational/Max.h +++ b/src/Rodin/Variational/Max.h @@ -49,7 +49,6 @@ namespace Rodin::Variational m_lhs(std::move(other.m_lhs)), m_rhs(std::move(other.m_rhs)) {} - inline constexpr Max& traceOf(Geometry::Attribute attrs) { @@ -58,7 +57,6 @@ namespace Rodin::Variational return *this; } - inline constexpr auto getValue(const Geometry::Point& p) const { @@ -70,14 +68,12 @@ namespace Rodin::Variational return lhs; } - inline const auto& getLHS() const { assert(m_lhs); return *m_lhs; } - inline const auto& getRHS() const { assert(m_rhs); @@ -129,7 +125,6 @@ namespace Rodin::Variational m_lhs(std::move(other.m_lhs)), m_rhs(std::move(other.m_rhs)) {} - inline constexpr Max& traceOf(Geometry::Attribute attrs) { @@ -137,7 +132,6 @@ namespace Rodin::Variational return *this; } - inline constexpr Real getValue(const Geometry::Point& p) const { @@ -149,14 +143,12 @@ namespace Rodin::Variational return lhs; } - inline const auto& getLHS() const { assert(m_lhs); return *m_lhs; } - inline const auto& getRHS() const { return m_rhs; diff --git a/src/Rodin/Variational/Min.h b/src/Rodin/Variational/Min.h index 7176f7579..ce1340b64 100644 --- a/src/Rodin/Variational/Min.h +++ b/src/Rodin/Variational/Min.h @@ -129,7 +129,6 @@ namespace Rodin::Variational m_lhs(std::move(other.m_lhs)), m_rhs(std::move(other.m_rhs)) {} - inline constexpr Min& traceOf(Geometry::Attribute attrs) { @@ -137,7 +136,6 @@ namespace Rodin::Variational return *this; } - inline constexpr ScalarType getValue(const Geometry::Point& p) const { @@ -149,14 +147,12 @@ namespace Rodin::Variational return rhs; } - inline const auto& getLHS() const { assert(m_lhs); return *m_lhs; } - inline const auto& getRHS() const { return m_rhs; diff --git a/src/Rodin/Variational/Mult.h b/src/Rodin/Variational/Mult.h index fa384f521..cacfcc9d9 100644 --- a/src/Rodin/Variational/Mult.h +++ b/src/Rodin/Variational/Mult.h @@ -505,7 +505,7 @@ namespace Rodin::Variational using RHSScalarType = typename FormLanguage::Traits::ScalarType; - using ScalarType = decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using Parent = ShapeFunctionBase, FES, SpaceType>; @@ -642,16 +642,16 @@ namespace Rodin::Variational return Mult(lhs, RealFunction(rhs)); } - template - class Mult> - : public LocalBilinearFormIntegratorBase() * std::declval())> + template + class Mult> + : public LocalBilinearFormIntegratorBase::Type> { public: using LHSType = LHS; - using RHSType = LocalBilinearFormIntegratorBase; + using RHSType = LocalBilinearFormIntegratorBase; - using ScalarType = decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using Parent = LocalBilinearFormIntegratorBase; @@ -725,16 +725,16 @@ namespace Rodin::Variational return Mult(lhs, rhs); } - template - class Mult> - : public LinearFormIntegratorBase() * std::declval())> + template + class Mult> + : public LinearFormIntegratorBase::Type> { public: using LHSType = LHS; - using RHSType = LinearFormIntegratorBase; + using RHSType = LinearFormIntegratorBase; - using ScalarType = decltype(std::declval() * std::declval()); + using ScalarType = typename FormLanguage::Mult::Type; using Parent = LinearFormIntegratorBase; diff --git a/src/Rodin/Variational/P1/Div.h b/src/Rodin/Variational/P1/Div.h index ed0eb41f7..3a87f8119 100644 --- a/src/Rodin/Variational/P1/Div.h +++ b/src/Rodin/Variational/P1/Div.h @@ -172,10 +172,10 @@ namespace Rodin::Variational */ template class Div, Mesh>, Space>> final - : public ShapeFunctionBase, Mesh>, Space>>> + : public ShapeFunctionBase, Mesh>, Space>>> { public: - using FESType = P1, Mesh>; + using FESType = P1, Mesh>; static constexpr ShapeFunctionSpaceType SpaceType = Space; using OperandType = ShapeFunction; diff --git a/src/Rodin/Variational/P1/GridFunction.h b/src/Rodin/Variational/P1/GridFunction.h index b8f138cc7..b040333d9 100644 --- a/src/Rodin/Variational/P1/GridFunction.h +++ b/src/Rodin/Variational/P1/GridFunction.h @@ -132,10 +132,11 @@ namespace Rodin::Variational const size_t vdim = fes.getVectorDimension(); res.resize(vdim); res.setZero(); + Math::Vector basis; for (Index local = 0; local < fe.getCount(); local++) { - const auto& basis = fe.getBasis(local); - res += getValue({d, i}, local).coeff(local % vdim) * basis(r); + fe.getBasis(local)(basis, r); + res += getValue({d, i}, local).coeff(local % vdim) * basis; } } else @@ -148,12 +149,12 @@ namespace Rodin::Variational { auto& data = this->getData(); auto& weights = this->getWeights().emplace(this->getFiniteElementSpace().getSize()); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(data.rows() == 1); std::copy(data.data(), data.data() + data.size(), weights.data()); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { const auto& fes = this->getFiniteElementSpace(); const size_t vdim = fes.getVectorDimension(); @@ -175,12 +176,12 @@ namespace Rodin::Variational assert(static_cast(weights.size()) == this->getFiniteElementSpace().getSize()); auto& data = this->getData(); const auto& w = this->getWeights().emplace(std::forward(weights)); - if constexpr (std::is_same_v) + if constexpr (std::is_same_v) { assert(data.rows() == 1); std::copy(w.data(), w.data() + w.size(), data.data()); } - else if constexpr (std::is_same_v>) + else if constexpr (std::is_same_v>) { const size_t sz = w.size(); const auto& fes = this->getFiniteElementSpace(); diff --git a/src/Rodin/Variational/Re.h b/src/Rodin/Variational/Re.h index fe382649e..189c839ea 100644 --- a/src/Rodin/Variational/Re.h +++ b/src/Rodin/Variational/Re.h @@ -15,8 +15,6 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; - using ScalarType = Real; - using Parent = RealFunctionBase>>; Re(const OperandType& f) @@ -41,7 +39,7 @@ namespace Rodin::Variational inline constexpr - ScalarType getValue(const Geometry::Point& p) const + Real getValue(const Geometry::Point& p) const { return getOperand().getValue(p).real(); } diff --git a/src/Rodin/Variational/ShapeFunction.h b/src/Rodin/Variational/ShapeFunction.h index 6ffb52787..928249bea 100644 --- a/src/Rodin/Variational/ShapeFunction.h +++ b/src/Rodin/Variational/ShapeFunction.h @@ -60,6 +60,49 @@ namespace Rodin::FormLanguage namespace Rodin::Variational { + namespace Internal + { + template + struct HasGetBasisMethod + { + template().getBasis(std::declval()...))> + static std::true_type Test(int); + + template + static std::false_type Test(...); + + using Type = decltype(Test(0)); + static constexpr bool Value = Type::value; + }; + + template + struct HasGetBasisMethod + { + template ().getBasis(std::declval()...))> + static std::true_type Test(int); + + template + static std::false_type Test(...); + + using Type = decltype(Test(0)); + static constexpr bool Value = Type::value; + }; + + template + struct HasGetBasisMethodR + { + template().getBasis(std::declval()...))> + static auto Test(int) -> + decltype(std::is_same::type, T>::value, std::true_type{}); + + template + static std::false_type Test(...); + + using Type = decltype(Test(0)); + static constexpr bool Value = Type::value; + }; + } + /** * @defgroup ShapeFunctionSpecializations ShapeFunction Template Specializations * @brief Template specializations of the ShapeFunction class. @@ -119,13 +162,11 @@ namespace Rodin::Variational m_fes(std::move(other.m_fes)) {} - inline Derived& getDerived() { return static_cast(*this); } - inline const Derived& getDerived() const { return static_cast(*this); @@ -135,7 +176,6 @@ namespace Rodin::Variational * @brief Indicates whether the shape function is part of a %Trial or %Test * function expression. */ - inline constexpr ShapeFunctionSpaceType getSpaceType() const { @@ -146,64 +186,27 @@ namespace Rodin::Variational * @brief Gets the shape of the range space. * @note CRTP function to be overriden in the Derived class. */ - inline constexpr RangeShape getRangeShape() const { return static_cast(*this).getRangeShape(); } - inline - constexpr - RangeType getRangeType() const - { - using R = typename FormLanguage::Traits>::RangeType; - if constexpr (std::is_same_v) - { - return RangeType::Boolean; - } - else if constexpr (std::is_same_v) - { - return RangeType::Integer; - } - else if constexpr (std::is_same_v) - { - return RangeType::Real; - } - else if constexpr (std::is_same_v>) - { - return RangeType::Vector; - } - else if constexpr (std::is_same_v>) - { - return RangeType::Matrix; - } - else - { - assert(false); - static_assert(Utility::DependentFalse::Value); - } - } - - inline auto x() const { return Component(*this, 0); } - inline auto y() const { return Component(*this, 1); } - inline auto z() const { return Component(*this, 2); } - inline constexpr auto T() const { @@ -214,7 +217,6 @@ namespace Rodin::Variational * @brief Gets the operand in the shape function expression. * @note CRTP function to be overriden in the Derived class. */ - inline constexpr const auto& getLeaf() const { @@ -226,20 +228,17 @@ namespace Rodin::Variational * @param[in] polytope Polytope * @note CRTP function to be overriden in the Derived class. */ - inline constexpr size_t getDOFs(const Geometry::Polytope& polytope) const { return static_cast(*this).getDOFs(polytope); } - inline const Geometry::Point& getPoint() const { return static_cast(*this).getPoint(); } - inline constexpr Derived& setPoint(const Geometry::Point& p) { @@ -252,32 +251,49 @@ namespace Rodin::Variational * @param[in] p Point where the shape function basis will be calculated * @note CRTP function to be overriden in the Derived class. */ - inline constexpr auto getBasis(size_t local) const { return static_cast(*this).getBasis(local); } + template + constexpr + void getBasis(T& basis, size_t local) const + { + if constexpr (Internal::HasGetBasisMethod::Value) + { + static_cast(*this).getBasis(basis, local); + } + else + { + basis = getBasis(local); + } + } + /** * @brief Call operator to get an expression which yields the shape * function basis at the given point. * - * Synonym to getBasis(). + * Synonym to getBasis(size_t). */ - template - inline constexpr - auto operator()(Args&&... args) const + auto operator()(size_t local) const + { + return getBasis(local); + } + + template + constexpr + void operator()(T& res, size_t local) const { - return this->getBasis(std::forward(args)...); + getBasis(res, local); } /** * @brief Gets the finite element space to which the shape function * belongs to. */ - inline constexpr const FES& getFiniteElementSpace() const { @@ -305,6 +321,8 @@ namespace Rodin::Variational using FESType = FES; static constexpr ShapeFunctionSpaceType SpaceType = Space; + using ScalarType = typename FormLanguage::Traits::ScalarType; + using RangeType = typename FormLanguage::Traits::RangeType; using Parent = ShapeFunctionBase, FESType, SpaceType>; @@ -326,7 +344,6 @@ namespace Rodin::Variational : Parent(std::move(other)) {} - inline constexpr auto& emplace() { @@ -334,14 +351,12 @@ namespace Rodin::Variational return *this; } - inline constexpr RangeShape getRangeShape() const { return { this->getFiniteElementSpace().getVectorDimension(), 1 }; } - inline constexpr GridFunction& getSolution() { @@ -349,7 +364,6 @@ namespace Rodin::Variational return m_gf.value(); } - inline constexpr const GridFunction& getSolution() const { @@ -357,7 +371,6 @@ namespace Rodin::Variational return m_gf.value(); } - inline constexpr size_t getDOFs(const Geometry::Polytope& polytope) const { @@ -366,21 +379,18 @@ namespace Rodin::Variational return this->getFiniteElementSpace().getFiniteElement(d, i).getCount(); } - inline const Geometry::Point& getPoint() const { assert(m_p.has_value()); return m_p.value().get(); } - inline ShapeFunction& setPoint(const Geometry::Point& p) { m_p = p; return *this; } - inline constexpr auto getBasis(size_t local) const { @@ -389,22 +399,9 @@ namespace Rodin::Variational const Index i = p.getPolytope().getIndex(); const auto& fes = this->getFiniteElementSpace(); const auto& fe = fes.getFiniteElement(d, i); - if constexpr (std::is_same_v) - { - return fes.getInverseMapping({ d, i }, fe.getBasis(local))(p); - } - else if constexpr (std::is_same_v>) - { - return this->object(fes.getInverseMapping({ d, i }, fe.getBasis(local))(p)); - } - else - { - assert(false); - return void(); - } + return this->object(fes.getInverseMapping({ d, i }, fe.getBasis(local))(p)); } - inline constexpr const auto& getLeaf() const { diff --git a/src/Rodin/Variational/Sine.h b/src/Rodin/Variational/Sine.h index 25a1e418d..3dc0ae94d 100644 --- a/src/Rodin/Variational/Sine.h +++ b/src/Rodin/Variational/Sine.h @@ -30,8 +30,6 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; - using ScalarType = Real; - using Parent = RealFunctionBase>>; Sin(const OperandType& v) @@ -65,7 +63,7 @@ namespace Rodin::Variational } inline - ScalarType getValue(const Geometry::Point& p) const + Real getValue(const Geometry::Point& p) const { return Math::sin(getOperand().getValue(p)); } diff --git a/src/Rodin/Variational/Sqrt.h b/src/Rodin/Variational/Sqrt.h index 5902ff074..51f60718c 100644 --- a/src/Rodin/Variational/Sqrt.h +++ b/src/Rodin/Variational/Sqrt.h @@ -57,7 +57,7 @@ namespace Rodin::Variational inline Real getValue(const Geometry::Point& p) const { - return std::sqrt(getOperand().getValue(p)); + return Math::sqrt(getOperand().getValue(p)); } inline diff --git a/src/Rodin/Variational/Sum.h b/src/Rodin/Variational/Sum.h index 4403df37f..b19928a2b 100644 --- a/src/Rodin/Variational/Sum.h +++ b/src/Rodin/Variational/Sum.h @@ -56,7 +56,8 @@ namespace Rodin::FormLanguage using RHSType = Variational::LinearFormIntegratorBase; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = + typename FormLanguage::Sum::Type; }; } @@ -361,7 +362,7 @@ namespace Rodin::Variational template class Sum, LinearFormIntegratorBase> : public FormLanguage::List< - LinearFormIntegratorBase() + std::declval())>> + LinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -372,7 +373,7 @@ namespace Rodin::Variational using RHSType = LinearFormIntegratorBase; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -409,7 +410,7 @@ namespace Rodin::Variational template class Sum, FormLanguage::List>> : public FormLanguage::List< - LinearFormIntegratorBase() + std::declval())>> + LinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -420,7 +421,7 @@ namespace Rodin::Variational using RHSType = FormLanguage::List>; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -458,7 +459,7 @@ namespace Rodin::Variational class Sum< FormLanguage::List>, LinearFormIntegratorBase> : public FormLanguage::List< - LinearFormIntegratorBase() + std::declval())>> + LinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -469,7 +470,7 @@ namespace Rodin::Variational using RHSType = LinearFormIntegratorBase; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -509,7 +510,7 @@ namespace Rodin::Variational FormLanguage::List>, FormLanguage::List>> : public FormLanguage::List< - LinearFormIntegratorBase() + std::declval())>> + LinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -520,7 +521,7 @@ namespace Rodin::Variational using RHSType = FormLanguage::List>; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -560,7 +561,7 @@ namespace Rodin::Variational template class Sum, LocalBilinearFormIntegratorBase> : public FormLanguage::List< - LocalBilinearFormIntegratorBase() + std::declval())>> + LocalBilinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -571,7 +572,7 @@ namespace Rodin::Variational using RHSType = LocalBilinearFormIntegratorBase; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -609,7 +610,7 @@ namespace Rodin::Variational template class Sum, FormLanguage::List>> : public FormLanguage::List< - LocalBilinearFormIntegratorBase() + std::declval())>> + LocalBilinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -620,7 +621,7 @@ namespace Rodin::Variational using RHSType = FormLanguage::List>; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -661,7 +662,7 @@ namespace Rodin::Variational class Sum< FormLanguage::List>, LocalBilinearFormIntegratorBase> : public FormLanguage::List< - LocalBilinearFormIntegratorBase() + std::declval())>> + LocalBilinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -672,7 +673,7 @@ namespace Rodin::Variational using RHSType = LocalBilinearFormIntegratorBase; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; @@ -714,7 +715,7 @@ namespace Rodin::Variational FormLanguage::List>, FormLanguage::List>> : public FormLanguage::List< - LocalBilinearFormIntegratorBase() + std::declval())>> + LocalBilinearFormIntegratorBase::Type>> { public: using LHSScalarType = LHSNumber; @@ -725,7 +726,7 @@ namespace Rodin::Variational using RHSType = FormLanguage::List>; - using ScalarType = decltype(std::declval() + std::declval()); + using ScalarType = typename FormLanguage::Sum::Type; using Parent = FormLanguage::List>; diff --git a/src/Rodin/Variational/Tangent.h b/src/Rodin/Variational/Tangent.h index cc9a100e1..23094c63b 100644 --- a/src/Rodin/Variational/Tangent.h +++ b/src/Rodin/Variational/Tangent.h @@ -30,8 +30,6 @@ namespace Rodin::Variational public: using OperandType = FunctionBase; - using ScalarType = Real; - using Parent = RealFunctionBase>>; constexpr @@ -51,7 +49,6 @@ namespace Rodin::Variational m_operand(std::move(other.m_v)) {} - inline constexpr Tan& traceOf(Geometry::Attribute attr) { @@ -59,7 +56,6 @@ namespace Rodin::Variational return *this; } - inline constexpr Tan& traceOf(const FlatSet& attrs) { @@ -67,20 +63,17 @@ namespace Rodin::Variational return *this; } - inline - ScalarType getValue(const Geometry::Point& p) const + Real getValue(const Geometry::Point& p) const { - return Math::tan(Real(getOperand().getValue(p))); + return Math::tan(getOperand().getValue(p)); } - inline const OperandType& getOperand() const { assert(m_operand); return *m_operand; } - inline Tan* copy() const noexcept override { diff --git a/src/Rodin/Variational/Trace.h b/src/Rodin/Variational/Trace.h index 84a57a3fa..61946f405 100644 --- a/src/Rodin/Variational/Trace.h +++ b/src/Rodin/Variational/Trace.h @@ -9,7 +9,6 @@ #include "ForwardDecls.h" -#include "RealFunction.h" #include "ShapeFunction.h" namespace Rodin::FormLanguage @@ -36,12 +35,12 @@ namespace Rodin::Variational */ template class Trace> final - : public RealFunctionBase>> + : public FunctionBase>> { public: using OperandType = FunctionBase; - using Parent = RealFunctionBase>; + using Parent = FunctionBase>; /** * @brief Constructs the Trace of the given matrix @@ -64,16 +63,12 @@ namespace Rodin::Variational m_operand(std::move(other.m_operand)) {} - inline constexpr auto getValue(const Geometry::Point& p) const { - using OperandRange = typename FormLanguage::Traits::RangeType; - static_assert(std::is_same_v>); return getOperand().getValue(p).trace(); } - inline constexpr const OperandType& getOperand() const { @@ -81,14 +76,13 @@ namespace Rodin::Variational return *m_operand; } - inline Trace& traceOf(Geometry::Attribute attrs) { m_operand.traceOf(attrs); return *this; } - inline Trace* copy() const noexcept override + Trace* copy() const noexcept override { return new Trace(*this); } @@ -130,48 +124,41 @@ namespace Rodin::Variational m_operand(std::move(other.m_operand)) {} - inline constexpr const OperandType& getOperand() const { return *m_operand; } - inline constexpr const auto& getLeaf() const { return getOperand().getLeaf(); } - inline constexpr RangeShape getRangeShape() const { return { 1, 1 }; } - inline constexpr size_t getDOFs(const Geometry::Polytope& element) const { return getOperand().getDOFs(element); } - inline const Geometry::Point& getPoint() const { return m_operand->getPoint(); } - inline Trace& setPoint(const Geometry::Point& p) { m_operand->setPoint(p); return *this; } - inline constexpr auto getBasis(size_t local) const { @@ -183,7 +170,7 @@ namespace Rodin::Variational return getOperand().getFiniteElementSpace(); } - inline Trace* copy() const noexcept override + Trace* copy() const noexcept override { return new Trace(*this); } diff --git a/src/Rodin/Variational/Traits.h b/src/Rodin/Variational/Traits.h index efedd7409..d4e10205f 100644 --- a/src/Rodin/Variational/Traits.h +++ b/src/Rodin/Variational/Traits.h @@ -25,164 +25,86 @@ namespace Rodin::FormLanguage { - namespace Internal - { - template - struct RangeOfSAT - {}; - - template <> - struct RangeOfSAT - { - using Type = Boolean; - Variational::RangeType Value = Variational::RangeType::Boolean; - }; - - template <> - struct RangeOfSAT - { - using Type = Integer; - Variational::RangeType Value = Variational::RangeType::Integer; - }; - - template <> - struct RangeOfSAT - { - using Type = Real; - Variational::RangeType Value = Variational::RangeType::Real; - }; - - template <> - struct RangeOfSAT - { - using Type = Complex; - Variational::RangeType Value = Variational::RangeType::Complex; - }; - - template <> - struct RangeOfSAT - { - using Type = Math::Vector; - Variational::RangeType Value = Variational::RangeType::Vector; - }; - - template <> - struct RangeOfSAT - { - using Type = Math::Matrix; - Variational::RangeType Value = Variational::RangeType::Matrix; - }; - } - - template - struct IsVectorAtCompileTime - { - static constexpr const bool Value = false; - }; - template - struct IsVectorAtCompileTime> - { - static constexpr const bool Value = T::IsVectorAtCompileTime; - }; + struct ResultOf; - template - struct IsBooleanRange + template + struct ResultOf> { - static constexpr const bool Value = std::is_same_v; + using Type = + std::invoke_result_t, const Geometry::Point&>; }; - template - struct IsIntegerRange + template + struct ResultOf> { - static constexpr const bool Value = std::is_same_v; + using Type = std::invoke_result_t, size_t>; }; template - struct IsRealRange - { - static constexpr const bool Value = std::is_same_v; - }; + struct RangeOf; - template - struct IsComplexRange + template <> + struct RangeOf { - static constexpr const bool Value = std::is_same_v; + using Type = Boolean; }; - template - struct IsVectorRange + template <> + struct RangeOf { - static constexpr const bool Value = false; + using Type = Integer; }; - template - struct IsVectorRange, std::void_t> + template <> + struct RangeOf { - static constexpr const bool Value = T::ColsAtCompileTime == 1; + using Type = Real; }; - template - struct IsMatrixRange + template <> + struct RangeOf { - static constexpr const bool Value = false; + using Type = Complex; }; - template - struct IsMatrixRange, std::void_t> + template + struct RangeOf> { - static constexpr const bool Value = - (T::ColsAtCompileTime == Eigen::Dynamic) || - (T::ColsAtCompileTime > 1); + using Type = Math::Vector; }; - template - struct ResultOf; - - template - struct ResultOf> + template + struct RangeOf> { - using Type = - std::invoke_result_t, const Geometry::Point&>; + using Type = Math::Matrix; }; - template - struct ResultOf> + template + struct RangeOf { - using Type = std::invoke_result_t, size_t>; + using Type = + std::conditional_t< + MatrixXpr::IsVectorAtCompileTime, Math::Vector, + std::conditional_t< + MatrixXpr::ColsAtCompileTime == 1, Math::Vector, + Math::Matrix + > + >; }; - template - struct RangeOf; - template struct RangeOf> { using ResultType = typename ResultOf>::Type; - using Type = - typename Internal::RangeOfSAT< - IsBooleanRange::Value, - IsIntegerRange::Value, - IsRealRange::Value, - IsComplexRange::Value, - IsVectorRange::Value, - IsMatrixRange::Value> - ::Type; + using Type = typename RangeOf::Type; }; template struct RangeOf> { using ResultType = typename ResultOf>::Type; - using Type = - typename Internal::RangeOfSAT< - IsBooleanRange::Value, - IsIntegerRange::Value, - IsRealRange::Value, - IsComplexRange::Value, - IsVectorRange::Value, - IsMatrixRange::Value> - ::Type; + using Type = typename RangeOf::Type; }; } diff --git a/src/Rodin/Variational/Transpose.h b/src/Rodin/Variational/Transpose.h index 8f2b088be..749f9a521 100644 --- a/src/Rodin/Variational/Transpose.h +++ b/src/Rodin/Variational/Transpose.h @@ -33,10 +33,6 @@ namespace Rodin::Variational using OperandRangeType = typename FormLanguage::Traits::RangeType; - static_assert( - std::is_same_v> - || std::is_same_v>); - /** * @brief Constructs the Transpose matrix of the given matrix. */ @@ -57,14 +53,12 @@ namespace Rodin::Variational m_operand(std::move(other.m_operand)) {} - inline constexpr RangeShape getRangeShape() const { return m_operand->getRangeShape().transpose(); } - inline constexpr const OperandType& getOperand() const { @@ -72,22 +66,21 @@ namespace Rodin::Variational return *m_operand; } - inline constexpr auto getValue(const Geometry::Point& p) const { return this->object(getOperand().getValue(p)).transpose(); } - inline + template constexpr - void getValue(Math::Matrix& out, const Geometry::Point& p) const + void getValue(T& out, const Geometry::Point& p) const { getOperand().getValue(out, p); out.transposeInPlace(); } - inline Transpose* copy() const noexcept override + Transpose* copy() const noexcept override { return new Transpose(*this); } @@ -132,28 +125,24 @@ namespace Rodin::Variational m_operand(std::move(other.m_operand)) {} - inline constexpr const auto& getLeaf() const { return getOperand().getLeaf(); } - inline constexpr RangeShape getRangeShape() const { return getOperand().getRangeShape().transpose(); } - inline constexpr size_t getDOFs(const Geometry::Polytope& simplex) const { return getOperand().getDOFs(simplex); } - inline constexpr const OperandType& getOperand() const { @@ -161,34 +150,30 @@ namespace Rodin::Variational return *m_operand; } - inline const Geometry::Point& getPoint() const { return m_operand->getPoint(); } - inline Transpose& setPoint(const Geometry::Point& p) { m_operand->setPoint(p); return *this; } - inline constexpr auto getBasis(size_t local) const { return this->object(getOperand().getBasis(local)).transpose(); } - inline constexpr const FES& getFiniteElementSpace() const { return m_operand.getFiniteElementSpace(); } - inline Transpose* copy() const noexcept override + Transpose* copy() const noexcept override { return new Transpose(*this); } diff --git a/src/Rodin/Variational/UnaryMinus.h b/src/Rodin/Variational/UnaryMinus.h index 23cf00462..e59132c47 100644 --- a/src/Rodin/Variational/UnaryMinus.h +++ b/src/Rodin/Variational/UnaryMinus.h @@ -13,7 +13,6 @@ #include "ForwardDecls.h" #include "Function.h" #include "ShapeFunction.h" -#include "RealFunction.h" #include "LinearFormIntegrator.h" #include "BilinearFormIntegrator.h" @@ -66,14 +65,12 @@ namespace Rodin::Variational m_op(std::move(other.m_op)) {} - inline constexpr RangeShape getRangeShape() const { return getOperand().getRangeShape(); } - inline constexpr const OperandType& getOperand() const { @@ -81,31 +78,19 @@ namespace Rodin::Variational return *m_op; } - inline auto getValue(const Geometry::Point& p) const { return -this->object(getOperand().getValue(p)); } - inline + template constexpr - void getValue(Math::Vector& res, const Geometry::Point& p) const + void getValue(T& res, const Geometry::Point& p) const { - static_assert(FormLanguage::IsVectorRange::Value); getOperand().getValue(res, p); res *= -1; } - inline - constexpr - void getValue(Math::Matrix& res, const Geometry::Point& p) const - { - static_assert(FormLanguage::IsMatrixRange::Value); - getOperand().getValue(res, p); - res *= -1; - } - - inline constexpr UnaryMinus& traceOf(Geometry::Attribute attr) { @@ -114,7 +99,6 @@ namespace Rodin::Variational return *this; } - inline constexpr UnaryMinus& traceOf(const FlatSet& attrs) { @@ -123,7 +107,7 @@ namespace Rodin::Variational return *this; } - inline UnaryMinus* copy() const noexcept override + UnaryMinus* copy() const noexcept override { return new UnaryMinus(*this); } diff --git a/src/Rodin/Variational/Zero.h b/src/Rodin/Variational/Zero.h index 9b58af965..1750e4d17 100644 --- a/src/Rodin/Variational/Zero.h +++ b/src/Rodin/Variational/Zero.h @@ -56,7 +56,7 @@ namespace Rodin::Variational constexpr auto getValue(const Geometry::Point&) const { - return Real(0); + return ScalarType(0); } inline Zero* copy() const noexcept override