From a1869717c37c13ead46d00cd7e4608848b96aeed Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 10:57:22 +0100 Subject: [PATCH 01/34] discrete: fix typo in {DiscreteField} --- lib/mito/discrete/DiscreteField.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mito/discrete/DiscreteField.h b/lib/mito/discrete/DiscreteField.h index 5ebebb933..152b0a5a1 100644 --- a/lib/mito/discrete/DiscreteField.h +++ b/lib/mito/discrete/DiscreteField.h @@ -99,7 +99,7 @@ namespace mito::discrete { inline auto cend() const { return std::cend(_map_entry_to_values); } private: - // the underlying mapping of entries to nodal values + // the underlying mapping of entries to values map_type _map_entry_to_values; // the name of the field From 0bc4785bde4db2ce7c0f412dae90e66ce0209386 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 15:29:01 +0100 Subject: [PATCH 02/34] fields: major redesign of fields Fields are now simply a concept of functions defined on a set of coordinates. This leads to major simplifications in the api. --- benchmarks/mito.lib/fields/laplacian.cc | 2 +- .../mito.lib/integration/integration.cc | 2 +- benchmarks/mito.lib/pdes/poisson.cc | 11 +- extensions/mito/mito.cc | 8 +- lib/mito/fem/DiscreteSystem.h | 22 ++-- .../elements/tri1/IsoparametricTriangleP1.h | 7 +- .../elements/tri2/IsoparametricTriangleP2.h | 6 +- lib/mito/fields/Field.h | 50 -------- lib/mito/fields/api.h | 16 --- lib/mito/fields/differential.h | 10 +- lib/mito/fields/factories.h | 23 +--- lib/mito/fields/fields_algebra.h | 51 +++----- lib/mito/fields/fields_arithmetic.h | 110 ------------------ lib/mito/fields/fields_functions_arithmetic.h | 69 ----------- lib/mito/fields/forward.h | 19 +-- lib/mito/fields/public.h | 6 - lib/mito/geometry/ReferenceSimplex.h | 3 +- lib/mito/geometry/metric_space.h | 12 +- lib/mito/manifolds/factories.h | 4 +- tests/mito.lib/constraints/dirichlet.cc | 2 +- tests/mito.lib/discrete/mesh_field.cc | 2 +- tests/mito.lib/fields/calculus_identities.cc | 2 +- .../mito.lib/fields/calculus_scalar_field.cc | 2 +- tests/mito.lib/fields/fields.cc | 4 +- tests/mito.lib/fields/fields_traits.cc | 7 +- tests/mito.lib/fields/polar_metric_field.cc | 8 +- .../mito.lib/fields/spherical_metric_field.cc | 10 +- .../integration/divergence_theorem.cc | 2 +- .../integration/quadrature_flip_segment_3D.cc | 3 +- .../integration/quadrature_load_mesh_2D.cc | 2 +- .../quadrature_load_mesh_2D_mpi.cc | 2 +- .../io/parallel_vtk_cloud_field_writer.cc | 2 +- .../io/parallel_vtk_mesh_field_writer.cc | 2 +- .../mito.lib/manifolds/euclidean_gradient.cc | 2 +- tests/mito.lib/manifolds/polar_gradient.cc | 2 +- .../mito.lib/manifolds/spherical_gradient.cc | 3 +- .../surface_half_sphere_cartesian.cc | 2 +- .../quadrature/quadrature_segment_1D.cc | 2 +- .../quadrature/quadrature_segment_3D.cc | 2 +- .../quadrature/quadrature_square_3D.cc | 3 +- .../quadrature/quadrature_triangle_2D.cc | 10 +- 41 files changed, 94 insertions(+), 413 deletions(-) delete mode 100644 lib/mito/fields/Field.h delete mode 100644 lib/mito/fields/fields_arithmetic.h delete mode 100644 lib/mito/fields/fields_functions_arithmetic.h diff --git a/benchmarks/mito.lib/fields/laplacian.cc b/benchmarks/mito.lib/fields/laplacian.cc index b61ed709d..86cf9bb13 100644 --- a/benchmarks/mito.lib/fields/laplacian.cc +++ b/benchmarks/mito.lib/fields/laplacian.cc @@ -47,7 +47,7 @@ laplacian_mito(const coordinates_t & x) constexpr auto x1 = mito::functions::component; // a scalar field - constexpr auto f = mito::fields::field(mito::functions::pow<4>(x0 * x1)); + constexpr auto f = mito::functions::pow<4>(x0 * x1); // the gradient of {f} constexpr auto gradient = mito::fields::gradient(f); diff --git a/benchmarks/mito.lib/integration/integration.cc b/benchmarks/mito.lib/integration/integration.cc index f21c564a4..363a07440 100644 --- a/benchmarks/mito.lib/integration/integration.cc +++ b/benchmarks/mito.lib/integration/integration.cc @@ -70,7 +70,7 @@ main() auto manifold = mito::manifolds::manifold(tetra_mesh, coord_system); // instantiate a scalar field - auto f = mito::fields::field(mito::functions::cos(x_0 * x_1)); + auto f = mito::functions::cos(x_0 * x_1); // instantiate a GAUSS integrator with degree of exactness equal to 2 auto integrator = mito::quadrature::integrator(manifold); diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc index f7e30148e..a1b3e6b00 100644 --- a/benchmarks/mito.lib/pdes/poisson.cc +++ b/benchmarks/mito.lib/pdes/poisson.cc @@ -63,7 +63,7 @@ main() auto boundary_mesh = mito::mesh::boundary(mesh); // the zero field - auto zero = mito::fields::field(mito::functions::zero); + auto zero = mito::functions::zero; // set homogeneous Dirichlet boundary condition auto constraints = mito::constraints::dirichlet_bc(boundary_mesh, zero); @@ -75,9 +75,8 @@ main() auto fem_lhs_block = mito::fem::blocks::grad_grad_block(); // the right hand side - auto f = mito::fields::field( - 2.0 * std::numbers::pi * std::numbers::pi * mito::functions::sin(std::numbers::pi * x) - * mito::functions::sin(std::numbers::pi * y)); + auto f = 2.0 * std::numbers::pi * std::numbers::pi * mito::functions::sin(std::numbers::pi * x) + * mito::functions::sin(std::numbers::pi * y); // channel << "Right hand side: " << f(coordinates_t{ 0.5, 0.5 }) << journal::endl; // a source term block @@ -106,8 +105,8 @@ main() solver.destroy(); // the exact solution field - auto u_ex = mito::fields::field( - mito::functions::sin(std::numbers::pi * x) * mito::functions::sin(std::numbers::pi * y)); + auto u_ex = + mito::functions::sin(std::numbers::pi * x) * mito::functions::sin(std::numbers::pi * y); // compute the L2 error auto error_L2 = discrete_system.compute_l2_error(u_ex); diff --git a/extensions/mito/mito.cc b/extensions/mito/mito.cc index bad5c83e1..74733fbb7 100644 --- a/extensions/mito/mito.cc +++ b/extensions/mito/mito.cc @@ -73,9 +73,7 @@ PYBIND11_MODULE(mito, m) // the mito scalar field 2D using scalar_function_2D_t = std::function; - using scalar_function_functor_2D_t = mito::functions::FunctionFromFunctor; - using scalar_field_2D_t = - decltype(mito::fields::field(std::declval())); + using scalar_field_2D_t = mito::functions::FunctionFromFunctor; mito::py::class_(m, "ScalarField2D") // the constructor .def(mito::py::init()) @@ -90,9 +88,7 @@ PYBIND11_MODULE(mito, m) // the mito scalar field 3D using scalar_function_3D_t = std::function; - using scalar_function_functor_3D_t = mito::functions::FunctionFromFunctor; - using scalar_field_3D_t = - decltype(mito::fields::field(std::declval())); + using scalar_field_3D_t = mito::functions::FunctionFromFunctor; mito::py::class_(m, "ScalarField3D") // the constructor .def(mito::py::init()) diff --git a/lib/mito/fem/DiscreteSystem.h b/lib/mito/fem/DiscreteSystem.h index de88210da..1c688788a 100644 --- a/lib/mito/fem/DiscreteSystem.h +++ b/lib/mito/fem/DiscreteSystem.h @@ -246,12 +246,11 @@ namespace mito::fem { auto u_numerical = _assemble_element_solution( element, _solution_field, tensor::make_integer_sequence()); // assemble the elementary error field as a function of the barycentric coordinates - auto u_error = - u_numerical - u_exact.function()(element.parametrization().function()); + auto u_error = u_numerical - u_exact(element.parametrization()); // compute the elementary contribution to the L2 norm - norm += blocks::l2_norm_block(u_error.function()) - .compute(element); + norm += + blocks::l2_norm_block(u_error).compute(element); } return std::sqrt(norm); @@ -297,22 +296,19 @@ namespace mito::fem { auto u_numerical = _assemble_element_solution( element, _solution_field, tensor::make_integer_sequence()); // assemble the elementary error field as a function of the barycentric coordinates - auto u_error = - u_numerical - u_exact.function()(element.parametrization().function()); + auto u_error = u_numerical - u_exact(element.parametrization()); // auto u_numerical_gradient = _assemble_element_solution_gradient( element, _solution_field, tensor::make_integer_sequence()); // assemble the elementary error gradient as a function of the barycentric // coordinates auto u_error_gradient = - u_numerical_gradient - - fields::gradient(u_exact).function()(element.parametrization().function()); + u_numerical_gradient - fields::gradient(u_exact)(element.parametrization()); // compute the elementary contributions to the H1 norm - norm += blocks::l2_norm_block(u_error.function()) - .compute(element) - + blocks::l2_norm_block( - u_error_gradient.function()) - .compute(element); + norm += + blocks::l2_norm_block(u_error).compute(element) + + blocks::l2_norm_block(u_error_gradient) + .compute(element); } return std::sqrt(norm); diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index fcd527d92..b57fcab3f 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -73,8 +73,7 @@ namespace mito::fem { constexpr auto phi_2 = shape_functions.shape<2>(); // return the isoparametric mapping from parametric to physical coordinates - return mito::fields::linear_combination( - std::array{ _x0, _x1, _x2 }, phi_0, phi_1, phi_2); + return functions::linear_combination(std::array{ _x0, _x1, _x2 }, phi_0, phi_1, phi_2); } // get the shape function associated with local node {a} @@ -113,8 +112,8 @@ namespace mito::fem { constexpr auto gradient() const { // assemble the gradient as a function of barycentric coordinates - auto gradient_function = - fields::field([&](const parametric_coordinates_type & xi) -> tensor::vector_t<2> { + auto gradient_function = functions::function( + [&](const parametric_coordinates_type & xi) -> tensor::vector_t<2> { // the jacobian of the mapping from the reference element to the physical // element evaluated at {xi} auto J = jacobian()(xi); diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 8248875df..2c8cf070e 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -74,7 +74,7 @@ namespace mito::fem { constexpr auto phi_5 = shape_functions.shape<5>(); // return the isoparametric mapping from barycentric to physical coordinates - return mito::fields::linear_combination( + return functions::linear_combination( std::array{ _x0, _x1, _x2, x3, x4, x5 }, phi_0, phi_1, phi_2, phi_3, phi_4, phi_5); } @@ -128,8 +128,8 @@ namespace mito::fem { constexpr auto gradient() const { // assemble the gradient as a function of barycentric coordinates - auto gradient_function = - fields::field([&](const parametric_coordinates_type & xi) -> tensor::vector_t<2> { + auto gradient_function = functions::function( + [&](const parametric_coordinates_type & xi) -> tensor::vector_t<2> { // the jacobian of the mapping from the reference element to the physical // element evaluated at {xi} auto J = jacobian()(xi); diff --git a/lib/mito/fields/Field.h b/lib/mito/fields/Field.h deleted file mode 100644 index 9425adb93..000000000 --- a/lib/mito/fields/Field.h +++ /dev/null @@ -1,50 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - - -// code guard -#pragma once - - -namespace mito::fields { - - /* - * A field of vectors or forms. The fields implements {operator()} for an input type - * {coordinates_t}, which returns the value of the field at that physical location. - */ - - template - requires(geometry::coordinates_c) - class Field { - private: - // the type of underlying function - using function_type = F; - - public: - // the type in input to the field (e.g. where the field is defined) - using input_type = function_type::input_type; - // the type of field (e.g. vector field, form field, ...) - using output_type = function_type::output_type; - // the type of coordinates - using coordinates_type = input_type; - - public: - // constructor - constexpr Field(function_type f) : _f{ std::move(f) } {} - - // the value of the field at position {x} - constexpr auto operator()(const coordinates_type & x) const -> output_type { return _f(x); } - - // accessor to the underlying function - constexpr auto function() const -> const function_type & { return _f; } - - private: - // the function assigning the value of the field at each point - function_type _f; - }; -} - - -// end of file diff --git a/lib/mito/fields/api.h b/lib/mito/fields/api.h index 1b7d7becb..4f9e06b91 100644 --- a/lib/mito/fields/api.h +++ b/lib/mito/fields/api.h @@ -9,22 +9,6 @@ namespace mito::fields { - // field alias - template - using field_t = Field; - - // factory for fields - template - constexpr auto field(const F & f) -> field_t; - - // factory for fields - template - constexpr auto field(F && f) -> field_t; - - // factory for fields (from functors) - template - constexpr auto field(F && f) -> field_t; - // uniform field template constexpr auto uniform_field(const Y & constant); diff --git a/lib/mito/fields/differential.h b/lib/mito/fields/differential.h index e18646475..1afbbef12 100644 --- a/lib/mito/fields/differential.h +++ b/lib/mito/fields/differential.h @@ -15,7 +15,7 @@ namespace mito::fields { constexpr auto derivative(const field_c auto & f) { // the {I...}-th first partial derivative - return field(mito::functions::derivative(f.function())); + return mito::functions::derivative(f); } // function to compute the gradient of a scalar field @@ -23,7 +23,7 @@ namespace mito::fields { constexpr auto gradient(const F & field) { // the type of coordinate - using coordinate_t = F::coordinates_type; + using coordinate_t = F::input_type; // the spatial dimension of the field constexpr int D = coordinate_t::dim; @@ -41,7 +41,7 @@ namespace mito::fields { constexpr auto gradient(const F & field) { // the type of coordinate - using coordinate_t = F::coordinates_type; + using coordinate_t = F::input_type; // the spatial dimension of the field constexpr int D = coordinate_t::dim; // the number of components of the vector field @@ -72,7 +72,7 @@ namespace mito::fields { constexpr auto divergence(const F & field) { // the type of coordinate - using coordinate_t = F::coordinates_type; + using coordinate_t = F::input_type; // the spatial dimension of the field constexpr int D = coordinate_t::dim; @@ -91,7 +91,7 @@ namespace mito::fields { constexpr auto divergence(const F & field) { // the type of coordinate - using coordinate_t = F::coordinates_type; + using coordinate_t = F::input_type; // the spatial dimension of the field constexpr int D = coordinate_t::dim; diff --git a/lib/mito/fields/factories.h b/lib/mito/fields/factories.h index 4360a2dc9..62cce8fed 100644 --- a/lib/mito/fields/factories.h +++ b/lib/mito/fields/factories.h @@ -9,32 +9,11 @@ namespace mito::fields { - // factory for fields - template - constexpr auto field(const F & f) -> field_t - { - return field_t(f); - } - - // factory for fields - template - constexpr auto field(F && f) -> field_t - { - return field_t(std::forward(f)); - } - - // factory for fields (from functors) - template - constexpr auto field(F && f) - { - return field(mito::functions::function(std::forward(f))); - } - // uniform field template constexpr auto uniform_field(const Y & constant) { - return field(mito::functions::constant(constant)); + return mito::functions::constant(constant); } } diff --git a/lib/mito/fields/fields_algebra.h b/lib/mito/fields/fields_algebra.h index b9c1dece7..b467fa3af 100644 --- a/lib/mito/fields/fields_algebra.h +++ b/lib/mito/fields/fields_algebra.h @@ -10,39 +10,14 @@ // Algebraic operations on Fields namespace mito::fields { - // f^(-1) - template - constexpr auto inverse(const F & f) - { - using coordinates_type = typename F::coordinates_type; - return field( - functions::function([f](const coordinates_type & x) { return inverse(f(x)); })); - } - // det(f) template constexpr auto determinant(const F & f) { - using coordinates_type = typename F::coordinates_type; - return field(functions::function([f](const coordinates_type & x) -> mito::tensor::scalar_t { + using coordinates_type = typename F::input_type; + return functions::function([f](const coordinates_type & x) -> mito::tensor::scalar_t { return determinant(f(x)); - })); - } - - // f transpose - template - constexpr auto transpose(const F & f) - { - return field(functions::transpose(f.function())); - } - - // sqrt(f) - template - constexpr auto sqrt(const F & f) - { - using coordinates_type = typename F::coordinates_type; - return field(functions::function( - [f](const coordinates_type & x) -> mito::tensor::scalar_t { return std::sqrt(f(x)); })); + }); } // the tensor product of two fields of one-forms @@ -50,9 +25,9 @@ namespace mito::fields { constexpr auto tensor(const F1 & fA, const F2 & fB) requires(compatible_fields_c) { - using coordinates_type = typename F1::coordinates_type; - return field(functions::function( - [fA, fB](const coordinates_type & x) -> auto { return tensor::tensor(fA(x), fB(x)); })); + using coordinates_type = typename F1::input_type; + return functions::function( + [fA, fB](const coordinates_type & x) -> auto { return tensor::tensor(fA(x), fB(x)); }); } // the tensor product of three fields of one-forms @@ -60,7 +35,7 @@ namespace mito::fields { constexpr auto tensor(const F1 & fA, const F2 & fB, const F3 & fC) requires(compatible_fields_c && compatible_fields_c) { - using coordinates_type = typename F1::coordinates_type; + using coordinates_type = typename F1::input_type; return field(functions::function([fA, fB, fC](const coordinates_type & x) -> auto { return tensor::tensor(fA(x), fB(x), fC(x)); })); @@ -78,9 +53,9 @@ namespace mito::fields { constexpr auto wedge(const F1 & fA, const F2 & fB) requires(compatible_fields_c) { - using coordinates_type = typename F1::coordinates_type; - return field(functions::function( - [fA, fB](const coordinates_type & x) { return tensor::wedge(fA(x), fB(x)); })); + using coordinates_type = typename F1::input_type; + return functions::function( + [fA, fB](const coordinates_type & x) { return tensor::wedge(fA(x), fB(x)); }); } // the wedge product of three fields of one-forms @@ -88,10 +63,10 @@ namespace mito::fields { constexpr auto wedge(const F1 & fA, const F2 & fB, const F3 & fC) requires(compatible_fields_c && compatible_fields_c) { - using coordinates_type = typename F1::coordinates_type; - return field(functions::function([fA, fB, fC](const coordinates_type & x) { + using coordinates_type = typename F1::input_type; + return functions::function([fA, fB, fC](const coordinates_type & x) { return tensor::wedge(fA(x), fB(x), fC(x)); - })); + }); } } diff --git a/lib/mito/fields/fields_arithmetic.h b/lib/mito/fields/fields_arithmetic.h deleted file mode 100644 index 01a73d102..000000000 --- a/lib/mito/fields/fields_arithmetic.h +++ /dev/null @@ -1,110 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -// Arithmetic operations on Fields -namespace mito::fields { - - // addition of N fields fa + ... + fn - template - constexpr auto linear_combination(std::array coeffs, const Fs &... fs) - requires(compatible_fields_c) - { - return field(functions::linear_combination(coeffs, fs.function()...)); - } - - // addition of fields fa + fb - template - constexpr auto operator+(const F1 & fA, const F2 & fB) - requires(compatible_fields_c) - { - return field(fA.function() + fB.function()); - } - - // scalar + field - template - constexpr auto operator+(const tensor::tensor_or_scalar_c auto & a, const F & f) - { - return field(a + f.function()); - } - - // field + scalar - constexpr auto operator+(const field_c auto & f, const tensor::tensor_or_scalar_c auto & a) - { - return a + f; - } - - // scalar * fields - template - constexpr auto operator*(const tensor::tensor_or_scalar_c auto & a, const F & f) - { - return field(a * f.function()); - } - - // field * scalar - constexpr auto operator*(const field_c auto & f, const tensor::tensor_or_scalar_c auto & a) - { - return a * f; - } - - // product of fields - template - constexpr auto operator*(const F1 & fA, const F2 & fB) - requires(compatible_fields_c) - { - return field(fA.function() * fB.function()); - } - - // field / a - constexpr auto operator/(const field_c auto & f, const tensor::scalar_t & a) - { - return (1.0 / a) * f; - } - - // unary operator- for fields - constexpr auto operator-(const field_c auto & f) - { - return -1.0 * f; - } - - // subtraction of fields fa - fb - constexpr auto operator-(const field_c auto & fA, const field_c auto & fB) - { - return fA + (-fB); - } - - // a - field - constexpr auto operator-(const tensor::tensor_or_scalar_c auto & a, const field_c auto & f) - { - return a + (-f); - } - - // field - a - constexpr auto operator-(const field_c auto & f, const tensor::tensor_or_scalar_c auto & a) - { - return f + (-a); - } - - // a / field - template - constexpr auto operator/(const tensor::tensor_or_scalar_c auto & a, const F & f) - { - return field(a / f.function()); - } - - // field1 / field2 - template - constexpr auto operator/(const F1 & f1, const F2 & f2) - requires(compatible_fields_c) - { - return field(f1.function() / f2.function()); - } -} - - -// end of file diff --git a/lib/mito/fields/fields_functions_arithmetic.h b/lib/mito/fields/fields_functions_arithmetic.h deleted file mode 100644 index bb479e484..000000000 --- a/lib/mito/fields/fields_functions_arithmetic.h +++ /dev/null @@ -1,69 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -// Arithmetic hybrid operations on Fields and Functions -namespace mito::fields { - - // addition of fields fa + fb - template - constexpr auto operator+(const F1 & fA, const F2 & fB) - { - return field(fA.function() + fB); - } - - // addition of fields fa + fb - template - constexpr auto operator+(const F1 & fA, const F2 & fB) - { - return fB + fA; - } - - // product of fields - template - constexpr auto operator*(const F1 & fA, const F2 & fB) - { - return field(fA.function() * fB); - } - - // product of fields - template - constexpr auto operator*(const F1 & fA, const F2 & fB) - { - return field(fA * fB.function()); - } - - // subtraction of fields fa - fb - constexpr auto operator-(const field_c auto & fA, const functions::function_c auto & fB) - { - return fA + (-fB); - } - - // subtraction of fields fa - fb - constexpr auto operator-(const functions::function_c auto & fA, const field_c auto & fB) - { - return fA + (-fB); - } - - // field1 / field2 - template - constexpr auto operator/(const F1 & f1, const F2 & f2) - { - return field(f1.function() / f2); - } - - // field1 / field2 - template - constexpr auto operator/(const F1 & f1, const F2 & f2) - { - return field(f1 / f2.function()); - } -} - - -// end of file diff --git a/lib/mito/fields/forward.h b/lib/mito/fields/forward.h index 09749055d..f133dee19 100644 --- a/lib/mito/fields/forward.h +++ b/lib/mito/fields/forward.h @@ -9,18 +9,9 @@ namespace mito::fields { - // class field - template - requires(geometry::coordinates_c) - class Field; - - // concept of a field - template - concept field_c = requires(FIELD c) { - // require that F only binds to {Field} specializations or derived classes - [](const Field &) { - }(c); - }; + // concept of a field: a field is a function of coordinates + template + concept field_c = functions::function_c && geometry::coordinates_c; // concept of a field {FIELD} being a scalar field template @@ -64,8 +55,8 @@ namespace mito::fields { concept compatible_fields_c = (sizeof...(FIELDS) <= 1) || // trivially true for 0 or 1 (std::same_as< - typename std::tuple_element_t<0, std::tuple>::coordinates_type, - typename FIELDS::coordinates_type> + typename std::tuple_element_t<0, std::tuple>::input_type, + typename FIELDS::input_type> && ...); } diff --git a/lib/mito/fields/public.h b/lib/mito/fields/public.h index fc412ce97..9aa3e76e2 100644 --- a/lib/mito/fields/public.h +++ b/lib/mito/fields/public.h @@ -16,15 +16,9 @@ // published type factories; this is the file you are looking for... #include "api.h" -// classes implementation -#include "Field.h" - // factories implementation #include "factories.h" -// arithmetic operations on fields -#include "fields_arithmetic.h" -#include "fields_functions_arithmetic.h" // algebraic operations on fields #include "fields_algebra.h" // differential calculus on fields diff --git a/lib/mito/geometry/ReferenceSimplex.h b/lib/mito/geometry/ReferenceSimplex.h index 6b94a0bd4..4c9c7b403 100644 --- a/lib/mito/geometry/ReferenceSimplex.h +++ b/lib/mito/geometry/ReferenceSimplex.h @@ -24,8 +24,7 @@ namespace mito::geometry { // the function extracting the I component of a parametric point template - static constexpr auto xi = - fields::field(functions::component); + static constexpr auto xi = functions::component; }; } // namespace mito diff --git a/lib/mito/geometry/metric_space.h b/lib/mito/geometry/metric_space.h index 75ef7c4f7..250cb8293 100644 --- a/lib/mito/geometry/metric_space.h +++ b/lib/mito/geometry/metric_space.h @@ -36,7 +36,7 @@ namespace mito::geometry { static constexpr auto g = metric::field(); // the inverse metric field in coordinates {coordinates_type} - static constexpr auto g_inv = fields::inverse(g); + static constexpr auto g_inv = functions::inverse(g); // get the I-th basis element for one-form fields template @@ -56,7 +56,7 @@ namespace mito::geometry { public: // the metric volume form static constexpr auto w = - fields::sqrt(fields::determinant(g)) * _wedge(tensor::make_integer_sequence{}); + functions::sqrt(fields::determinant(g)) * _wedge(tensor::make_integer_sequence{}); // get the metric equivalent vector field to a given one-form field static constexpr auto metric_equivalent(const fields::one_form_field_c auto & one_form); @@ -71,7 +71,7 @@ namespace mito::geometry { const fields::one_form_field_c auto & one_form) { // return a vector field that, once evaluated at {x}... - return fields::field(functions::function([one_form](const coordinates_type & x) { + return functions::function([one_form](const coordinates_type & x) { // returns the contraction of the inverse metric with the components of the one // form auto _one_form_components = [one_form, x](tensor::integer_sequence) { @@ -80,7 +80,7 @@ namespace mito::geometry { // all done return _one_form_components(tensor::make_integer_sequence{}); - })); + }); } template @@ -98,10 +98,10 @@ namespace mito::geometry { requires(fields::compatible_fields_c) { // return a one form field that, once evaluated at {x}... - return fields::field(functions::function([vector, matrix](const coordinates_type & x) { + return functions::function([vector, matrix](const coordinates_type & x) { // returns the metric equivalent form to {vector(x)} return tensor::one_form(vector(x), matrix(x)); - })); + }); } } // namespace mito diff --git a/lib/mito/manifolds/factories.h b/lib/mito/manifolds/factories.h index 063bcc537..eed61c477 100644 --- a/lib/mito/manifolds/factories.h +++ b/lib/mito/manifolds/factories.h @@ -72,7 +72,7 @@ namespace mito::manifolds { // if the submanifold is a line (either in 2D or 3D) then use one placeholder {_} if constexpr (mesh_type::order == 1) { // the restriction of the metric volume form to the submanifold - auto wS = fields::field( + auto wS = functions::function( [w, normal_field...](const coordsT & x) { return w(x)(normal_field(x)..., _); }); // return the manifold @@ -82,7 +82,7 @@ namespace mito::manifolds { // if the submanifold is a surface (in 3D) then use two placeholders {_, _} if constexpr (mesh_type::order == 2) { // the restriction of the metric volume element to the submanifold - auto wS = fields::field( + auto wS = functions::function( [w, normal_field...](const coordsT & x) { return w(x)(normal_field(x)..., _, _); }); // return the manifold diff --git a/tests/mito.lib/constraints/dirichlet.cc b/tests/mito.lib/constraints/dirichlet.cc index 871ef62a6..cbd778a1f 100644 --- a/tests/mito.lib/constraints/dirichlet.cc +++ b/tests/mito.lib/constraints/dirichlet.cc @@ -34,7 +34,7 @@ TEST(Constraints, Dirichlet) auto mesh = mito::io::summit::reader(fileStream, coord_system); // the value of the Dirichlet boundary condition - auto bc_value = mito::fields::field(mito::functions::sqrt(x * x + y * y)); + auto bc_value = mito::functions::sqrt(x * x + y * y); // get all the nodes on the mesh boundary auto boundary_mesh = mito::mesh::boundary(mesh); diff --git a/tests/mito.lib/discrete/mesh_field.cc b/tests/mito.lib/discrete/mesh_field.cc index 5eee2910d..ad31e3f19 100644 --- a/tests/mito.lib/discrete/mesh_field.cc +++ b/tests/mito.lib/discrete/mesh_field.cc @@ -26,7 +26,7 @@ TEST(Discretization, NodalFieldSphere) auto mesh_field = mito::discrete::mesh_field>(mesh, "normal"); // the normal field to the submanifold - constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { + constexpr auto normal_field = mito::functions::function([](const coordinates_t & x) -> auto { mito::tensor::scalar_t phi = std::atan2(x[1], x[0]); mito::tensor::scalar_t theta = std::atan2(std::hypot(x[1], x[0]), x[2]); return std::sin(theta) * std::cos(phi) * mito::tensor::e_0<3> diff --git a/tests/mito.lib/fields/calculus_identities.cc b/tests/mito.lib/fields/calculus_identities.cc index ba7458e20..e395dd7ec 100644 --- a/tests/mito.lib/fields/calculus_identities.cc +++ b/tests/mito.lib/fields/calculus_identities.cc @@ -49,7 +49,7 @@ TEST(Identities, DivGrad) // the divergence of the gradient transposed of {f} constexpr auto div_grad_T = - mito::fields::divergence(mito::fields::transpose(mito::fields::gradient(f))); + mito::fields::divergence(mito::functions::transpose(mito::fields::gradient(f))); // the gradient of the divergence of {f} constexpr auto grad_div = mito::fields::gradient(mito::fields::divergence(f)); diff --git a/tests/mito.lib/fields/calculus_scalar_field.cc b/tests/mito.lib/fields/calculus_scalar_field.cc index bd49d1852..c46377a20 100644 --- a/tests/mito.lib/fields/calculus_scalar_field.cc +++ b/tests/mito.lib/fields/calculus_scalar_field.cc @@ -31,7 +31,7 @@ TEST(Laplacian, ScalarFields) constexpr auto x1 = mito::functions::component; // a scalar field - constexpr auto f = mito::fields::field(sin(x0 * x1)); + constexpr auto f = sin(x0 * x1); // the gradient of {f} constexpr auto gradient = mito::fields::gradient(f); diff --git a/tests/mito.lib/fields/fields.cc b/tests/mito.lib/fields/fields.cc index a5a8793dd..8a17eafed 100644 --- a/tests/mito.lib/fields/fields.cc +++ b/tests/mito.lib/fields/fields.cc @@ -21,7 +21,7 @@ TEST(Fields, VectorFields) constexpr auto x1 = mito::functions::component; // a scalar field - constexpr auto f = mito::fields::field(cos(x0 * x1)); + constexpr auto f = cos(x0 * x1); // a point in space constexpr auto x = mito::geometry::coordinates({ 0.0, 0.0 }); @@ -36,7 +36,7 @@ TEST(Fields, VectorFields) constexpr auto e1 = mito::functions::constant(mito::tensor::e_1<2>); // a vector field - constexpr auto g = mito::fields::field(cos(x0 * x1) * (e0 + e1)); + constexpr auto g = cos(x0 * x1) * (e0 + e1); // check value of field at {x} static_assert(g(x) == mito::tensor::vector_t<2>{ 1.0, 1.0 }); diff --git a/tests/mito.lib/fields/fields_traits.cc b/tests/mito.lib/fields/fields_traits.cc index 1fa549f8e..ca682301f 100644 --- a/tests/mito.lib/fields/fields_traits.cc +++ b/tests/mito.lib/fields/fields_traits.cc @@ -17,11 +17,8 @@ static constexpr auto x1 = mito::functions::component; TEST(Fields, Traits) { - // cos(x0 * x1) - constexpr auto cos = mito::functions::cos(x0 * x1); - // a scalar field - constexpr auto f = mito::fields::field(cos); + constexpr auto f = mito::functions::cos(x0 * x1); // assert that {f} is a scalar field static_assert(mito::fields::scalar_field_c); @@ -35,7 +32,7 @@ TEST(Fields, Traits) mito::fields::uniform_field(mito::tensor::e_1<2>); // a vector field - constexpr auto g = mito::fields::field(cos * (e0 + e1)); + constexpr auto g = f * (e0 + e1); // assert that {g} is a vector field static_assert(mito::fields::vector_field_c); diff --git a/tests/mito.lib/fields/polar_metric_field.cc b/tests/mito.lib/fields/polar_metric_field.cc index 1bf087a30..ca65c7c2c 100644 --- a/tests/mito.lib/fields/polar_metric_field.cc +++ b/tests/mito.lib/fields/polar_metric_field.cc @@ -31,13 +31,13 @@ TEST(Fields, PolarCoordinates) constexpr auto g = (e_r * e_r) * e_rr + (e_t * e_t) * e_tt; // the inverse metric field - constexpr auto g_inv = mito::fields::inverse(g); + constexpr auto g_inv = mito::functions::inverse(g); // the basis one-forms - constexpr auto dr = mito::fields::field( + constexpr auto dr = mito::functions::function( [e_r, g_inv](const coordinates_t & x) { return mito::tensor::one_form(e_r(x), g_inv(x)); }); - constexpr auto dt = mito::fields::field( + constexpr auto dt = mito::functions::function( [e_t, g_inv](const coordinates_t & x) { return mito::tensor::one_form(e_t(x), g_inv(x)); }); // a point in space @@ -53,7 +53,7 @@ TEST(Fields, PolarCoordinates) // the metric volume element constexpr auto w = - mito::fields::sqrt(mito::fields::determinant(g)) * mito::fields::wedge(dr, dt); + mito::functions::sqrt(mito::fields::determinant(g)) * mito::fields::wedge(dr, dt); constexpr auto dr_scalar = mito::tensor::scalar_t{ 0.01 }; constexpr auto dt_scalar = mito::tensor::scalar_t{ 0.01 }; diff --git a/tests/mito.lib/fields/spherical_metric_field.cc b/tests/mito.lib/fields/spherical_metric_field.cc index 841175d71..8581aac90 100644 --- a/tests/mito.lib/fields/spherical_metric_field.cc +++ b/tests/mito.lib/fields/spherical_metric_field.cc @@ -35,14 +35,14 @@ TEST(Manifolds, SphericalCoordinates) constexpr auto g = (e_r * e_r) * e_rr + (e_t * e_t) * e_tt + (e_p * e_p) * e_pp; // the inverse metric field - constexpr auto g_inv = mito::fields::inverse(g); + constexpr auto g_inv = mito::functions::inverse(g); // the basis one-forms - constexpr auto dr = mito::fields::field( + constexpr auto dr = mito::functions::function( [e_r, g_inv](const coordinates_t & x) { return mito::tensor::one_form(e_r(x), g_inv(x)); }); - constexpr auto dt = mito::fields::field( + constexpr auto dt = mito::functions::function( [e_t, g_inv](const coordinates_t & x) { return mito::tensor::one_form(e_t(x), g_inv(x)); }); - constexpr auto dp = mito::fields::field( + constexpr auto dp = mito::functions::function( [e_p, g_inv](const coordinates_t & x) { return mito::tensor::one_form(e_p(x), g_inv(x)); }); // a point in space @@ -64,7 +64,7 @@ TEST(Manifolds, SphericalCoordinates) // the metric volume element constexpr auto w = - mito::fields::sqrt(mito::fields::determinant(g)) * mito::fields::wedge(dr, dt, dp); + mito::functions::sqrt(mito::fields::determinant(g)) * mito::fields::wedge(dr, dt, dp); constexpr auto dr_scalar = mito::tensor::scalar_t{ 0.01 }; constexpr auto dt_scalar = mito::tensor::scalar_t{ 0.01 }; diff --git a/tests/mito.lib/integration/divergence_theorem.cc b/tests/mito.lib/integration/divergence_theorem.cc index 36d3cb105..654588d84 100644 --- a/tests/mito.lib/integration/divergence_theorem.cc +++ b/tests/mito.lib/integration/divergence_theorem.cc @@ -87,7 +87,7 @@ TEST(DivergenceTheorem, Mesh2D) // the normals calculations // the normal to the boundary - constexpr auto n = mito::fields::field([](const coordinates_t & x) { + constexpr auto n = mito::functions::function([](const coordinates_t & x) { // the left, right, bottom and top normals return (x[0] == 0.0) * (-e0(x)) + (x[0] == 1.0) * e0(x) + (x[1] == 0.0) * (-e1(x)) + (x[1] == 1.0) * e1(x); diff --git a/tests/mito.lib/integration/quadrature_flip_segment_3D.cc b/tests/mito.lib/integration/quadrature_flip_segment_3D.cc index b68be6d97..af5185ff0 100644 --- a/tests/mito.lib/integration/quadrature_flip_segment_3D.cc +++ b/tests/mito.lib/integration/quadrature_flip_segment_3D.cc @@ -43,7 +43,8 @@ TEST(Quadrature, FlipSegment) constexpr auto normal_field_2 = mito::fields::uniform_field(v_3); // the integrand - auto f = mito::fields::field([](const coordinates_t & x) { return std::cos(x[0] * x[1]); }); + auto f = + mito::functions::function([](const coordinates_t & x) { return std::cos(x[0] * x[1]); }); // a mesh with {segment0} auto mesh = mito::mesh::mesh>(); diff --git a/tests/mito.lib/integration/quadrature_load_mesh_2D.cc b/tests/mito.lib/integration/quadrature_load_mesh_2D.cc index fab877789..009fd9b23 100644 --- a/tests/mito.lib/integration/quadrature_load_mesh_2D.cc +++ b/tests/mito.lib/integration/quadrature_load_mesh_2D.cc @@ -30,7 +30,7 @@ TEST(Quadrature, LoadMeshTriangles) auto manifold = mito::manifolds::manifold(mesh, coord_system); // instantiate a scalar field - auto f = mito::fields::field(mito::functions::cos(x_0 * x_1)); + auto f = mito::functions::cos(x_0 * x_1); // instantiate a GAUSS integrator with degree of exactness equal to 2 auto integrator = mito::quadrature::integrator(manifold); diff --git a/tests/mito.lib/integration/quadrature_load_mesh_2D_mpi.cc b/tests/mito.lib/integration/quadrature_load_mesh_2D_mpi.cc index 5022dc60d..fa997c5b8 100644 --- a/tests/mito.lib/integration/quadrature_load_mesh_2D_mpi.cc +++ b/tests/mito.lib/integration/quadrature_load_mesh_2D_mpi.cc @@ -49,7 +49,7 @@ TEST(Quadrature, LoadMeshTrianglesMPI) auto manifold = mito::manifolds::manifold(mesh_partition, coord_system); // instantiate a scalar field - auto f = mito::fields::field(mito::functions::cos(x_0 * x_1)); + auto f = mito::functions::cos(x_0 * x_1); // instantiate a GAUSS integrator with degree of exactness equal to 2 auto integrator = mito::quadrature::integrator(manifold); diff --git a/tests/mito.lib/io/parallel_vtk_cloud_field_writer.cc b/tests/mito.lib/io/parallel_vtk_cloud_field_writer.cc index 4043095e0..20efb3db3 100644 --- a/tests/mito.lib/io/parallel_vtk_cloud_field_writer.cc +++ b/tests/mito.lib/io/parallel_vtk_cloud_field_writer.cc @@ -62,7 +62,7 @@ TEST(ParallelVtkWriter, CloudField) auto point_field = mito::discrete::point_field>(cloud, "normal"); // the normal field to the submanifold - constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { + constexpr auto normal_field = mito::functions::function([](const coordinates_t & x) -> auto { mito::tensor::scalar_t phi = std::atan2(x[1], x[0]); mito::tensor::scalar_t theta = std::atan2(std::hypot(x[1], x[0]), x[2]); return std::sin(theta) * std::cos(phi) * mito::tensor::e_0<3> diff --git a/tests/mito.lib/io/parallel_vtk_mesh_field_writer.cc b/tests/mito.lib/io/parallel_vtk_mesh_field_writer.cc index bc2ee42af..e468aa4a0 100644 --- a/tests/mito.lib/io/parallel_vtk_mesh_field_writer.cc +++ b/tests/mito.lib/io/parallel_vtk_mesh_field_writer.cc @@ -42,7 +42,7 @@ TEST(ParallelVtkWriter, MeshField) auto mesh_field = mito::discrete::mesh_field>(mesh, "normal"); // the normal field to the ball - constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { + constexpr auto normal_field = mito::functions::function([](const coordinates_t & x) -> auto { mito::tensor::scalar_t phi = std::atan2(x[1], x[0]); mito::tensor::scalar_t theta = std::atan2(std::hypot(x[1], x[0]), x[2]); return std::sin(theta) * std::cos(phi) * mito::tensor::e_0<3> diff --git a/tests/mito.lib/manifolds/euclidean_gradient.cc b/tests/mito.lib/manifolds/euclidean_gradient.cc index 42a54fa00..d29d7ea3b 100644 --- a/tests/mito.lib/manifolds/euclidean_gradient.cc +++ b/tests/mito.lib/manifolds/euclidean_gradient.cc @@ -29,7 +29,7 @@ TEST(Manifolds, CartesianGradient) constexpr auto x_1 = mito::geometry::cartesian::x_1<2>; // a scalar field - constexpr auto f = mito::fields::field(x_0 * x_1); + constexpr auto f = x_0 * x_1; // df/dx[0] = x1 constexpr auto df0 = mito::fields::derivative<0>(f); diff --git a/tests/mito.lib/manifolds/polar_gradient.cc b/tests/mito.lib/manifolds/polar_gradient.cc index 4e009804f..e800ecc93 100644 --- a/tests/mito.lib/manifolds/polar_gradient.cc +++ b/tests/mito.lib/manifolds/polar_gradient.cc @@ -33,7 +33,7 @@ TEST(Manifolds, PolarGradient) constexpr auto theta = mito::geometry::polar::theta; // a scalar field - constexpr auto f = mito::fields::field(r * mito::functions::sin(theta)); + constexpr auto f = r * mito::functions::sin(theta); // df/dr = sin(theta) constexpr auto df0 = mito::fields::derivative<0>(f); diff --git a/tests/mito.lib/manifolds/spherical_gradient.cc b/tests/mito.lib/manifolds/spherical_gradient.cc index 84e4e8995..668854a6e 100644 --- a/tests/mito.lib/manifolds/spherical_gradient.cc +++ b/tests/mito.lib/manifolds/spherical_gradient.cc @@ -38,8 +38,7 @@ TEST(Manifolds, SphericalGradient) constexpr auto phi = mito::geometry::spherical::phi; // a scalar field - constexpr auto f = - mito::fields::field(r * mito::functions::sin(theta) * mito::functions::cos(phi)); + constexpr auto f = r * mito::functions::sin(theta) * mito::functions::cos(phi); // df/dr = sin(theta) * cos(phi) constexpr auto df0 = mito::fields::derivative<0>(f); diff --git a/tests/mito.lib/manifolds/surface_half_sphere_cartesian.cc b/tests/mito.lib/manifolds/surface_half_sphere_cartesian.cc index f40e54c8f..48534d42d 100644 --- a/tests/mito.lib/manifolds/surface_half_sphere_cartesian.cc +++ b/tests/mito.lib/manifolds/surface_half_sphere_cartesian.cc @@ -22,7 +22,7 @@ TEST(Manifolds, HalfSphereCartesian) auto mesh = mito::io::summit::reader>(fileStream, coord_system); // the normal field to the submanifold - constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { + constexpr auto normal_field = mito::functions::function([](const coordinates_t & x) -> auto { mito::tensor::scalar_t phi = std::atan2(x[1], x[0]); mito::tensor::scalar_t theta = std::atan2(std::hypot(x[1], x[0]), x[2]); return std::sin(theta) * std::cos(phi) * mito::tensor::e_0<3> diff --git a/tests/mito.lib/quadrature/quadrature_segment_1D.cc b/tests/mito.lib/quadrature/quadrature_segment_1D.cc index efa557a49..cec2e10ba 100644 --- a/tests/mito.lib/quadrature/quadrature_segment_1D.cc +++ b/tests/mito.lib/quadrature/quadrature_segment_1D.cc @@ -38,7 +38,7 @@ TEST(Quadrature, Segment) auto integrator = mito::quadrature::integrator(manifold); // a scalar function - auto f_exp = mito::fields::field(mito::functions::exp(-x_0)); + auto f_exp = mito::functions::exp(-x_0); // integrate exp(-x) on (0, 1) auto integral = integrator.integrate(f_exp); diff --git a/tests/mito.lib/quadrature/quadrature_segment_3D.cc b/tests/mito.lib/quadrature/quadrature_segment_3D.cc index 670645e67..3559d84ec 100644 --- a/tests/mito.lib/quadrature/quadrature_segment_3D.cc +++ b/tests/mito.lib/quadrature/quadrature_segment_3D.cc @@ -42,7 +42,7 @@ TEST(Quadrature, Segment3D) constexpr auto normal_field_2 = mito::fields::uniform_field(v_3); // the integrand - auto f = mito::fields::field([](const coordinates_t & x) { return x[0] * x[1]; }); + auto f = mito::functions::function([](const coordinates_t & x) { return x[0] * x[1]; }); // a mesh with {segment0} auto mesh = mito::mesh::mesh>(); diff --git a/tests/mito.lib/quadrature/quadrature_square_3D.cc b/tests/mito.lib/quadrature/quadrature_square_3D.cc index 8c4967007..aeff27399 100644 --- a/tests/mito.lib/quadrature/quadrature_square_3D.cc +++ b/tests/mito.lib/quadrature/quadrature_square_3D.cc @@ -77,7 +77,8 @@ TEST(Quadrature, Square) auto integrator = mito::quadrature::integrator(manifold); // a scalar field - auto f_xy = mito::fields::field([](const coordinates_t & x) -> real { return x[0] * x[1]; }); + auto f_xy = + mito::functions::function([](const coordinates_t & x) -> real { return x[0] * x[1]; }); // integrate the field real result = integrator.integrate(f_xy); diff --git a/tests/mito.lib/quadrature/quadrature_triangle_2D.cc b/tests/mito.lib/quadrature/quadrature_triangle_2D.cc index 3322cc9c1..ad3d6c438 100644 --- a/tests/mito.lib/quadrature/quadrature_triangle_2D.cc +++ b/tests/mito.lib/quadrature/quadrature_triangle_2D.cc @@ -68,7 +68,7 @@ TEST(Quadrature, Square) mito::quadrature::integrator(bodyManifold); // a scalar field - auto f = mito::fields::field(mito::functions::cos(x_0 * x_1)); + auto f = mito::functions::cos(x_0 * x_1); // integrate the field real result = bodyIntegrator.integrate(f); // the exact solution @@ -80,7 +80,7 @@ TEST(Quadrature, Square) EXPECT_NEAR(result, exact, 1.e-3); // a scalar field - auto f_one = mito::fields::field(mito::functions::constant(1.0)); + auto f_one = mito::functions::constant(1.0); // integrate the field result = bodyIntegrator.integrate(f_one); // exact 1.0 // report @@ -90,7 +90,7 @@ TEST(Quadrature, Square) EXPECT_DOUBLE_EQ(result, 1.0); // a scalar field - auto f_linear = mito::fields::field(x_0); + auto f_linear = x_0; // integrate the field result = bodyIntegrator.integrate(f_linear); // exact 0.5 // report @@ -100,7 +100,7 @@ TEST(Quadrature, Square) EXPECT_DOUBLE_EQ(result, 0.5); // a scalar function - auto f_xy = mito::fields::field(x_0 * x_1); + auto f_xy = x_0 * x_1; // integrate the field result = bodyIntegrator.integrate(f_xy); // exact 0.25 // report @@ -110,7 +110,7 @@ TEST(Quadrature, Square) EXPECT_DOUBLE_EQ(result, 0.25); // a scalar function - auto f_xx = mito::fields::field(x_0 * x_0); + auto f_xx = x_0 * x_0; // integrate the field result = bodyIntegrator.integrate(f_xx); // exact 1.0 / 3.0 // report From 3db6c2d2e29ec8e923206d175efc09d7965458c6 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 15:34:39 +0100 Subject: [PATCH 03/34] fields: remove {derivative} from {fields} namespace It was redundant with respect to {derivative} from the {functions} namespace. --- lib/mito/fields/differential.h | 20 ++++++++----------- .../mito.lib/manifolds/euclidean_gradient.cc | 4 ++-- tests/mito.lib/manifolds/polar_gradient.cc | 4 ++-- .../mito.lib/manifolds/spherical_gradient.cc | 6 +++--- 4 files changed, 15 insertions(+), 19 deletions(-) diff --git a/lib/mito/fields/differential.h b/lib/mito/fields/differential.h index 1afbbef12..b83c2ec20 100644 --- a/lib/mito/fields/differential.h +++ b/lib/mito/fields/differential.h @@ -10,14 +10,6 @@ // Differential operators on Fields namespace mito::fields { - // the {I...}-th first partial derivative of a field - template - constexpr auto derivative(const field_c auto & f) - { - // the {I...}-th first partial derivative - return mito::functions::derivative(f); - } - // function to compute the gradient of a scalar field template constexpr auto gradient(const F & field) @@ -30,7 +22,9 @@ namespace mito::fields { // helper function to compute the gradient of a scalar constexpr auto _grad = [](const F & field, std::index_sequence) { // the vector of the partial derivatives - return ((derivative(field) * uniform_field(tensor::e)) + ...); + return ( + (functions::derivative(field) * uniform_field(tensor::e)) + + ...); }; return _grad(field, std::make_index_sequence{}); @@ -54,7 +48,7 @@ namespace mito::fields { // the tensor of the partial derivatives // (\partial field_I / \partial x_K ) * e_IK, I = 0, ..., N-1, K = 0, ..., D-1 return ( - (derivative(field * uniform_field(tensor::e)) + (functions::derivative(field * uniform_field(tensor::e)) * uniform_field(tensor::unit, I, K>)) + ...); }; @@ -79,7 +73,9 @@ namespace mito::fields { // helper function to compute the divergence of a vector field constexpr auto _div = [](const F & field, std::index_sequence) { // the summation of the I-th partial derivative of the I-th component - return ((derivative(field * uniform_field(tensor::e))) + ...); + return ( + (functions::derivative(field * uniform_field(tensor::e))) + + ...); }; // all done @@ -102,7 +98,7 @@ namespace mito::fields { // the vector of the partial derivatives // (\partial field_KI / \partial x_I) * e_K return ( - (derivative( + (functions::derivative( (field * uniform_field(tensor::e)) * uniform_field(tensor::e)) * uniform_field(tensor::e)) diff --git a/tests/mito.lib/manifolds/euclidean_gradient.cc b/tests/mito.lib/manifolds/euclidean_gradient.cc index d29d7ea3b..e2cda7bcc 100644 --- a/tests/mito.lib/manifolds/euclidean_gradient.cc +++ b/tests/mito.lib/manifolds/euclidean_gradient.cc @@ -32,10 +32,10 @@ TEST(Manifolds, CartesianGradient) constexpr auto f = x_0 * x_1; // df/dx[0] = x1 - constexpr auto df0 = mito::fields::derivative<0>(f); + constexpr auto df0 = mito::functions::derivative<0>(f); // df/dx[1] = x0 - constexpr auto df1 = mito::fields::derivative<1>(f); + constexpr auto df1 = mito::functions::derivative<1>(f); // a point in space constexpr auto x = mito::geometry::cartesian::coordinates({ 1.0, 1.0 }); diff --git a/tests/mito.lib/manifolds/polar_gradient.cc b/tests/mito.lib/manifolds/polar_gradient.cc index e800ecc93..1301af678 100644 --- a/tests/mito.lib/manifolds/polar_gradient.cc +++ b/tests/mito.lib/manifolds/polar_gradient.cc @@ -36,10 +36,10 @@ TEST(Manifolds, PolarGradient) constexpr auto f = r * mito::functions::sin(theta); // df/dr = sin(theta) - constexpr auto df0 = mito::fields::derivative<0>(f); + constexpr auto df0 = mito::functions::derivative<0>(f); // df/dtheta = r * cos(theta) - constexpr auto df1 = mito::fields::derivative<1>(f); + constexpr auto df1 = mito::functions::derivative<1>(f); // a point in space {r = 2.0, theta = pi / 6.0} constexpr auto x = mito::geometry::polar::coordinates({ 2.0, pi_sixth }); diff --git a/tests/mito.lib/manifolds/spherical_gradient.cc b/tests/mito.lib/manifolds/spherical_gradient.cc index 668854a6e..dc39edfec 100644 --- a/tests/mito.lib/manifolds/spherical_gradient.cc +++ b/tests/mito.lib/manifolds/spherical_gradient.cc @@ -41,13 +41,13 @@ TEST(Manifolds, SphericalGradient) constexpr auto f = r * mito::functions::sin(theta) * mito::functions::cos(phi); // df/dr = sin(theta) * cos(phi) - constexpr auto df0 = mito::fields::derivative<0>(f); + constexpr auto df0 = mito::functions::derivative<0>(f); // df/dtheta = r * cos(theta) * cos(phi) - constexpr auto df1 = mito::fields::derivative<1>(f); + constexpr auto df1 = mito::functions::derivative<1>(f); // df/dphi = - r * cos(theta) * sin(phi) - constexpr auto df2 = mito::fields::derivative<2>(f); + constexpr auto df2 = mito::functions::derivative<2>(f); // a point in space {r = 2.0, theta = pi / 6.0, phi = pi / 6.0} constexpr auto x = mito::geometry::spherical::coordinates({ 2.0, pi_sixth, pi_fourth }); From 309018ea7d89feaf5ff213f86eba2f98eb01a4e8 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 15:41:07 +0100 Subject: [PATCH 04/34] functions: add {determinant} function --- lib/mito/functions/algebra.h | 6 ++++++ lib/mito/functions/function.h | 39 +++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/lib/mito/functions/algebra.h b/lib/mito/functions/algebra.h index 2bc6c72ab..6297e078b 100644 --- a/lib/mito/functions/algebra.h +++ b/lib/mito/functions/algebra.h @@ -97,6 +97,12 @@ namespace mito::functions { return Transpose(f); } + // determinant of f + constexpr auto determinant(const function_c auto & f) + { + return Determinant(f); + } + // f inverse constexpr auto inverse(const function_c auto & f) { diff --git a/lib/mito/functions/function.h b/lib/mito/functions/function.h index 0f7e58c86..0978f153e 100644 --- a/lib/mito/functions/function.h +++ b/lib/mito/functions/function.h @@ -636,6 +636,45 @@ namespace mito::functions { const F _f; }; + // function computing the inverse of a second order tensor + template + requires(tensor::square_matrix_c) + class Determinant : + public Function { + + public: + // the type of the function (what I derive from) + using function_type = + Function; + // the input type + using input_type = function_type::input_type; + // the output type + using output_type = function_type::output_type; + + public: + // constructor + constexpr Determinant(const F & f) : _f(f) {} + + // call operator for function composition + template + constexpr auto operator()(const H & f) const + { + return Composition(*this, f); + } + + // call operator + constexpr auto operator()(const input_type & x) const -> output_type + { + return tensor::determinant(_f(x)); + } + + // the matrix function to invert + constexpr auto f() const -> const F & { return _f; } + + private: + const F _f; + }; + // function computing the inverse of a second order tensor template requires(tensor::square_matrix_c) From 6ae0b7204dd6aa8879f94fc99dbbc6a61cc69662 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 15:41:58 +0100 Subject: [PATCH 05/34] fields: remove {determinant} from {fields} namespace It was redundant with respect to {determinant} from the {functions} namespace. --- lib/mito/fields/fields_algebra.h | 10 ---------- lib/mito/geometry/metric_space.h | 2 +- tests/mito.lib/fields/polar_metric_field.cc | 2 +- tests/mito.lib/fields/spherical_metric_field.cc | 2 +- 4 files changed, 3 insertions(+), 13 deletions(-) diff --git a/lib/mito/fields/fields_algebra.h b/lib/mito/fields/fields_algebra.h index b467fa3af..d9734c32f 100644 --- a/lib/mito/fields/fields_algebra.h +++ b/lib/mito/fields/fields_algebra.h @@ -10,16 +10,6 @@ // Algebraic operations on Fields namespace mito::fields { - // det(f) - template - constexpr auto determinant(const F & f) - { - using coordinates_type = typename F::input_type; - return functions::function([f](const coordinates_type & x) -> mito::tensor::scalar_t { - return determinant(f(x)); - }); - } - // the tensor product of two fields of one-forms template constexpr auto tensor(const F1 & fA, const F2 & fB) diff --git a/lib/mito/geometry/metric_space.h b/lib/mito/geometry/metric_space.h index 250cb8293..b7ef420ce 100644 --- a/lib/mito/geometry/metric_space.h +++ b/lib/mito/geometry/metric_space.h @@ -56,7 +56,7 @@ namespace mito::geometry { public: // the metric volume form static constexpr auto w = - functions::sqrt(fields::determinant(g)) * _wedge(tensor::make_integer_sequence{}); + functions::sqrt(functions::determinant(g)) * _wedge(tensor::make_integer_sequence{}); // get the metric equivalent vector field to a given one-form field static constexpr auto metric_equivalent(const fields::one_form_field_c auto & one_form); diff --git a/tests/mito.lib/fields/polar_metric_field.cc b/tests/mito.lib/fields/polar_metric_field.cc index ca65c7c2c..e906f0a07 100644 --- a/tests/mito.lib/fields/polar_metric_field.cc +++ b/tests/mito.lib/fields/polar_metric_field.cc @@ -53,7 +53,7 @@ TEST(Fields, PolarCoordinates) // the metric volume element constexpr auto w = - mito::functions::sqrt(mito::fields::determinant(g)) * mito::fields::wedge(dr, dt); + mito::functions::sqrt(mito::functions::determinant(g)) * mito::fields::wedge(dr, dt); constexpr auto dr_scalar = mito::tensor::scalar_t{ 0.01 }; constexpr auto dt_scalar = mito::tensor::scalar_t{ 0.01 }; diff --git a/tests/mito.lib/fields/spherical_metric_field.cc b/tests/mito.lib/fields/spherical_metric_field.cc index 8581aac90..1b4a389bf 100644 --- a/tests/mito.lib/fields/spherical_metric_field.cc +++ b/tests/mito.lib/fields/spherical_metric_field.cc @@ -64,7 +64,7 @@ TEST(Manifolds, SphericalCoordinates) // the metric volume element constexpr auto w = - mito::functions::sqrt(mito::fields::determinant(g)) * mito::fields::wedge(dr, dt, dp); + mito::functions::sqrt(mito::functions::determinant(g)) * mito::fields::wedge(dr, dt, dp); constexpr auto dr_scalar = mito::tensor::scalar_t{ 0.01 }; constexpr auto dt_scalar = mito::tensor::scalar_t{ 0.01 }; From 27f34a48cf48b278eada9c9d34c9ae37f488b02a Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 15:57:51 +0100 Subject: [PATCH 06/34] fields: remove {identity_tensor_field} from {fields} namespace Can be replaced by {functions::identity}, no need to bother fields. --- lib/mito/fields/api.h | 3 --- lib/mito/functions/factories.h | 10 ++++++++++ lib/mito/geometry/cartesian/metric.h | 2 +- tests/mito.lib/fields/fields_traits.cc | 2 +- 4 files changed, 12 insertions(+), 5 deletions(-) diff --git a/lib/mito/fields/api.h b/lib/mito/fields/api.h index 4f9e06b91..8b61f60c5 100644 --- a/lib/mito/fields/api.h +++ b/lib/mito/fields/api.h @@ -13,9 +13,6 @@ namespace mito::fields { template constexpr auto uniform_field(const Y & constant); - // the order N identity tensor in D dimensions - template - constexpr auto identity_tensor_field = uniform_field(tensor::identity); } diff --git a/lib/mito/functions/factories.h b/lib/mito/functions/factories.h index c17c358f7..40c1d9347 100644 --- a/lib/mito/functions/factories.h +++ b/lib/mito/functions/factories.h @@ -16,6 +16,16 @@ namespace mito::functions { return Constant(c); } + // returns the order N identity tensor as a function of X + template + constexpr auto identity() + { + // the order N identity tensor + constexpr auto id = tensor::identity; + // return the identity ad a constant function + return constant(id); + } + // returns the function wrapping the functor {f} template constexpr auto function(F && f) diff --git a/lib/mito/geometry/cartesian/metric.h b/lib/mito/geometry/cartesian/metric.h index f27444ffc..f6f9c1703 100644 --- a/lib/mito/geometry/cartesian/metric.h +++ b/lib/mito/geometry/cartesian/metric.h @@ -15,7 +15,7 @@ namespace mito::geometry { static constexpr auto field() { // return the identity field - return fields::identity_tensor_field, D>; + return functions::identity, D>(); } }; } diff --git a/tests/mito.lib/fields/fields_traits.cc b/tests/mito.lib/fields/fields_traits.cc index ca682301f..a9580ec46 100644 --- a/tests/mito.lib/fields/fields_traits.cc +++ b/tests/mito.lib/fields/fields_traits.cc @@ -37,7 +37,7 @@ TEST(Fields, Traits) static_assert(mito::fields::vector_field_c); // a second-order tensor field (2x2 identity tensor field in 3 dimensional space) - constexpr auto i = mito::fields::identity_tensor_field; + constexpr auto i = mito::functions::identity(); // assert that {i} is a tensor field static_assert(mito::fields::tensor_field_c); // assert that {i} is a symmetric tensor field From 5dcc4620cfeae50bf4206c7bd89bd2cdb5ba4889 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 16:13:19 +0100 Subject: [PATCH 07/34] tests: remove {uniform_field} from all tests In most cases, a 'literal constant' can be directly used. In other cases, it can be replaced by {functions::constant}. No need to bother fields. --- tests/mito.lib/fields/calculus_identities.cc | 8 +------- tests/mito.lib/fields/calculus_vector_field.cc | 8 +------- tests/mito.lib/fields/fields_traits.cc | 14 +++++--------- tests/mito.lib/fields/gradient_non_square.cc | 8 +------- tests/mito.lib/fields/polar_metric_field.cc | 8 ++++---- tests/mito.lib/fields/spherical_metric_field.cc | 12 ++++++------ .../integration/quadrature_flip_segment_3D.cc | 4 ++-- .../manifolds/surface_half_sphere_spherical.cc | 2 +- tests/mito.lib/manifolds/triangle_3D.cc | 2 +- tests/mito.lib/quadrature/quadrature_segment_3D.cc | 4 ++-- tests/mito.lib/quadrature/quadrature_square_3D.cc | 2 +- 11 files changed, 25 insertions(+), 47 deletions(-) diff --git a/tests/mito.lib/fields/calculus_identities.cc b/tests/mito.lib/fields/calculus_identities.cc index e395dd7ec..ef0bf7ccd 100644 --- a/tests/mito.lib/fields/calculus_identities.cc +++ b/tests/mito.lib/fields/calculus_identities.cc @@ -23,12 +23,6 @@ constexpr auto pi_fourth = std::numbers::pi / 4.0; TEST(Identities, DivGrad) { - // the field returning the constant {e_0} unit vector in 2D - constexpr auto e0 = mito::fields::uniform_field(e_0); - - // the field returning the constant {e_1} unit vector in 2D - constexpr auto e1 = mito::fields::uniform_field(e_1); - // the sine function constexpr auto sin = mito::functions::sin; @@ -42,7 +36,7 @@ TEST(Identities, DivGrad) constexpr auto x1 = mito::functions::component; // a vector field - constexpr auto f = sin(x0 * x1) * e0 + cos(x0 * x1) * e1; + constexpr auto f = sin(x0 * x1) * e_0 + cos(x0 * x1) * e_1; // a point in space constexpr auto x = mito::geometry::coordinates({ pi_sixth, pi_fourth }); diff --git a/tests/mito.lib/fields/calculus_vector_field.cc b/tests/mito.lib/fields/calculus_vector_field.cc index a31ec88c9..54426f1ef 100644 --- a/tests/mito.lib/fields/calculus_vector_field.cc +++ b/tests/mito.lib/fields/calculus_vector_field.cc @@ -29,12 +29,6 @@ constexpr auto pi_fourth = std::numbers::pi / 4.0; TEST(Laplacian, VectorFields) { - // the field returning the constant {e_0} unit vector in 2D - constexpr auto e0 = mito::fields::uniform_field(e_0); - - // the field returning the constant {e_1} unit vector in 2D - constexpr auto e1 = mito::fields::uniform_field(e_1); - // the sine function constexpr auto sin = mito::functions::sin; @@ -48,7 +42,7 @@ TEST(Laplacian, VectorFields) constexpr auto x1 = mito::functions::component; // a vector field - constexpr auto g = sin(x0 * x1) * e0 + cos(x0 * x1) * e1; + constexpr auto g = sin(x0 * x1) * e_0 + cos(x0 * x1) * e_1; // a point in space constexpr auto x = mito::geometry::coordinates({ pi_sixth, pi_fourth }); diff --git a/tests/mito.lib/fields/fields_traits.cc b/tests/mito.lib/fields/fields_traits.cc index a9580ec46..f30785929 100644 --- a/tests/mito.lib/fields/fields_traits.cc +++ b/tests/mito.lib/fields/fields_traits.cc @@ -10,6 +10,10 @@ // the type of coordinates using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; +// the basis for vectors in 2D +static constexpr auto e_0 = mito::tensor::e_0<2>; +static constexpr auto e_1 = mito::tensor::e_1<2>; + // the function extracting the components of 2D vector static constexpr auto x0 = mito::functions::component; static constexpr auto x1 = mito::functions::component; @@ -23,16 +27,8 @@ TEST(Fields, Traits) // assert that {f} is a scalar field static_assert(mito::fields::scalar_field_c); - // the function returning the constant e0 unit vector in 2D - [[maybe_unused]] constexpr auto e0 = - mito::fields::uniform_field(mito::tensor::e_0<2>); - - // the function returning the constant e1 unit vector in 2D - [[maybe_unused]] constexpr auto e1 = - mito::fields::uniform_field(mito::tensor::e_1<2>); - // a vector field - constexpr auto g = f * (e0 + e1); + constexpr auto g = f * (e_0 + e_1); // assert that {g} is a vector field static_assert(mito::fields::vector_field_c); diff --git a/tests/mito.lib/fields/gradient_non_square.cc b/tests/mito.lib/fields/gradient_non_square.cc index ef37f2845..970a29094 100644 --- a/tests/mito.lib/fields/gradient_non_square.cc +++ b/tests/mito.lib/fields/gradient_non_square.cc @@ -31,12 +31,6 @@ constexpr auto pi_fourth = std::numbers::pi / 4.0; TEST(Gradient, NonSquare) { - // the field returning the constant {e_0} unit vector in 2D - constexpr auto e0 = mito::fields::uniform_field(e_0); - - // the field returning the constant {e_1} unit vector in 2D - constexpr auto e1 = mito::fields::uniform_field(e_1); - // the sine function constexpr auto sin = mito::functions::sin; @@ -53,7 +47,7 @@ TEST(Gradient, NonSquare) constexpr auto x2 = mito::functions::component; // a vector field - constexpr auto f = (sin(x0 * x1) + x2) * e0 + (cos(x0 * x1) - x2) * e1; + constexpr auto f = (sin(x0 * x1) + x2) * e_0 + (cos(x0 * x1) - x2) * e_1; // a point in space constexpr auto x = mito::geometry::coordinates({ pi_sixth, pi_fourth, 1.0 }); diff --git a/tests/mito.lib/fields/polar_metric_field.cc b/tests/mito.lib/fields/polar_metric_field.cc index e906f0a07..2d5992f19 100644 --- a/tests/mito.lib/fields/polar_metric_field.cc +++ b/tests/mito.lib/fields/polar_metric_field.cc @@ -20,12 +20,12 @@ static constexpr auto r = mito::functions::component; TEST(Fields, PolarCoordinates) { // the basis for vector fields (e_r and e_theta) - constexpr auto e_r = mito::fields::uniform_field(mito::tensor::e_0<2>); - constexpr auto e_t = r * mito::fields::uniform_field(mito::tensor::e_1<2>); + constexpr auto e_r = mito::functions::constant(mito::tensor::e_0<2>); + constexpr auto e_t = r * mito::functions::constant(mito::tensor::e_1<2>); // the basis for diagonal second-order tensor fields (e_rr and e_thetatheta) - constexpr auto e_rr = mito::fields::uniform_field(mito::tensor::e_00<2>); - constexpr auto e_tt = mito::fields::uniform_field(mito::tensor::e_11<2>); + constexpr auto e_rr = mito::functions::constant(mito::tensor::e_00<2>); + constexpr auto e_tt = mito::functions::constant(mito::tensor::e_11<2>); // the metric field constexpr auto g = (e_r * e_r) * e_rr + (e_t * e_t) * e_tt; diff --git a/tests/mito.lib/fields/spherical_metric_field.cc b/tests/mito.lib/fields/spherical_metric_field.cc index 1b4a389bf..9bb762d83 100644 --- a/tests/mito.lib/fields/spherical_metric_field.cc +++ b/tests/mito.lib/fields/spherical_metric_field.cc @@ -21,15 +21,15 @@ static constexpr auto t = mito::functions::component; TEST(Manifolds, SphericalCoordinates) { // the basis for vector fields (e_r, e_theta, e_phi) - constexpr auto e_r = mito::fields::uniform_field(mito::tensor::e_0<3>); - constexpr auto e_t = r * mito::fields::uniform_field(mito::tensor::e_1<3>); + constexpr auto e_r = mito::functions::constant(mito::tensor::e_0<3>); + constexpr auto e_t = r * mito::functions::constant(mito::tensor::e_1<3>); constexpr auto e_p = r * mito::functions::sin(t) - * mito::fields::uniform_field(mito::tensor::e_2<3>); + * mito::functions::constant(mito::tensor::e_2<3>); // the basis for diagonal second-order tensor fields (e_rr, e_thetatheta, e_phiphi) - constexpr auto e_rr = mito::fields::uniform_field(mito::tensor::e_00<3>); - constexpr auto e_tt = mito::fields::uniform_field(mito::tensor::e_11<3>); - constexpr auto e_pp = mito::fields::uniform_field(mito::tensor::e_22<3>); + constexpr auto e_rr = mito::functions::constant(mito::tensor::e_00<3>); + constexpr auto e_tt = mito::functions::constant(mito::tensor::e_11<3>); + constexpr auto e_pp = mito::functions::constant(mito::tensor::e_22<3>); // the metric field constexpr auto g = (e_r * e_r) * e_rr + (e_t * e_t) * e_tt + (e_p * e_p) * e_pp; diff --git a/tests/mito.lib/integration/quadrature_flip_segment_3D.cc b/tests/mito.lib/integration/quadrature_flip_segment_3D.cc index af5185ff0..6de6bdd25 100644 --- a/tests/mito.lib/integration/quadrature_flip_segment_3D.cc +++ b/tests/mito.lib/integration/quadrature_flip_segment_3D.cc @@ -39,8 +39,8 @@ TEST(Quadrature, FlipSegment) constexpr auto v_3 = cross / mito::tensor::norm(cross); // the normal and binormal fields to the segment - constexpr auto normal_field_1 = mito::fields::uniform_field(v_2); - constexpr auto normal_field_2 = mito::fields::uniform_field(v_3); + constexpr auto normal_field_1 = mito::functions::constant(v_2); + constexpr auto normal_field_2 = mito::functions::constant(v_3); // the integrand auto f = diff --git a/tests/mito.lib/manifolds/surface_half_sphere_spherical.cc b/tests/mito.lib/manifolds/surface_half_sphere_spherical.cc index 1ab3e327f..c15bb6fd6 100644 --- a/tests/mito.lib/manifolds/surface_half_sphere_spherical.cc +++ b/tests/mito.lib/manifolds/surface_half_sphere_spherical.cc @@ -19,7 +19,7 @@ using cartesian_coordinates_t = mito::geometry::coordinates_t<3, CARTESIAN>; // the basis for vector fields static constexpr auto e_r = - mito::fields::uniform_field(mito::tensor::e_0<3>); + mito::functions::constant(mito::tensor::e_0<3>); TEST(Manifolds, HalfSphereSpherical) diff --git a/tests/mito.lib/manifolds/triangle_3D.cc b/tests/mito.lib/manifolds/triangle_3D.cc index c53ab9dc0..067e43e56 100644 --- a/tests/mito.lib/manifolds/triangle_3D.cc +++ b/tests/mito.lib/manifolds/triangle_3D.cc @@ -27,7 +27,7 @@ TEST(Manifolds, Triangle3D) // the normal vector to the submanifold constexpr auto cross = mito::tensor::cross(x_1 - x_0, x_2 - x_0); constexpr auto normal_vector = cross / mito::tensor::norm(cross); - constexpr auto normal_field = mito::fields::uniform_field(normal_vector); + constexpr auto normal_field = mito::functions::constant(normal_vector); // create a submanifold on {mesh} with the appropriate normal field auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field); diff --git a/tests/mito.lib/quadrature/quadrature_segment_3D.cc b/tests/mito.lib/quadrature/quadrature_segment_3D.cc index 3559d84ec..3702d2bcf 100644 --- a/tests/mito.lib/quadrature/quadrature_segment_3D.cc +++ b/tests/mito.lib/quadrature/quadrature_segment_3D.cc @@ -38,8 +38,8 @@ TEST(Quadrature, Segment3D) constexpr auto cross = mito::tensor::cross(v_1, v_2); constexpr auto v_3 = cross / mito::tensor::norm(cross); - constexpr auto normal_field_1 = mito::fields::uniform_field(v_2); - constexpr auto normal_field_2 = mito::fields::uniform_field(v_3); + constexpr auto normal_field_1 = mito::functions::constant(v_2); + constexpr auto normal_field_2 = mito::functions::constant(v_3); // the integrand auto f = mito::functions::function([](const coordinates_t & x) { return x[0] * x[1]; }); diff --git a/tests/mito.lib/quadrature/quadrature_square_3D.cc b/tests/mito.lib/quadrature/quadrature_square_3D.cc index aeff27399..273d696f0 100644 --- a/tests/mito.lib/quadrature/quadrature_square_3D.cc +++ b/tests/mito.lib/quadrature/quadrature_square_3D.cc @@ -68,7 +68,7 @@ TEST(Quadrature, Square) // the normal vector to the square constexpr auto cross = mito::tensor::cross(x_1 - x_0, x_2 - x_0); constexpr auto normal_vector = cross / mito::tensor::norm(cross); - constexpr auto normal_field = mito::fields::uniform_field(normal_vector); + constexpr auto normal_field = mito::functions::constant(normal_vector); // create a submanifold on {mesh} with the appropriate normal fields auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field); From ec9d5122431c63a4cf685831cd79fb10e07ce897 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 16:16:41 +0100 Subject: [PATCH 08/34] fields: remove {uniform_field} from {fields} namespace Can be replaced by {functions::constant}, no need to bother fields. --- lib/mito/fields/api.h | 4 ---- lib/mito/fields/differential.h | 18 +++++++++++------- lib/mito/fields/factories.h | 6 ------ lib/mito/geometry/cartesian/basis.h | 2 +- lib/mito/geometry/polar/basis.h | 4 ++-- lib/mito/geometry/polar/metric.h | 4 ++-- lib/mito/geometry/spherical/basis.h | 6 +++--- lib/mito/geometry/spherical/metric.h | 6 +++--- 8 files changed, 22 insertions(+), 28 deletions(-) diff --git a/lib/mito/fields/api.h b/lib/mito/fields/api.h index 8b61f60c5..3c0234bcb 100644 --- a/lib/mito/fields/api.h +++ b/lib/mito/fields/api.h @@ -9,10 +9,6 @@ namespace mito::fields { - // uniform field - template - constexpr auto uniform_field(const Y & constant); - } diff --git a/lib/mito/fields/differential.h b/lib/mito/fields/differential.h index b83c2ec20..9720ce8fc 100644 --- a/lib/mito/fields/differential.h +++ b/lib/mito/fields/differential.h @@ -23,7 +23,8 @@ namespace mito::fields { constexpr auto _grad = [](const F & field, std::index_sequence) { // the vector of the partial derivatives return ( - (functions::derivative(field) * uniform_field(tensor::e)) + (functions::derivative(field) + * functions::constant(tensor::e)) + ...); }; @@ -48,8 +49,10 @@ namespace mito::fields { // the tensor of the partial derivatives // (\partial field_I / \partial x_K ) * e_IK, I = 0, ..., N-1, K = 0, ..., D-1 return ( - (functions::derivative(field * uniform_field(tensor::e)) - * uniform_field(tensor::unit, I, K>)) + (functions::derivative( + field * functions::constant(tensor::e)) + * functions::constant( + tensor::unit, I, K>)) + ...); }; @@ -74,7 +77,8 @@ namespace mito::fields { constexpr auto _div = [](const F & field, std::index_sequence) { // the summation of the I-th partial derivative of the I-th component return ( - (functions::derivative(field * uniform_field(tensor::e))) + (functions::derivative( + field * functions::constant(tensor::e))) + ...); }; @@ -99,9 +103,9 @@ namespace mito::fields { // (\partial field_KI / \partial x_I) * e_K return ( (functions::derivative( - (field * uniform_field(tensor::e)) - * uniform_field(tensor::e)) - * uniform_field(tensor::e)) + (field * functions::constant(tensor::e)) + * functions::constant(tensor::e)) + * functions::constant(tensor::e)) + ...); }; diff --git a/lib/mito/fields/factories.h b/lib/mito/fields/factories.h index 62cce8fed..3c0234bcb 100644 --- a/lib/mito/fields/factories.h +++ b/lib/mito/fields/factories.h @@ -9,12 +9,6 @@ namespace mito::fields { - // uniform field - template - constexpr auto uniform_field(const Y & constant) - { - return mito::functions::constant(constant); - } } diff --git a/lib/mito/geometry/cartesian/basis.h b/lib/mito/geometry/cartesian/basis.h index a4812ceb9..dd4079811 100644 --- a/lib/mito/geometry/cartesian/basis.h +++ b/lib/mito/geometry/cartesian/basis.h @@ -19,7 +19,7 @@ namespace mito::geometry { requires(I >= 0 && I < D) static constexpr auto e() { - return fields::uniform_field>(tensor::e); + return functions::constant>(tensor::e); } }; } diff --git a/lib/mito/geometry/polar/basis.h b/lib/mito/geometry/polar/basis.h index cc31679b5..a3d88622b 100644 --- a/lib/mito/geometry/polar/basis.h +++ b/lib/mito/geometry/polar/basis.h @@ -25,12 +25,12 @@ namespace mito::geometry { { if constexpr (I == 0) { // return e_r - return fields::uniform_field(mito::tensor::e_0<2>); + return functions::constant(mito::tensor::e_0<2>); } if constexpr (I == 1) { // e_theta - return _r * fields::uniform_field(mito::tensor::e_1<2>); + return _r * functions::constant(mito::tensor::e_1<2>); } } }; diff --git a/lib/mito/geometry/polar/metric.h b/lib/mito/geometry/polar/metric.h index 25ed02600..3d96fdfc5 100644 --- a/lib/mito/geometry/polar/metric.h +++ b/lib/mito/geometry/polar/metric.h @@ -18,9 +18,9 @@ namespace mito::geometry { constexpr auto r = functions::component; // the function returning the constant e_00 tensor in 2D - constexpr auto e_rr = fields::uniform_field(tensor::e_00<2>); + constexpr auto e_rr = functions::constant(tensor::e_00<2>); // the function returning the constant e_11 tensor in 2D - constexpr auto e_tt = fields::uniform_field(tensor::e_11<2>); + constexpr auto e_tt = functions::constant(tensor::e_11<2>); // return the metric field in polar coordinates return e_rr + functions::pow<2>(r) * e_tt; diff --git a/lib/mito/geometry/spherical/basis.h b/lib/mito/geometry/spherical/basis.h index 6300ada64..fd4743b8f 100644 --- a/lib/mito/geometry/spherical/basis.h +++ b/lib/mito/geometry/spherical/basis.h @@ -28,18 +28,18 @@ namespace mito::geometry { { if constexpr (I == 0) { // return e_r - return fields::uniform_field(mito::tensor::e_0<3>); + return functions::constant(mito::tensor::e_0<3>); } if constexpr (I == 1) { // return e_theta - return _r * fields::uniform_field(mito::tensor::e_1<3>); + return _r * functions::constant(mito::tensor::e_1<3>); } if constexpr (I == 2) { // return e_phi return _r * functions::sin(_theta) - * fields::uniform_field(mito::tensor::e_2<3>); + * functions::constant(mito::tensor::e_2<3>); } } }; diff --git a/lib/mito/geometry/spherical/metric.h b/lib/mito/geometry/spherical/metric.h index 1bf4b22eb..72c959fc8 100644 --- a/lib/mito/geometry/spherical/metric.h +++ b/lib/mito/geometry/spherical/metric.h @@ -21,11 +21,11 @@ namespace mito::geometry { constexpr auto t = functions::component; // the function returning the constant e_00 tensor in 3D - constexpr auto e_rr = fields::uniform_field(tensor::e_00<3>); + constexpr auto e_rr = functions::constant(tensor::e_00<3>); // the function returning the constant e_11 tensor in 3D - constexpr auto e_tt = fields::uniform_field(tensor::e_11<3>); + constexpr auto e_tt = functions::constant(tensor::e_11<3>); // the function returning the constant e_11 tensor in 3D - constexpr auto e_pp = fields::uniform_field(tensor::e_22<3>); + constexpr auto e_pp = functions::constant(tensor::e_22<3>); // return the metric field in spherical coordinates return e_rr + functions::pow<2>(r) * e_tt From 03ec84ed29faecebe0e060c69cb51b9ca977bd2c Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 1 Dec 2025 17:31:45 +0100 Subject: [PATCH 09/34] tests: replace test of quadrature field with equivalent test of point field The test is equivalent but now easier to read. Also, quadrature fields will be removed eventually. --- .cmake/mito_tests_mito_lib.cmake | 2 +- tests/mito.lib/discrete/point_field.cc | 64 ++++++++++++++++++++ tests/mito.lib/discrete/quadrature_field.cc | 67 --------------------- 3 files changed, 65 insertions(+), 68 deletions(-) create mode 100644 tests/mito.lib/discrete/point_field.cc delete mode 100644 tests/mito.lib/discrete/quadrature_field.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 88ab480e6..8ac154ed2 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -41,7 +41,7 @@ mito_test_driver(tests/mito.lib/geometry/spherical_metric_space.cc) mito_test_driver(tests/mito.lib/constraints/dirichlet.cc) # discrete -mito_test_driver(tests/mito.lib/discrete/quadrature_field.cc) +mito_test_driver(tests/mito.lib/discrete/point_field.cc) mito_test_driver(tests/mito.lib/discrete/mesh_field.cc) # fem diff --git a/tests/mito.lib/discrete/point_field.cc b/tests/mito.lib/discrete/point_field.cc new file mode 100644 index 000000000..5e3bd75fd --- /dev/null +++ b/tests/mito.lib/discrete/point_field.cc @@ -0,0 +1,64 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// pick a dimension in space +constexpr auto D = 3; +// a type for the field +using vector_type = mito::tensor::vector_t<2>; +// assemble the quadrature field +using point_field_type = mito::discrete::point_field_t; + +// fetch the point cloud +auto & cloud = mito::geometry::point_cloud(); +// build a point +auto point = cloud.point(); + + +void +constFunction(const point_field_type & field) +{ + // get a (deep) copy of the vector field at {point} + auto vector = field(point); + + // edit the copy + vector[1] = 10; + + // all done + return; +} + +TEST(Discretization, QuadratureFields) +{ + // a quadrature field + auto field = point_field_type("field"); + + // insert one entry using the default value + field.insert(point); + + // get a reference the value of the vector field at {point} + auto & vector = field(point); + + // expect that the vector size is 2 + EXPECT_TRUE(vector.size == 2); + + // modify the content of {vector} + vector[0] = 0.0; + vector[1] = 1.0; + + // expect that editing {vector} has edited the {field} + EXPECT_DOUBLE_EQ(0.0, field(point)[0]); + EXPECT_DOUBLE_EQ(1.0, field(point)[1]); + + // call a function that modifies a (deep) copy of the field + constFunction(field); + + // expect that the {field} is unchanged (deep copy vs. shallow copy) + EXPECT_DOUBLE_EQ(0.0, field(point)[0]); + EXPECT_DOUBLE_EQ(1.0, field(point)[1]); +} \ No newline at end of file diff --git a/tests/mito.lib/discrete/quadrature_field.cc b/tests/mito.lib/discrete/quadrature_field.cc deleted file mode 100644 index d76880dc7..000000000 --- a/tests/mito.lib/discrete/quadrature_field.cc +++ /dev/null @@ -1,67 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -#include -#include - - -// pick a cell type -using cell_type = mito::topology::oriented_simplex_t<1>; -// one quadrature points per cell -constexpr int Q = 1; -// a type for the field -using vector_type = mito::tensor::vector_t<2>; -// assemble the quadrature field -using quadrature_field_type = mito::discrete::quadrature_field_t; - -// an empty topology -auto & topology = mito::topology::topology(); - -// build a segment -const auto segment0 = topology.segment({ topology.vertex(), topology.vertex() }); - - -void -constFunction(const quadrature_field_type & field) -{ - // get a (deep) copy of the vector at (0) - auto vector = field(segment0)[0]; - - // edit the copy - vector[1] = 10; - - // all done - return; -} - -TEST(Discretization, QuadratureFields) -{ - // a quadrature field - auto field = quadrature_field_type("field"); - - // insert one entry using the default value - field.insert(segment0); - - // get a reference the 0-th vector in the field - auto & vector = field(segment0)[0]; - - // expect that the vector size is 2 - EXPECT_TRUE(vector.size == 2); - - // modify the content of {vector} - vector[0] = 0.0; - vector[1] = 1.0; - - // expect that editing {vector} has edited the {field} - EXPECT_DOUBLE_EQ(0.0, field(segment0)[0][0]); - EXPECT_DOUBLE_EQ(1.0, field(segment0)[0][1]); - - // call a function that modifies a (deep) copy of the field - constFunction(field); - - // expect that the {field} is unchanged (deep copy vs. shallow copy) - EXPECT_DOUBLE_EQ(0.0, field(segment0)[0][0]); - EXPECT_DOUBLE_EQ(1.0, field(segment0)[0][1]); -} \ No newline at end of file From 301dded3527ad500115c9be40bfcb8209b77cc11 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 21:16:09 +0100 Subject: [PATCH 10/34] functions: remove {functions::x} Better to use the more explanatory {component} equivalent. --- lib/mito/functions/api.h | 4 -- tests/mito.lib/functions/algebra.cc | 39 ++++++++++++------- .../mito.lib/functions/derivative_product.cc | 26 ++++++++----- tests/mito.lib/functions/derivative_sum.cc | 9 +++-- .../mito.lib/functions/partial_derivatives.cc | 18 ++++++--- .../mito.lib/functions/tensor_derivatives.cc | 9 +++-- 6 files changed, 65 insertions(+), 40 deletions(-) diff --git a/lib/mito/functions/api.h b/lib/mito/functions/api.h index bb0162b8d..869e74c99 100644 --- a/lib/mito/functions/api.h +++ b/lib/mito/functions/api.h @@ -67,10 +67,6 @@ namespace mito::functions { template [[maybe_unused]] constexpr auto component = Component(); - // the function associating to a D-dimensional vector its N-th component - template - [[maybe_unused]] constexpr auto x = component, N>; - // the linear combination template constexpr auto linear_combination(std::array coeffs, Funcs... funcs); diff --git a/tests/mito.lib/functions/algebra.cc b/tests/mito.lib/functions/algebra.cc index 07f75a95a..ce173b86d 100644 --- a/tests/mito.lib/functions/algebra.cc +++ b/tests/mito.lib/functions/algebra.cc @@ -17,20 +17,23 @@ using std::numbers::pi; TEST(Algebra, ScalarValuedFunctions) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a vector - constexpr mito::tensor::vector_t<2> x = { 1.0, pi }; + constexpr vector_t x = { 1.0, pi }; // the function extracting the x_0 component - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; // the function extracting the x_1 component - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // a scalar-valued function constexpr auto function1 = mito::functions::cos(x0 * x1); // a scalar-valued function - constexpr auto function2 = 5.0 * mito::functions::one>; + constexpr auto function2 = 5.0 * mito::functions::one; // the sum of the two scalar functions constexpr auto function3 = function1 + function2; @@ -61,9 +64,9 @@ TEST(Algebra, ScalarValuedFunctions) constexpr auto function12 = (function1 + function2) / function1; static_assert((function1(x) + function2(x)) / function1(x) == function12(x)); - constexpr auto e0 = mito::functions::constant>(mito::tensor::e_0<3>); - constexpr auto e1 = mito::functions::constant>(mito::tensor::e_1<3>); - constexpr auto e2 = mito::functions::constant>(mito::tensor::e_2<3>); + constexpr auto e0 = mito::functions::constant(mito::tensor::e_0<3>); + constexpr auto e1 = mito::functions::constant(mito::tensor::e_1<3>); + constexpr auto e2 = mito::functions::constant(mito::tensor::e_2<3>); // a vector function constexpr auto function16 = function1 * e0 + function1 * e1 + function1 * e2; @@ -82,12 +85,16 @@ TEST(Algebra, ScalarValuedFunctions) TEST(Algebra, TensorValuedFunctions) { + // vectors in 3D + using vector_t = mito::tensor::vector_t<3>; + // the function extracting the x_0 component of a 3D vector - constexpr auto x0 = mito::functions::x<0, 3>; + constexpr auto x0 = mito::functions::component; // the function extracting the x_1 component of a 3D vector - constexpr auto x1 = mito::functions::x<1, 3>; + constexpr auto x1 = mito::functions::component; // the function extracting the x_2 component of a 3D vector - constexpr auto x2 = mito::functions::x<2, 3>; + constexpr auto x2 = mito::functions::component; + // the cosine function constexpr auto cos = mito::functions::cos; @@ -122,19 +129,23 @@ TEST(Algebra, TensorValuedFunctions) TEST(Algebra, ScalarProduct) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a 2D vector - constexpr auto x = mito::tensor::vector_t<2>{ 0.1, 1.0 }; + constexpr auto x = vector_t{ 0.1, 1.0 }; // the function returning the constant e0 unit vector in 2D - constexpr auto e0 = mito::functions::constant>(mito::tensor::e_0<2>); + constexpr auto e0 = mito::functions::constant(mito::tensor::e_0<2>); // check result static_assert(1.0 == (e0 * e0)(x)); // the function extracting the x_0 component - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; + // the function extracting the x_1 component - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // check result static_assert(0.1 == (x0 * x1)(x)); diff --git a/tests/mito.lib/functions/derivative_product.cc b/tests/mito.lib/functions/derivative_product.cc index 643879e8e..479b671d5 100644 --- a/tests/mito.lib/functions/derivative_product.cc +++ b/tests/mito.lib/functions/derivative_product.cc @@ -101,23 +101,26 @@ TEST(Derivatives, Power) TEST(Derivatives, ScalarProduct) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a 2D vector - constexpr auto x = mito::tensor::vector_t<2>{ 0.1, 1.0 }; + constexpr auto x = vector_t{ 0.1, 1.0 }; // a 2D vector - constexpr auto a = mito::tensor::vector_t<2>{ -1.0, 1.0 }; + constexpr auto a = vector_t{ -1.0, 1.0 }; // the function returning the constant e0 unit vector in 2D - constexpr auto e0 = mito::functions::constant>(mito::tensor::e_0<2>); + constexpr auto e0 = mito::functions::constant(mito::tensor::e_0<2>); // the function returning the constant e1 unit vector in 2D - constexpr auto e1 = mito::functions::constant>(mito::tensor::e_1<2>); + constexpr auto e1 = mito::functions::constant(mito::tensor::e_1<2>); // the function extracting the x_0 component of a 2D vector - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; // the function extracting the x_1 component of a 2D vector - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // a function {f} constexpr auto f = a * (x0 * e0 + x1 * e1); @@ -134,11 +137,14 @@ TEST(Derivatives, ScalarProduct) TEST(Derivatives, ScalarProductConstant) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a 2D vector - constexpr auto x = mito::tensor::vector_t<2>{ 0.1, 1.0 }; + constexpr auto x = vector_t{ 0.1, 1.0 }; // a 2D vector - constexpr auto a = mito::tensor::vector_t<2>{ -1.0, 1.0 }; + constexpr auto a = vector_t{ -1.0, 1.0 }; // the constant e0 unit vector in 2D constexpr auto e0 = mito::tensor::e_0<2>; @@ -147,10 +153,10 @@ TEST(Derivatives, ScalarProductConstant) constexpr auto e1 = mito::tensor::e_1<2>; // the function extracting the x_0 component of a 2D vector - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; // the function extracting the x_1 component of a 2D vector - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // a function {f} constexpr auto f = a * (x0 * e0 + x1 * e1); diff --git a/tests/mito.lib/functions/derivative_sum.cc b/tests/mito.lib/functions/derivative_sum.cc index b2fbd9b3e..741452dc4 100644 --- a/tests/mito.lib/functions/derivative_sum.cc +++ b/tests/mito.lib/functions/derivative_sum.cc @@ -83,13 +83,16 @@ TEST(Derivatives, Subtraction) TEST(Derivatives, VectorSum) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a 2D vector - constexpr auto x = mito::tensor::vector_t<2>{ 0.1, 1.0 }; + constexpr auto x = vector_t{ 0.1, 1.0 }; // the function extracting the x_0 component - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; // the function extracting the x_1 component - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // check result static_assert(1.1 == (x0 + x1)(x)); diff --git a/tests/mito.lib/functions/partial_derivatives.cc b/tests/mito.lib/functions/partial_derivatives.cc index 2fe7fa353..2236771e3 100644 --- a/tests/mito.lib/functions/partial_derivatives.cc +++ b/tests/mito.lib/functions/partial_derivatives.cc @@ -9,16 +9,19 @@ TEST(VectorFunctions, Components) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a 2D vector - constexpr auto x = mito::tensor::vector_t<2>{ 0.1, 1.0 }; + constexpr auto x = vector_t{ 0.1, 1.0 }; // the function extracting the x_0 component - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; // check result static_assert(0.1 == x0(x)); // the function extracting the x_1 component - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // check result static_assert(1.0 == x1(x)); @@ -56,13 +59,16 @@ TEST(VectorFunctions, Components) TEST(Derivatives, PartialDerivatives) { + // vectors in 2D + using vector_t = mito::tensor::vector_t<2>; + // a 2D vector - constexpr auto x = mito::tensor::vector_t<2>{ 0.1, 1.0 }; + constexpr auto x = vector_t{ 0.1, 1.0 }; // the function extracting the x_0 component - constexpr auto x0 = mito::functions::x<0, 2>; + constexpr auto x0 = mito::functions::component; // the function extracting the x_1 component - constexpr auto x1 = mito::functions::x<1, 2>; + constexpr auto x1 = mito::functions::component; // sin(x0 * x1) constexpr auto sin = mito::functions::sin(x0 * x1); diff --git a/tests/mito.lib/functions/tensor_derivatives.cc b/tests/mito.lib/functions/tensor_derivatives.cc index 0ce8991c2..23bcc139b 100644 --- a/tests/mito.lib/functions/tensor_derivatives.cc +++ b/tests/mito.lib/functions/tensor_derivatives.cc @@ -7,9 +7,12 @@ #include -// the coordinates functions in 2D space -static constexpr auto x0 = mito::functions::x<0, 2>; -static constexpr auto x1 = mito::functions::x<1, 2>; +// vectors in 2D +using vector_t = mito::tensor::vector_t<2>; + +// the components of vectors in 2D space +static constexpr auto x0 = mito::functions::component; +static constexpr auto x1 = mito::functions::component; // unit vectors in 3D static constexpr auto e0 = mito::tensor::e_0<3>; From d7059abb75a60b6dd17ee296928a923764ef77fc Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 30 Nov 2025 21:17:36 +0100 Subject: [PATCH 11/34] functions: make {Component} available only for subscriptable types --- lib/mito/functions/api.h | 2 +- lib/mito/functions/library_functions.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/mito/functions/api.h b/lib/mito/functions/api.h index 869e74c99..024545570 100644 --- a/lib/mito/functions/api.h +++ b/lib/mito/functions/api.h @@ -64,7 +64,7 @@ namespace mito::functions { constexpr auto function(F && f); // the function extracting the N-th component from an input argument of type {T} - template + template [[maybe_unused]] constexpr auto component = Component(); // the linear combination diff --git a/lib/mito/functions/library_functions.h b/lib/mito/functions/library_functions.h index bd92f4a2a..ceff68b37 100644 --- a/lib/mito/functions/library_functions.h +++ b/lib/mito/functions/library_functions.h @@ -13,7 +13,7 @@ namespace mito::functions { // function extracting the I-th component of a tensor - template + template class Component : public Function { public: From 2c976b0e6182e3cacc028c9ea0e5145ddd358a42 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Fri, 5 Dec 2025 09:17:17 +0100 Subject: [PATCH 12/34] tests/discrete: fix typos --- tests/mito.lib/discrete/point_field.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/mito.lib/discrete/point_field.cc b/tests/mito.lib/discrete/point_field.cc index 5e3bd75fd..ba3de1e8d 100644 --- a/tests/mito.lib/discrete/point_field.cc +++ b/tests/mito.lib/discrete/point_field.cc @@ -10,9 +10,9 @@ // pick a dimension in space constexpr auto D = 3; // a type for the field -using vector_type = mito::tensor::vector_t<2>; +using vector_t = mito::tensor::vector_t<2>; // assemble the quadrature field -using point_field_type = mito::discrete::point_field_t; +using point_field_t = mito::discrete::point_field_t; // fetch the point cloud auto & cloud = mito::geometry::point_cloud(); @@ -21,7 +21,7 @@ auto point = cloud.point(); void -constFunction(const point_field_type & field) +constFunction(const point_field_t & field) { // get a (deep) copy of the vector field at {point} auto vector = field(point); @@ -33,10 +33,10 @@ constFunction(const point_field_type & field) return; } -TEST(Discretization, QuadratureFields) +TEST(Discretization, PointFields) { - // a quadrature field - auto field = point_field_type("field"); + // a point field + auto field = point_field_t("field"); // insert one entry using the default value field.insert(point); From a27bc1a95b93b9c0cff0999adcdd74aa5a4c524f Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Fri, 5 Dec 2025 10:38:40 +0100 Subject: [PATCH 13/34] utilities: add utility function to print a type into a string --- lib/mito/utilities/externals.h | 3 +++ lib/mito/utilities/utilities.h | 31 +++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/lib/mito/utilities/externals.h b/lib/mito/utilities/externals.h index f3d4a1e47..aaa6a97a3 100644 --- a/lib/mito/utilities/externals.h +++ b/lib/mito/utilities/externals.h @@ -15,6 +15,9 @@ #include #include #include +#include +#include +#include // support #include "../journal.h" diff --git a/lib/mito/utilities/utilities.h b/lib/mito/utilities/utilities.h index 6f195f998..c7f9a1aad 100644 --- a/lib/mito/utilities/utilities.h +++ b/lib/mito/utilities/utilities.h @@ -23,6 +23,37 @@ namespace mito::utilities { // strip constant/volatile and reference from type template using base_type = typename std::remove_cvref_t; + + // function returning the type name of {T} as a string + template + std::string type_name() + { + std::string result; + + // Add qualifiers manually + if (std::is_const_v) + result += "const "; + if (std::is_volatile_v) + result += "volatile "; + if (std::is_lvalue_reference_v) + result += "& "; + if (std::is_rvalue_reference_v) + result += "&& "; + + using Base = std::remove_cv_t>; + + int status = 0; + char * dem = abi::__cxa_demangle(typeid(Base).name(), nullptr, nullptr, &status); + + if (status == 0) { + result += dem; + } else { + result += typeid(Base).name(); // fallback + } + std::free(dem); + + return result; + } } From 26d853458cf0d6c9f290f69fc57a0c1965a99110 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 7 Dec 2025 09:14:19 +0100 Subject: [PATCH 14/34] discrete: insertion in {DiscreteField} now happens through {emplace} instead of {operator[]} {emplace} is more appropriate here because: 1) it does not need the value type to implement a default constructor; 2) it does not overwrite existing entries. --- lib/mito/discrete/DiscreteField.h | 6 +++++- lib/mito/discrete/externals.h | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/lib/mito/discrete/DiscreteField.h b/lib/mito/discrete/DiscreteField.h index 152b0a5a1..f16080667 100644 --- a/lib/mito/discrete/DiscreteField.h +++ b/lib/mito/discrete/DiscreteField.h @@ -73,8 +73,12 @@ namespace mito::discrete { */ inline auto insert(const key_type & key, const value_type & entry = value_type()) -> void { + // check that the entry is not already present + // (emplace does not overwrite existing entries) + assert(_map_entry_to_values.find(key) == std::cend(_map_entry_to_values)); + // add a new entry to the field - _map_entry_to_values[key] = entry; + _map_entry_to_values.emplace(key, entry); // all done return; diff --git a/lib/mito/discrete/externals.h b/lib/mito/discrete/externals.h index 3f27eafc0..18e7afef1 100644 --- a/lib/mito/discrete/externals.h +++ b/lib/mito/discrete/externals.h @@ -9,6 +9,7 @@ // externals #include +#include // support #include "../mesh.h" From b826646d8921f9e02ae618f672509a8480759b85 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 7 Dec 2025 19:56:49 +0100 Subject: [PATCH 15/34] geometry: implement method {CoordinateSystem::origin()} --- lib/mito/geometry/CoordinateSystem.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lib/mito/geometry/CoordinateSystem.h b/lib/mito/geometry/CoordinateSystem.h index 3e9b3d947..025848a52 100644 --- a/lib/mito/geometry/CoordinateSystem.h +++ b/lib/mito/geometry/CoordinateSystem.h @@ -88,6 +88,9 @@ namespace mito::geometry { return _coordinates_map.at(point); } + // get the zero coordinates + constexpr auto origin() const -> coordinates_type { return coordinates_type{}; } + // get the coordinates of the midpoint between {point_a} and {point_b} constexpr auto midpoint(const point_type & point_a, const point_type & point_b) const -> coordinates_type From d8bb6a28b2e8c3d02b1c8c527a1c863b1432e5d0 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 7 Dec 2025 19:58:07 +0100 Subject: [PATCH 16/34] geometry: add (N+1)-th coordinate to {ReferenceSimplex} This coordinate resolves to 1-xi0-xi1-..., consistently with what expected for barycentric coordinates. --- lib/mito/geometry/ReferenceSimplex.h | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/lib/mito/geometry/ReferenceSimplex.h b/lib/mito/geometry/ReferenceSimplex.h index 4c9c7b403..253f7d092 100644 --- a/lib/mito/geometry/ReferenceSimplex.h +++ b/lib/mito/geometry/ReferenceSimplex.h @@ -18,13 +18,29 @@ namespace mito::geometry { // the area of the reference simplex constexpr static double area = 1.0 / mito::tensor::factorial(); + private: + // helper to compute 1 - xi0 - xi1 - ... - xi(N-1) + template + static constexpr auto _one_minus_xis(tensor::integer_sequence) + { + return (1.0 - ... - xi); + } + public: // the type of coordinates in the parametric space using parametric_coordinates_type = coordinates_t; // the function extracting the I component of a parametric point template - static constexpr auto xi = functions::component; + static constexpr auto xi = [] { + if constexpr (I < N) { + // I-th parametric coordinate + return functions::component; + } else { + // xi = 1 - xi<0> - xi<1> - ... - xi + return _one_minus_xis(tensor::make_integer_sequence{}); + } + }(); }; } // namespace mito From 2bc22791256d8fb13ae88f4f970938f4b4be4a17 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 7 Dec 2025 20:04:12 +0100 Subject: [PATCH 17/34] geometry: compute the {parametrization} of {GeometricSimplex} using barycentric coordinates of the reference simplex --- lib/mito/geometry/GeometricSimplex.h | 49 +++++++++++----------------- lib/mito/quadrature/Integrator.h | 11 +++---- 2 files changed, 24 insertions(+), 36 deletions(-) diff --git a/lib/mito/geometry/GeometricSimplex.h b/lib/mito/geometry/GeometricSimplex.h index 5e6555e51..7432b6270 100644 --- a/lib/mito/geometry/GeometricSimplex.h +++ b/lib/mito/geometry/GeometricSimplex.h @@ -134,38 +134,27 @@ namespace mito::geometry { // return the composition of this simplex in terms of its vertices constexpr auto nodes() const -> const nodes_type & { return _nodes; } - // get the coordinates of a parametric point in physical space + // get the parametrization of the geometric simplex in physical space template - constexpr auto parametrization( - const parametric_coordinates_type & point, - const coordinateSystemT & coordinate_system) const - -> coordinateSystemT::coordinates_type + constexpr auto parametrization(const coordinateSystemT & coordinate_system) const -> auto { - // get the coordinates of the first node - auto coord_0 = coordinate_system.coordinates(_nodes[0]->point()); - - // the vector going from {coord_0} to the position of the parametric point (initialize - // with the zero vector) - auto result = coord_0 - coord_0; - // the sum of the parametric coordinates (initialize with zero) - auto xi = 0.0; - // loop on the vertices of the element - tensor::constexpr_for_1([&]() { - // get the coordinates of the i-th vertex - const auto & coord_i = coordinate_system.coordinates(_nodes[i]->point()); - // assemble the contribution of the i-th vertex - result += (coord_i - coord_0) * point[i]; - // accumulate the parametric coordinate - xi += point[i]; - }); - // get the coordinates of the last vertex - const auto & coord_i = coordinate_system.coordinates(_nodes[n_vertices - 1]->point()); - // add the contribution of the last vertex - // (note that the last barycentric coordinate is 1 - sum of the others) - result += (coord_i - coord_0) * (1.0 - xi); - - // return the coordinates of the parametric point - return coord_0 + result; + // helper to assemble the parametrization on this simplex + constexpr auto _assemble = []( + const auto & nodes, const auto & coordinate_system, + tensor::integer_sequence) { + // get the origin of the coordinate system + constexpr auto origin = coordinate_system.origin(); + + // assemble the parametrization as x0 * xi<0> + ... + // where {xi} are the barycentric coordinates on the reference simplex and the + // {xa} are the position vectors of the nodes + return ( + ((reference_simplex_type::template xi + * (coordinate_system.coordinates(nodes[a]->point()) - origin))) + + ...); + }; + return _assemble( + _nodes, coordinate_system, tensor::make_integer_sequence{}); } private: diff --git a/lib/mito/quadrature/Integrator.h b/lib/mito/quadrature/Integrator.h index 4490763cc..83bc62b79 100644 --- a/lib/mito/quadrature/Integrator.h +++ b/lib/mito/quadrature/Integrator.h @@ -37,13 +37,12 @@ namespace mito::quadrature { { // loop on elements for (const auto & element : _manifold.elements()) { - // use element parametrization and manifold's coordinate systemto map the position - // of quadrature points in the canonical element to the coordinate of the quadrature - // point + // get element parametrization under the manifold's coordinate system + const auto parametrization = element.parametrization(_manifold.coordinate_system()); + // populate the field with the coordinates of the quadrature points in physical + // space _coordinates.insert( - element.simplex(), - { element.parametrization( - _quadratureRule.point(q), _manifold.coordinate_system())... }); + element.simplex(), { parametrization(_quadratureRule.point(q))... }); } // all done From 1164fdbfc3b8cfb2ad5539a7446fdf0aeec154ff Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 7 Dec 2025 20:11:51 +0100 Subject: [PATCH 18/34] fem/elements: isoparametric elements now inherit the parametrization from the underlying {GeometricSimplex} --- lib/mito/fem/elements/IsoparametricTriangle.h | 8 +++++++- .../elements/tri1/IsoparametricTriangleP1.h | 12 ----------- .../elements/tri2/IsoparametricTriangleP2.h | 20 ------------------- 3 files changed, 7 insertions(+), 33 deletions(-) diff --git a/lib/mito/fem/elements/IsoparametricTriangle.h b/lib/mito/fem/elements/IsoparametricTriangle.h index 7c8ed5bcb..94fb54e21 100644 --- a/lib/mito/fem/elements/IsoparametricTriangle.h +++ b/lib/mito/fem/elements/IsoparametricTriangle.h @@ -35,6 +35,7 @@ namespace mito::fem { constexpr IsoparametricTriangle( const cell_type & cell, const coordinate_system_type & coord_system) : _cell(cell), + _coord_system(coord_system), _x0{ coord_system.coordinates(cell.nodes()[0]->point()) - coordinates_type{} }, _x1{ coord_system.coordinates(cell.nodes()[1]->point()) - coordinates_type{} }, _x2{ coord_system.coordinates(cell.nodes()[2]->point()) - coordinates_type{} } @@ -59,11 +60,16 @@ namespace mito::fem { // get the geometric simplex constexpr auto cell() const noexcept -> const cell_type & { return _cell; } + // get the mapping from parametric coordinates to physical coordinates + constexpr auto parametrization() const { return _cell.parametrization(_coord_system); } + protected: - // QUESTION: do we need to maintain a reference to the geometric simplex? // a const reference to the geometric simplex const cell_type & _cell; + // a const reference to the coordinate system + const coordinate_system_type & _coord_system; + // the coordinates of the discretization nodes of the triangle const vector_type _x0; const vector_type _x1; diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index b57fcab3f..ab58d5f29 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -64,18 +64,6 @@ namespace mito::fem { return _connectivity; } - // get the isoparametric mapping from barycentric coordinates to physical coordinates - constexpr auto parametrization() const - { - // get the shape functions - constexpr auto phi_0 = shape_functions.shape<0>(); - constexpr auto phi_1 = shape_functions.shape<1>(); - constexpr auto phi_2 = shape_functions.shape<2>(); - - // return the isoparametric mapping from parametric to physical coordinates - return functions::linear_combination(std::array{ _x0, _x1, _x2 }, phi_0, phi_1, phi_2); - } - // get the shape function associated with local node {a} template requires(a >= 0 && a < n_nodes) diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 2c8cf070e..5db4dd804 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -58,26 +58,6 @@ namespace mito::fem { constexpr IsoparametricTriangleP2 & operator=(IsoparametricTriangleP2 &&) noexcept = delete; public: - // get the isoparametric mapping from barycentric coordinates to physical coordinates - constexpr auto parametrization() const - { - auto x3 = 0.5 * (_x0 + _x1); - auto x4 = 0.5 * (_x1 + _x2); - auto x5 = 0.5 * (_x2 + _x0); - - // get the shape functions - constexpr auto phi_0 = shape_functions.shape<0>(); - constexpr auto phi_1 = shape_functions.shape<1>(); - constexpr auto phi_2 = shape_functions.shape<2>(); - constexpr auto phi_3 = shape_functions.shape<3>(); - constexpr auto phi_4 = shape_functions.shape<4>(); - constexpr auto phi_5 = shape_functions.shape<5>(); - - // return the isoparametric mapping from barycentric to physical coordinates - return functions::linear_combination( - std::array{ _x0, _x1, _x2, x3, x4, x5 }, phi_0, phi_1, phi_2, phi_3, phi_4, phi_5); - } - // get the discretization nodes constexpr auto connectivity() const noexcept -> const connectivity_type & { From e22fe9d769b12f736369015e635af48471ec6a3c Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 10:31:40 +0100 Subject: [PATCH 19/34] fem: implement class {DomainField} Instances of this class can be localized on finite elements via the {localize} method. --- .cmake/mito_tests_mito_lib.cmake | 1 + lib/mito/fem/DomainField.h | 61 ++++++++++++++++++++++++++++++ lib/mito/fem/api.h | 8 ++++ lib/mito/fem/factories.h | 8 ++++ lib/mito/fem/forward.h | 5 +++ lib/mito/fem/public.h | 1 + tests/mito.lib/fem/domain_field.cc | 58 ++++++++++++++++++++++++++++ 7 files changed, 142 insertions(+) create mode 100644 lib/mito/fem/DomainField.h create mode 100644 tests/mito.lib/fem/domain_field.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 8ac154ed2..2f027acac 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -51,6 +51,7 @@ mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_construction.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p2.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc) +mito_test_driver(tests/mito.lib/fem/domain_field.cc) # io mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc) diff --git a/lib/mito/fem/DomainField.h b/lib/mito/fem/DomainField.h new file mode 100644 index 000000000..5d3621feb --- /dev/null +++ b/lib/mito/fem/DomainField.h @@ -0,0 +1,61 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {DomainField} represents a field defined over a domain. It is built from a continuous +// field defined in the physical space. The field can be localized on finite elements via the +// {localize} method, which composes the underlying field with the element parametrization. + + +namespace mito::fem { + + template + class DomainField { + + private: + // the underlying field type + using field_type = F; + + public: + // default constructor + constexpr DomainField(const F & field) : _field(field) {} + + // default destructor + constexpr ~DomainField() = default; + + // delete copy constructor + constexpr DomainField(const DomainField &) = delete; + + // default move constructor + constexpr DomainField(DomainField &&) = delete; + + // delete copy assignment + constexpr auto operator=(const DomainField &) -> DomainField & = delete; + + // delete move assignment + constexpr auto operator=(DomainField &&) -> DomainField & = delete; + + public: + // localize the field on {element} + template + constexpr auto localize(const elementT & element) const -> auto + { + // return the composition of the underlying field with the element parametrization + return _field(element.parametrization()); + } + + private: + // the underlying global field + const field_type _field; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/api.h b/lib/mito/fem/api.h index b0222bb33..fdc9df4c4 100644 --- a/lib/mito/fem/api.h +++ b/lib/mito/fem/api.h @@ -13,6 +13,14 @@ namespace mito::fem { template constexpr auto nodal_field(const functionSpaceT & function_space, std::string name); + // domain field alias + template + using domain_field_t = DomainField; + + // domain field factory + template + constexpr auto domain_field(const F &); + // the possible discretization types: continuous Galerking (CG) vs. discontinuous Galerkin (DG) enum class discretization_t { CG, DG }; diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index 6148588a8..45ee0ac25 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -24,6 +24,14 @@ namespace mito::fem { return discrete::nodal_field_t(nodes, name); } + // domain field factory + template + constexpr auto domain_field(const F & field) + { + // build a domain field from the given function + return domain_field_t(field); + } + // TOFIX: create a constructor that takes no constraints // TOFIX: {constraints} should be a collection of constraints as opposed to an instance of a diff --git a/lib/mito/fem/forward.h b/lib/mito/fem/forward.h index d212e8f4a..aa2e4773d 100644 --- a/lib/mito/fem/forward.h +++ b/lib/mito/fem/forward.h @@ -29,6 +29,11 @@ namespace mito::fem { // class discrete system template class DiscreteSystem; + + // class domain field + template + class DomainField; + } diff --git a/lib/mito/fem/public.h b/lib/mito/fem/public.h index 1ce8000ab..461f024f0 100644 --- a/lib/mito/fem/public.h +++ b/lib/mito/fem/public.h @@ -21,6 +21,7 @@ // classes implementation #include "FunctionSpace.h" +#include "DomainField.h" #include "DiscreteSystem.h" #include "Weakform.h" diff --git a/tests/mito.lib/fem/domain_field.cc b/tests/mito.lib/fem/domain_field.cc new file mode 100644 index 000000000..a4dd1c31d --- /dev/null +++ b/tests/mito.lib/fem/domain_field.cc @@ -0,0 +1,58 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// cartesian coordinates in 2D +using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; +// the type of coordinate system +using coord_system_t = mito::geometry::coordinate_system_t; +// the x scalar field in 2D +constexpr auto x = mito::functions::component; +// the y scalar field in 2D +constexpr auto y = mito::functions::component; + + +TEST(Fem, DomainField) +{ + // create a channel + journal::info_t channel("tests.domain_field"); + + // the coordinate system + auto coord_system = coord_system_t(); + + // create a domain field + auto field = mito::fem::domain_field(x * y); + + // create some nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0, 0.0 }); + auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0 }); + + // create a geometric simplex + auto geometric_simplex = mito::geometry::triangle<2>({ node_0, node_1, node_2 }); + + // an isoparametric triangle + auto element = mito::fem::IsoparametricTriangle(geometric_simplex, coord_system); + + // TOFIX: This syntax should also be allowed: field.localize(geometric_simplex). + // However, in case, the coordinate system should be passed somehow to the function. + // localize the field on the simplex + auto localized_field = field.localize(element); + + // evaluate the localized field at the center of the triangle + auto value = localized_field({ 1.0 / 3.0, 1.0 / 3.0 }); + + // check the value + EXPECT_DOUBLE_EQ(value, (x * y)({ 1.0 / 3.0, 1.0 / 3.0 })); + + // all done + return; +} + + +// end of file \ No newline at end of file From 7262b82d6ac7c7b6cea3cc1e8faa8f435160cb12 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 10:50:16 +0100 Subject: [PATCH 20/34] discrete: rename {key_type} to {input_type} in class {DiscreteField} --- lib/mito/discrete/DiscreteField.h | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/lib/mito/discrete/DiscreteField.h b/lib/mito/discrete/DiscreteField.h index f16080667..e9d653609 100644 --- a/lib/mito/discrete/DiscreteField.h +++ b/lib/mito/discrete/DiscreteField.h @@ -12,15 +12,16 @@ namespace mito::discrete { template class DiscreteField { - private: - // the key type - using key_type = keyT; - // the value type + public: + // the input type of the field + using input_type = keyT; + // the value type returned by the field using value_type = valueT; - // a map from {key_type} to {value_type} values + private: + // a map from {input_type} to {value_type} values using map_type = - std::unordered_map>; + std::unordered_map>; public: // constructor @@ -55,7 +56,7 @@ namespace mito::discrete { /** * accessor for the value of a given entry */ - inline auto operator()(const key_type & key) const -> const value_type & + inline auto operator()(const input_type & key) const -> const value_type & { return _map_entry_to_values.at(key); } @@ -63,7 +64,7 @@ namespace mito::discrete { /** * mutator for the value of a given entry */ - inline auto operator()(const key_type & key) -> value_type & + inline auto operator()(const input_type & key) -> value_type & { return _map_entry_to_values.at(key); } @@ -71,7 +72,7 @@ namespace mito::discrete { /** * insert a new entry to the field */ - inline auto insert(const key_type & key, const value_type & entry = value_type()) -> void + inline auto insert(const input_type & key, const value_type & entry = value_type()) -> void { // check that the entry is not already present // (emplace does not overwrite existing entries) From a18cbf013bb4147a0db2d7984090b69cbb209539 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 10:57:39 +0100 Subject: [PATCH 21/34] fem: move utilities header to parent directory --- lib/mito/fem/elements/public.h | 3 --- lib/mito/fem/public.h | 4 ++++ lib/mito/fem/{elements => }/utilities.h | 0 3 files changed, 4 insertions(+), 3 deletions(-) rename lib/mito/fem/{elements => }/utilities.h (100%) diff --git a/lib/mito/fem/elements/public.h b/lib/mito/fem/elements/public.h index 40f144579..cf056e39c 100644 --- a/lib/mito/fem/elements/public.h +++ b/lib/mito/fem/elements/public.h @@ -26,7 +26,4 @@ // factories implementation #include "factories.h" -// utilities implementation -#include "utilities.h" - // end of file diff --git a/lib/mito/fem/public.h b/lib/mito/fem/public.h index 461f024f0..6ec9a959d 100644 --- a/lib/mito/fem/public.h +++ b/lib/mito/fem/public.h @@ -19,8 +19,12 @@ // blocks implementation #include "blocks.h" +// utilities implementation +#include "utilities.h" + // classes implementation #include "FunctionSpace.h" +#include "FemField.h" #include "DomainField.h" #include "DiscreteSystem.h" #include "Weakform.h" diff --git a/lib/mito/fem/elements/utilities.h b/lib/mito/fem/utilities.h similarity index 100% rename from lib/mito/fem/elements/utilities.h rename to lib/mito/fem/utilities.h From 78fb3a4f20017816c8bb510753f4e398559ef717 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 11:24:20 +0100 Subject: [PATCH 22/34] fem: add class {FemField} --- lib/mito/fem/FemField.h | 93 ++++++++++++++++++++++++++++++++++++++++ lib/mito/fem/api.h | 8 ++++ lib/mito/fem/factories.h | 7 +++ lib/mito/fem/forward.h | 3 ++ 4 files changed, 111 insertions(+) create mode 100644 lib/mito/fem/FemField.h diff --git a/lib/mito/fem/FemField.h b/lib/mito/fem/FemField.h new file mode 100644 index 000000000..0ca83380e --- /dev/null +++ b/lib/mito/fem/FemField.h @@ -0,0 +1,93 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {FemField} represents a finite element field defined via its nodal values on a set of +// discretization nodes. The field can be localized on finite elements via the {localize} method, +// which assembles the field values on the element from the nodal values and the element shape +// functions. + +namespace mito::fem { + + // TODO: implement higher-dimensional fields (e.g. vector fields, tensor fields, ...) + + template + class FemField { + + private: + // the field value type + using field_value_type = fieldValueT; + // the nodal field type + using nodal_field_type = discrete::nodal_field_t; + // the node type + using node_type = typename nodal_field_type::input_type; + + public: + // constructor from temporary nodal field + FemField(nodal_field_type && nodal_field) : _nodal_field(nodal_field) {} + + // default destructor + ~FemField() = default; + + // delete copy constructor + FemField(const FemField &) = delete; + + // default move constructor + FemField(FemField &&) = default; + + // delete copy assignment + auto operator=(const FemField &) -> FemField & = delete; + + // delete move assignment + auto operator=(FemField &&) -> FemField & = delete; + + public: + // get a read-only nodal value at {node} + auto nodal_value(const node_type & node) const -> const field_value_type & + { + return _nodal_field(node); + } + + // get a mutable nodal value at {node} + auto nodal_value(const node_type & node) -> field_value_type & + { + return _nodal_field(node); + } + + // QUESTION: how do we guarantee that {element} is compatible with the function space that + // built this field? + // TODO: add concept that constrains {elementT} to be a finite element type + // localize the function on {element} + template + auto localize(const elementT & element) const -> auto + { + // the element type + using element_type = elementT; + + // helper lambda to assemble the solution on {element} + constexpr auto _assemble = []( + const element_type & element, + const nodal_field_type & field, + tensor::integer_sequence) { + // assemble the solution field from the shape functions + return ((element.template shape() * field(element.connectivity()[a])) + ...); + }; + return _assemble( + element, _nodal_field, tensor::make_integer_sequence{}); + } + + private: + // the field with the nodal values of the field + nodal_field_type _nodal_field; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/api.h b/lib/mito/fem/api.h index fdc9df4c4..f14fcdd40 100644 --- a/lib/mito/fem/api.h +++ b/lib/mito/fem/api.h @@ -21,6 +21,14 @@ namespace mito::fem { template constexpr auto domain_field(const F &); + // finite element field alias + template + using fem_field_t = FemField; + + // fem field factory + template + constexpr auto fem_field(const fieldValueT &); + // the possible discretization types: continuous Galerking (CG) vs. discontinuous Galerkin (DG) enum class discretization_t { CG, DG }; diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index 45ee0ac25..ceddee131 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -32,6 +32,13 @@ namespace mito::fem { return domain_field_t(field); } + // fem field factory + template + constexpr auto fem_field(discrete::nodal_field_t && nodal_field) + { + return fem_field_t(nodal_field); + } + // TOFIX: create a constructor that takes no constraints // TOFIX: {constraints} should be a collection of constraints as opposed to an instance of a diff --git a/lib/mito/fem/forward.h b/lib/mito/fem/forward.h index aa2e4773d..f5c9b3bdb 100644 --- a/lib/mito/fem/forward.h +++ b/lib/mito/fem/forward.h @@ -34,6 +34,9 @@ namespace mito::fem { template class DomainField; + // class finite element field + template + class FemField; } From fe8fa0c39714a6ee9e6549bb2e74dba9c991dba2 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 11:25:30 +0100 Subject: [PATCH 23/34] fem: class {FunctionSpace} is now able to hand out {FemField}s --- lib/mito/fem/FemField.h | 2 +- lib/mito/fem/FunctionSpace.h | 19 ++++++++++++++++++- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/lib/mito/fem/FemField.h b/lib/mito/fem/FemField.h index 0ca83380e..f653ee0cb 100644 --- a/lib/mito/fem/FemField.h +++ b/lib/mito/fem/FemField.h @@ -30,7 +30,7 @@ namespace mito::fem { public: // constructor from temporary nodal field - FemField(nodal_field_type && nodal_field) : _nodal_field(nodal_field) {} + FemField(nodal_field_type && nodal_field) : _nodal_field(std::move(nodal_field)) {} // default destructor ~FemField() = default; diff --git a/lib/mito/fem/FunctionSpace.h b/lib/mito/fem/FunctionSpace.h index f4e7b8246..90a5e734f 100644 --- a/lib/mito/fem/FunctionSpace.h +++ b/lib/mito/fem/FunctionSpace.h @@ -35,7 +35,9 @@ namespace mito::fem { // the type of a map between the mesh nodes and discretization nodes using map_type = std::unordered_map< mesh_node_type, discretization_node_type, utilities::hash_function>; - + // a finite element field type + template + using fem_field_type = fem_field_t; public: // the constructor @@ -90,6 +92,21 @@ namespace mito::fem { // accessor for the node map constexpr auto node_map() const noexcept -> const map_type & { return _node_map; } + // hand out an preallocated fem field for all the discretization nodes with name {name} + template + constexpr auto fem_field(std::string name) const -> fem_field_type + { + // assemble the node type + using node_type = discretization_node_type; + + // get the discretization nodes + std::unordered_set> nodes; + get_discretization_nodes(*this, nodes); + + // build a nodal field on the discretization nodes collected from the function space + return fem_field_t(discrete::nodal_field_t(nodes, name)); + } + private: // a collection of finite elements elements_type _elements; From 34f7721d9df63fe9c4683ab9e291927e6fccb508 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 13:38:09 +0100 Subject: [PATCH 24/34] functions: treat functions defined on coordinates with the same rules as functions of tensors --- lib/mito/functions/externals.h | 1 + lib/mito/functions/forward.h | 6 ++++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/lib/mito/functions/externals.h b/lib/mito/functions/externals.h index 6a73919ff..0356f08f0 100644 --- a/lib/mito/functions/externals.h +++ b/lib/mito/functions/externals.h @@ -12,6 +12,7 @@ // support #include "../journal.h" #include "../tensor.h" +#include "../coordinates.h" // end of file diff --git a/lib/mito/functions/forward.h b/lib/mito/functions/forward.h index 8a47519d2..3c0b66f73 100644 --- a/lib/mito/functions/forward.h +++ b/lib/mito/functions/forward.h @@ -49,9 +49,11 @@ namespace mito::functions { template concept vector_function_c = vector_domain_function_c and vector_valued_function_c; - // concept of a function defined on tensors + // concept of a function defined on tensors or coordinates template - concept tensor_domain_function_c = function_c and tensor::tensor_c; + concept tensor_domain_function_c = function_c + and (tensor::tensor_c + or geometry::coordinates_c); // concept of a tensor-valued function template From 032db7eeea4bb3c00ea305d3698d02b5d3ffca06 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 13:39:17 +0100 Subject: [PATCH 25/34] tests/fem: add gradient calculation to test of domain fields --- tests/mito.lib/fem/domain_field.cc | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tests/mito.lib/fem/domain_field.cc b/tests/mito.lib/fem/domain_field.cc index a4dd1c31d..ddd2564b3 100644 --- a/tests/mito.lib/fem/domain_field.cc +++ b/tests/mito.lib/fem/domain_field.cc @@ -29,9 +29,9 @@ TEST(Fem, DomainField) auto field = mito::fem::domain_field(x * y); // create some nodes - auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); - auto node_1 = mito::geometry::node(coord_system, { 1.0, 0.0 }); - auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0 }); + auto node_0 = mito::geometry::node(coord_system, { 1.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 0.0, 1.0 }); + auto node_2 = mito::geometry::node(coord_system, { 0.0, 0.0 }); // create a geometric simplex auto geometric_simplex = mito::geometry::triangle<2>({ node_0, node_1, node_2 }); @@ -50,6 +50,16 @@ TEST(Fem, DomainField) // check the value EXPECT_DOUBLE_EQ(value, (x * y)({ 1.0 / 3.0, 1.0 / 3.0 })); + // compute the gradient of the localized field with respect to the barycentric coordinates + auto gradient = mito::fields::gradient(localized_field); + + // evaluate the localized field gradient at the center of the triangle + auto value_gradient = gradient({ 1.0 / 3.0, 1.0 / 3.0 }); + + // check the value of the gradient at the center of the triangle + EXPECT_DOUBLE_EQ(value_gradient[0], y({ 1.0 / 3.0, 1.0 / 3.0 })); + EXPECT_DOUBLE_EQ(value_gradient[1], x({ 1.0 / 3.0, 1.0 / 3.0 })); + // all done return; } From b8ffe09b30c668208be93d3d6ba13dcc18405042 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 13:43:13 +0100 Subject: [PATCH 26/34] fem: remove {fem_field} factory Fem fields should be handed out by the function space. --- lib/mito/fem/api.h | 4 ---- lib/mito/fem/factories.h | 7 ------- 2 files changed, 11 deletions(-) diff --git a/lib/mito/fem/api.h b/lib/mito/fem/api.h index f14fcdd40..1ae53c9c3 100644 --- a/lib/mito/fem/api.h +++ b/lib/mito/fem/api.h @@ -25,10 +25,6 @@ namespace mito::fem { template using fem_field_t = FemField; - // fem field factory - template - constexpr auto fem_field(const fieldValueT &); - // the possible discretization types: continuous Galerking (CG) vs. discontinuous Galerkin (DG) enum class discretization_t { CG, DG }; diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index ceddee131..45ee0ac25 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -32,13 +32,6 @@ namespace mito::fem { return domain_field_t(field); } - // fem field factory - template - constexpr auto fem_field(discrete::nodal_field_t && nodal_field) - { - return fem_field_t(nodal_field); - } - // TOFIX: create a constructor that takes no constraints // TOFIX: {constraints} should be a collection of constraints as opposed to an instance of a From 59780d7aa1a68753b7874626d07d1457f0a8ca85 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 13:44:08 +0100 Subject: [PATCH 27/34] fem: rename {nodal_value} as {operator()} in class {FemField} --- lib/mito/fem/FemField.h | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/lib/mito/fem/FemField.h b/lib/mito/fem/FemField.h index f653ee0cb..d4d495b58 100644 --- a/lib/mito/fem/FemField.h +++ b/lib/mito/fem/FemField.h @@ -49,16 +49,13 @@ namespace mito::fem { public: // get a read-only nodal value at {node} - auto nodal_value(const node_type & node) const -> const field_value_type & + auto operator()(const node_type & node) const -> const field_value_type & { return _nodal_field(node); } // get a mutable nodal value at {node} - auto nodal_value(const node_type & node) -> field_value_type & - { - return _nodal_field(node); - } + auto operator()(const node_type & node) -> field_value_type & { return _nodal_field(node); } // QUESTION: how do we guarantee that {element} is compatible with the function space that // built this field? From 106a18c266fb702c3bbd3780e81a5247563f0099 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 13:51:22 +0100 Subject: [PATCH 28/34] fem: class {DiscreteSystem} now stores the solution as a {FemField} instance as opposed to a nodal field --- benchmarks/mito.lib/pdes/poisson.cc | 6 +- lib/mito/fem/DiscreteSystem.h | 95 +++++++++-------------------- 2 files changed, 32 insertions(+), 69 deletions(-) diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc index a1b3e6b00..d716d34c1 100644 --- a/benchmarks/mito.lib/pdes/poisson.cc +++ b/benchmarks/mito.lib/pdes/poisson.cc @@ -109,12 +109,14 @@ main() mito::functions::sin(std::numbers::pi * x) * mito::functions::sin(std::numbers::pi * y); // compute the L2 error - auto error_L2 = discrete_system.compute_l2_error(u_ex); + auto error_L2 = + discrete_system.compute_l2_error(mito::fem::domain_field(u_ex)); // report channel << "L2 error: " << error_L2 << journal::endl; // compute the H1 error - auto error_H1 = discrete_system.compute_h1_error(u_ex); + auto error_H1 = + discrete_system.compute_h1_error(mito::fem::domain_field(u_ex)); // report channel << "H1 error: " << error_H1 << journal::endl; diff --git a/lib/mito/fem/DiscreteSystem.h b/lib/mito/fem/DiscreteSystem.h index 1c688788a..c3bd22f0b 100644 --- a/lib/mito/fem/DiscreteSystem.h +++ b/lib/mito/fem/DiscreteSystem.h @@ -36,8 +36,8 @@ namespace mito::fem { // solutions // the solution field type using solution_field_type = tensor::scalar_t; - // the nodal field type - using nodal_field_type = discrete::nodal_field_t; + // the fem field type + using fem_field_type = fem_field_t; // the number of nodes per element static constexpr int n_element_nodes = element_type::n_nodes; @@ -49,7 +49,8 @@ namespace mito::fem { _function_space(function_space), _weakform(weakform), _equation_map(), - _solution_field(nodal_field(function_space, label + ".solution")), + _solution_field( + function_space.template fem_field(label + ".solution")), _linear_system(label) { // make a channel @@ -198,13 +199,11 @@ namespace mito::fem { // get the node map from the function space auto node_map = _function_space.node_map(); - // fill information in nodal field - for (auto & [node, value] : _solution_field) { - // get the equation number of {node} - int eq = _equation_map.at(node); + // fill information in finite element field + for (auto & [node, eq] : _equation_map) { if (eq != -1) { - // read the solution at {eq} - value = u[eq]; + // note the solution on the solution field + _solution_field(node) = u[eq]; } } @@ -212,8 +211,8 @@ namespace mito::fem { return; } - // accessor to the solution nodal field - constexpr auto solution() const noexcept -> const nodal_field_type & + // accessor to the solution finite element field + constexpr auto solution() const noexcept -> const fem_field_type & { return _solution_field; } @@ -223,30 +222,19 @@ namespace mito::fem { // compute the L2 norm of the solution template - constexpr auto compute_l2_error(const fields::field_c auto & u_exact) const - -> tensor::scalar_t + constexpr auto compute_l2_error(const auto /*TOFIX*/ & u_exact) const -> tensor::scalar_t { // initialize the norm auto norm = tensor::scalar_t{ 0.0 }; - // helper lambda to assemble the solution on a given element - constexpr auto _assemble_element_solution = []( - const auto & element, - const auto & _solution_field, - tensor::integer_sequence) { - // assemble the solution field from the shape functions - return ( - (element.template shape() * _solution_field(element.connectivity()[a])) - + ...); - }; - // loop on all the cells of the mesh for (const auto & element : _function_space.elements()) { // assemble the solution field on this element - auto u_numerical = _assemble_element_solution( - element, _solution_field, tensor::make_integer_sequence()); + auto u_numerical = _solution_field.localize(element); + // localize the exact solution on this element + auto u_exact_local = u_exact.localize(element); // assemble the elementary error field as a function of the barycentric coordinates - auto u_error = u_numerical - u_exact(element.parametrization()); + auto u_error = u_numerical - u_exact_local; // compute the elementary contribution to the L2 norm norm += @@ -258,52 +246,25 @@ namespace mito::fem { // compute the H1 norm of the solution template - constexpr auto compute_h1_error(const fields::field_c auto & u_exact) const - -> tensor::scalar_t + constexpr auto compute_h1_error(const auto /*TOFIX*/ & u_exact) const -> tensor::scalar_t { // initialize the norm auto norm = tensor::scalar_t{ 0.0 }; - // QUESTION: is there a proper place to put this method? Perhaps this belongs to the - // isoparametric element class? - // helper lambda to assemble the solution on a given element - constexpr auto _assemble_element_solution = []( - const auto & element, - const auto & _solution_field, - tensor::integer_sequence) { - // assemble the solution field from the shape functions - return ( - (element.template shape() * _solution_field(element.connectivity()[a])) - + ...); - }; - - // TOFIX: we should deduce this from {_assemble_element_solution} - // helper lambda to assemble the solution on a given element - constexpr auto _assemble_element_solution_gradient = - []( - const auto & element, const auto & _solution_field, - tensor::integer_sequence) { - // assemble the solution field from the shape functions - return ( - (element.template gradient() - * _solution_field(element.connectivity()[a])) - + ...); - }; - // loop on all the cells of the mesh for (const auto & element : _function_space.elements()) { // assemble the solution field on this element - auto u_numerical = _assemble_element_solution( - element, _solution_field, tensor::make_integer_sequence()); + auto u_numerical = _solution_field.localize(element); + // localize the exact solution on this element + auto u_exact_local = u_exact.localize(element); + // assemble the elementary error field as a function of the barycentric coordinates + auto u_error = u_numerical - u_exact_local; + // assemble the gradient of the solution field on this element + auto u_numerical_gradient = fields::gradient(u_numerical); + // localize the gradient of the exact solution on this element + auto u_exact_gradient = fields::gradient(u_exact_local); // assemble the elementary error field as a function of the barycentric coordinates - auto u_error = u_numerical - u_exact(element.parametrization()); - // - auto u_numerical_gradient = _assemble_element_solution_gradient( - element, _solution_field, tensor::make_integer_sequence()); - // assemble the elementary error gradient as a function of the barycentric - // coordinates - auto u_error_gradient = - u_numerical_gradient - fields::gradient(u_exact)(element.parametrization()); + auto u_error_gradient = u_numerical_gradient - u_exact_gradient; // compute the elementary contributions to the H1 norm norm += blocks::l2_norm_block(u_error).compute(element) @@ -324,8 +285,8 @@ namespace mito::fem { // the equation map equation_map_type _equation_map; - // the solution nodal field - nodal_field_type _solution_field; + // the solution finite element field + fem_field_type _solution_field; // the linear system of equations linear_system_type _linear_system; From 69dd569738060a798015ec024628e17d301a7ef0 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 14:27:45 +0100 Subject: [PATCH 29/34] fem: remove unused {nodal_field} factory Nodal fields are now handed out by the function space. --- lib/mito/fem/api.h | 4 ---- lib/mito/fem/factories.h | 15 --------------- 2 files changed, 19 deletions(-) diff --git a/lib/mito/fem/api.h b/lib/mito/fem/api.h index 1ae53c9c3..2a4f7215e 100644 --- a/lib/mito/fem/api.h +++ b/lib/mito/fem/api.h @@ -9,10 +9,6 @@ namespace mito::fem { - // nodal field factory - template - constexpr auto nodal_field(const functionSpaceT & function_space, std::string name); - // domain field alias template using domain_field_t = DomainField; diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index 45ee0ac25..9ac305241 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -9,21 +9,6 @@ namespace mito::fem { - // nodal field factory - template - constexpr auto nodal_field(const functionSpaceT & function_space, std::string name) - { - // assemble the node type - using node_type = functionSpaceT::discretization_node_type; - - // get the nodes in the mesh - std::unordered_set> nodes; - get_discretization_nodes(function_space, nodes); - - // build a nodal field on the discretization nodes collected from the function space - return discrete::nodal_field_t(nodes, name); - } - // domain field factory template constexpr auto domain_field(const F & field) From 1862f9ca725434e09db9e0c23d62d1a1cbf3019a Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 14:28:03 +0100 Subject: [PATCH 30/34] fem: add concept for localizable fields --- lib/mito/fem/forward.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/mito/fem/forward.h b/lib/mito/fem/forward.h index f5c9b3bdb..a5f50008c 100644 --- a/lib/mito/fem/forward.h +++ b/lib/mito/fem/forward.h @@ -37,6 +37,12 @@ namespace mito::fem { // class finite element field template class FemField; + + // concept of a localizable field + template + concept localizable_field_c = requires(const F & f, const E & e) { + { f.localize(e) }; + }; } From 19ef51bdbbbc2a670e1c9e4db6c254f36d6a94fa Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 14:29:15 +0100 Subject: [PATCH 31/34] fem: calculation of the H1 and L2 norms now happen through free functions These functions should not be tied to the {DiscreteSystem}. --- benchmarks/mito.lib/pdes/poisson.cc | 11 ++-- lib/mito/fem/DiscreteSystem.h | 55 -------------------- lib/mito/fem/norms.h | 80 +++++++++++++++++++++++++++++ lib/mito/fem/public.h | 3 ++ 4 files changed, 90 insertions(+), 59 deletions(-) create mode 100644 lib/mito/fem/norms.h diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc index d716d34c1..847c6e6ad 100644 --- a/benchmarks/mito.lib/pdes/poisson.cc +++ b/benchmarks/mito.lib/pdes/poisson.cc @@ -104,19 +104,22 @@ main() // free the solver solver.destroy(); + // get the solution field + auto & solution = discrete_system.solution(); + // the exact solution field auto u_ex = mito::functions::sin(std::numbers::pi * x) * mito::functions::sin(std::numbers::pi * y); // compute the L2 error - auto error_L2 = - discrete_system.compute_l2_error(mito::fem::domain_field(u_ex)); + auto error_L2 = mito::fem::compute_l2_norm( + function_space, solution, mito::fem::domain_field(u_ex)); // report channel << "L2 error: " << error_L2 << journal::endl; // compute the H1 error - auto error_H1 = - discrete_system.compute_h1_error(mito::fem::domain_field(u_ex)); + auto error_H1 = mito::fem::compute_h1_norm( + function_space, solution, mito::fem::domain_field(u_ex)); // report channel << "H1 error: " << error_H1 << journal::endl; diff --git a/lib/mito/fem/DiscreteSystem.h b/lib/mito/fem/DiscreteSystem.h index c3bd22f0b..babd6605d 100644 --- a/lib/mito/fem/DiscreteSystem.h +++ b/lib/mito/fem/DiscreteSystem.h @@ -220,61 +220,6 @@ namespace mito::fem { // accessor to the number of equations constexpr auto n_equations() const noexcept -> int { return _n_equations; } - // compute the L2 norm of the solution - template - constexpr auto compute_l2_error(const auto /*TOFIX*/ & u_exact) const -> tensor::scalar_t - { - // initialize the norm - auto norm = tensor::scalar_t{ 0.0 }; - - // loop on all the cells of the mesh - for (const auto & element : _function_space.elements()) { - // assemble the solution field on this element - auto u_numerical = _solution_field.localize(element); - // localize the exact solution on this element - auto u_exact_local = u_exact.localize(element); - // assemble the elementary error field as a function of the barycentric coordinates - auto u_error = u_numerical - u_exact_local; - - // compute the elementary contribution to the L2 norm - norm += - blocks::l2_norm_block(u_error).compute(element); - } - - return std::sqrt(norm); - } - - // compute the H1 norm of the solution - template - constexpr auto compute_h1_error(const auto /*TOFIX*/ & u_exact) const -> tensor::scalar_t - { - // initialize the norm - auto norm = tensor::scalar_t{ 0.0 }; - - // loop on all the cells of the mesh - for (const auto & element : _function_space.elements()) { - // assemble the solution field on this element - auto u_numerical = _solution_field.localize(element); - // localize the exact solution on this element - auto u_exact_local = u_exact.localize(element); - // assemble the elementary error field as a function of the barycentric coordinates - auto u_error = u_numerical - u_exact_local; - // assemble the gradient of the solution field on this element - auto u_numerical_gradient = fields::gradient(u_numerical); - // localize the gradient of the exact solution on this element - auto u_exact_gradient = fields::gradient(u_exact_local); - // assemble the elementary error field as a function of the barycentric coordinates - auto u_error_gradient = u_numerical_gradient - u_exact_gradient; - // compute the elementary contributions to the H1 norm - norm += - blocks::l2_norm_block(u_error).compute(element) - + blocks::l2_norm_block(u_error_gradient) - .compute(element); - } - - return std::sqrt(norm); - } - private: // a const reference to the function space const function_space_type & _function_space; diff --git a/lib/mito/fem/norms.h b/lib/mito/fem/norms.h new file mode 100644 index 000000000..e826620ab --- /dev/null +++ b/lib/mito/fem/norms.h @@ -0,0 +1,80 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + // compute L2 norm on a given function space of the difference between two localizable fields + template + requires( + localizable_field_c + && localizable_field_c) + constexpr auto compute_l2_norm( + const functionSpaceT & function_space, const F1 & u1, const F2 & u2) -> tensor::scalar_t + { + // get the element type + using element_type = typename functionSpaceT::element_type; + + // initialize the norm + auto norm = tensor::scalar_t{ 0.0 }; + + // loop on all the elements of the function space + for (const auto & element : function_space.elements()) { + // localize {u1} on this element + auto u1_local = u1.localize(element); + // localize {u2} on this element + auto u2_local = u2.localize(element); + // compute the elementary contribution to the norm + norm += blocks::l2_norm_block(u1_local - u2_local) + .compute(element); + } + + // take the square root of the accumulated norm + return std::sqrt(norm); + } + + // compute H1 norm on a given function space of the difference between two localizable fields + template + requires( + localizable_field_c + && localizable_field_c) + constexpr auto compute_h1_norm( + const functionSpaceT & function_space, const F1 & u1, const F2 & u2) -> tensor::scalar_t + { + // get the element type + using element_type = typename functionSpaceT::element_type; + + // initialize the norm + auto norm = tensor::scalar_t{ 0.0 }; + + // loop on all the elements of the function space + for (const auto & element : function_space.elements()) { + // localize {u1} on this element + auto u1_local = u1.localize(element); + // localize {u2} on this element + auto u2_local = u2.localize(element); + // assemble the gradient of the solution field on this element + auto u1_local_gradient = fields::gradient(u1_local); + // localize the gradient of the exact solution on this element + auto u2_local_gradient = fields::gradient(u2_local); + // compute the elementary contributions to the H1 norm + norm += blocks::l2_norm_block(u1_local - u2_local) + .compute(element) + + blocks::l2_norm_block( + u1_local_gradient - u2_local_gradient) + .compute(element); + } + + // take the square root of the accumulated norm + return std::sqrt(norm); + } + +} + + +// end of file diff --git a/lib/mito/fem/public.h b/lib/mito/fem/public.h index 6ec9a959d..81c5c810f 100644 --- a/lib/mito/fem/public.h +++ b/lib/mito/fem/public.h @@ -35,4 +35,7 @@ // factories implementation #include "factories.h" +// norms implementation +#include "norms.h" + // end of file From b65ad232ac10ad89c220d9f6d62090cd33afc96d Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 16:12:17 +0100 Subject: [PATCH 32/34] fem: implement {FemField} iterators --- lib/mito/fem/FemField.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/mito/fem/FemField.h b/lib/mito/fem/FemField.h index d4d495b58..0bcd500b8 100644 --- a/lib/mito/fem/FemField.h +++ b/lib/mito/fem/FemField.h @@ -79,6 +79,12 @@ namespace mito::fem { element, _nodal_field, tensor::make_integer_sequence{}); } + // iterators on the nodal field + constexpr auto begin() noexcept { return _nodal_field.begin(); } + constexpr auto end() noexcept { return _nodal_field.end(); } + constexpr auto begin() const noexcept { return _nodal_field.begin(); } + constexpr auto end() const noexcept { return _nodal_field.end(); } + private: // the field with the nodal values of the field nodal_field_type _nodal_field; From 35085724053c45109129b28429e43585a8a22618 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 8 Dec 2025 16:13:01 +0100 Subject: [PATCH 33/34] tests/fem: add unit test for {FemField} --- .cmake/mito_tests_mito_lib.cmake | 1 + tests/mito.lib/fem/fem_field.cc | 84 ++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+) create mode 100644 tests/mito.lib/fem/fem_field.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 2f027acac..8ee95ba5b 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -52,6 +52,7 @@ mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p2.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc) mito_test_driver(tests/mito.lib/fem/domain_field.cc) +mito_test_driver(tests/mito.lib/fem/fem_field.cc) # io mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc) diff --git a/tests/mito.lib/fem/fem_field.cc b/tests/mito.lib/fem/fem_field.cc new file mode 100644 index 000000000..7afaa4954 --- /dev/null +++ b/tests/mito.lib/fem/fem_field.cc @@ -0,0 +1,84 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// cartesian coordinates in 2D +using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; + +// simplicial cells in 2D +using cell_t = mito::geometry::triangle_t<2>; +// first degree finite elements +constexpr int degree = 1; +// assemble the finite element type +using finite_element_t = mito::fem::isoparametric_simplex_t; +// the x scalar field in 2D +constexpr auto x = mito::functions::component; +// the y scalar field in 2D +constexpr auto y = mito::functions::component; + + +TEST(Fem, FemField) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system(); + + // read the mesh of a square in 2D + std::ifstream fileStream("square.summit"); + auto mesh = mito::io::summit::reader(fileStream, coord_system); + + // create the body manifold + auto manifold = mito::manifolds::manifold(mesh, coord_system); + + // TOFIX: it should not be mandatory to set constraints to create a function space, let's remove + // this bit once we implement constraints properly + // get the boundary mesh + auto boundary_mesh = mito::mesh::boundary(mesh); + + // the zero field + auto zero = mito::functions::zero; + + // set homogeneous Dirichlet boundary condition + auto constraints = mito::constraints::dirichlet_bc(boundary_mesh, zero); + + // the function space (linear elements on the manifold) + auto function_space = mito::fem::function_space(manifold, constraints); + + // get the map between mesh nodes and discretization nodes + const auto & node_map = function_space.node_map(); + + // get a scalar-valued fem field on the function space + auto fem_field = function_space.fem_field("linear field"); + + // create a continuous linear field + auto field = x + y; + + // loop on all the elements of the functions space + for (const auto & element : function_space.elements()) { + // loop on all the nodes of the element + for (const auto & node : element.cell().nodes()) { + // compute the coordinates of the node + auto coords = coord_system.coordinates(node->point()); + // set the field value at {node} + fem_field(node_map.at(node)) = field(coords); + } + } + + // loop on all the elements of the functions space + for (const auto & element : function_space.elements()) { + // localize the field on {element} + auto local = fem_field.localize(element); + // get the coordinates of the center of the element + auto center_coords = mito::geometry::barycenter(element.cell(), coord_system); + // check that the results is the same regardless of whether we evaluate the localized field + // or the original continuous field + EXPECT_DOUBLE_EQ(local({ 1.0 / 3.0, 1.0 / 3.0 }), field(center_coords)); + } + + // all done + return; +} \ No newline at end of file From c8e9e707b458c0d2cf7762162b76283c4b0b850c Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Tue, 16 Dec 2025 15:07:11 +0100 Subject: [PATCH 34/34] benchmarks: fix minor compiler error in {poisson} benchmark --- benchmarks/mito.lib/pdes/poisson.cc | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc index 847c6e6ad..08d6e8bcc 100644 --- a/benchmarks/mito.lib/pdes/poisson.cc +++ b/benchmarks/mito.lib/pdes/poisson.cc @@ -105,7 +105,7 @@ main() solver.destroy(); // get the solution field - auto & solution = discrete_system.solution(); + const auto & solution = discrete_system.solution(); // the exact solution field auto u_ex = @@ -137,8 +137,6 @@ main() // write output file writer.write(); - // get the solution field - const auto & solution = discrete_system.solution(); // write mesh to vtk file auto writer_solution = mito::io::vtk::field_writer("poisson_square_solution", function_space, coord_system);