diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake
index 88ab480e6..8ee95ba5b 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
@@ -51,6 +51,8 @@ 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)
+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/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<coordinates_t, 1>;
 
     // 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<mito::quadrature::GAUSS, 2>(manifold);
diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc
index f7e30148e..08d6e8bcc 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<coordinates_t>);
+    auto zero = mito::functions::zero<coordinates_t>;
 
     // 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<finite_element_t, quadrature_rule_t>();
 
     // 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
@@ -105,17 +104,22 @@ main()
     // free the solver
     solver.destroy();
 
+    // get the solution field
+    const auto & solution = discrete_system.solution();
+
     // 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<quadrature_rule_t>(u_ex);
+    auto error_L2 = mito::fem::compute_l2_norm<quadrature_rule_t>(
+        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<quadrature_rule_t>(u_ex);
+    auto error_H1 = mito::fem::compute_h1_norm<quadrature_rule_t>(
+        function_space, solution, mito::fem::domain_field(u_ex));
     // report
     channel << "H1 error: " << error_H1 << journal::endl;
 
@@ -133,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);
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<mito::tensor::scalar_t(const coordinates_2D_t &)>;
-    using scalar_function_functor_2D_t = mito::functions::FunctionFromFunctor<scalar_function_2D_t>;
-    using scalar_field_2D_t =
-        decltype(mito::fields::field(std::declval<scalar_function_functor_2D_t>()));
+    using scalar_field_2D_t = mito::functions::FunctionFromFunctor<scalar_function_2D_t>;
     mito::py::class_<scalar_field_2D_t>(m, "ScalarField2D")
         // the constructor
         .def(mito::py::init<scalar_function_2D_t>())
@@ -90,9 +88,7 @@ PYBIND11_MODULE(mito, m)
 
     // the mito scalar field 3D
     using scalar_function_3D_t = std::function<mito::tensor::scalar_t(const coordinates_3D_t &)>;
-    using scalar_function_functor_3D_t = mito::functions::FunctionFromFunctor<scalar_function_3D_t>;
-    using scalar_field_3D_t =
-        decltype(mito::fields::field(std::declval<scalar_function_functor_3D_t>()));
+    using scalar_field_3D_t = mito::functions::FunctionFromFunctor<scalar_function_3D_t>;
     mito::py::class_<scalar_field_3D_t>(m, "ScalarField3D")
         // the constructor
         .def(mito::py::init<scalar_function_3D_t>())
diff --git a/lib/mito/discrete/DiscreteField.h b/lib/mito/discrete/DiscreteField.h
index 5ebebb933..e9d653609 100644
--- a/lib/mito/discrete/DiscreteField.h
+++ b/lib/mito/discrete/DiscreteField.h
@@ -12,15 +12,16 @@ namespace mito::discrete {
     template <class keyT, class valueT>
     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<key_type, value_type, utilities::hash_function<key_type>>;
+            std::unordered_map<input_type, value_type, utilities::hash_function<input_type>>;
 
       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,10 +72,14 @@ 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)
+            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;
@@ -99,7 +104,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
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 <string>
+#include <cassert>
 
 // support
 #include "../mesh.h"
diff --git a/lib/mito/fem/DiscreteSystem.h b/lib/mito/fem/DiscreteSystem.h
index de88210da..babd6605d 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<solution_field_type>;
+        // the fem field type
+        using fem_field_type = fem_field_t<solution_field_type>;
         // 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<solution_field_type>(function_space, label + ".solution")),
+            _solution_field(
+                function_space.template fem_field<solution_field_type>(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;
         }
@@ -221,103 +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 <class quadratureRuleT>
-        constexpr auto compute_l2_error(const fields::field_c auto & 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 = []<int... a>(
-                                                            const auto & element,
-                                                            const auto & _solution_field,
-                                                            tensor::integer_sequence<a...>) {
-                // assemble the solution field from the shape functions
-                return (
-                    (element.template shape<a>() * _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<n_element_nodes>());
-                // assemble the elementary error field as a function of the barycentric coordinates
-                auto u_error =
-                    u_numerical - u_exact.function()(element.parametrization().function());
-
-                // compute the elementary contribution to the L2 norm
-                norm += blocks::l2_norm_block<element_type, quadratureRuleT>(u_error.function())
-                            .compute(element);
-            }
-
-            return std::sqrt(norm);
-        }
-
-        // compute the H1 norm of the solution
-        template <class quadratureRuleT>
-        constexpr auto compute_h1_error(const fields::field_c auto & 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 = []<int... a>(
-                                                            const auto & element,
-                                                            const auto & _solution_field,
-                                                            tensor::integer_sequence<a...>) {
-                // assemble the solution field from the shape functions
-                return (
-                    (element.template shape<a>() * _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 =
-                []<int... a>(
-                    const auto & element, const auto & _solution_field,
-                    tensor::integer_sequence<a...>) {
-                    // assemble the solution field from the shape functions
-                    return (
-                        (element.template gradient<a>()
-                         * _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<n_element_nodes>());
-                // 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_numerical_gradient = _assemble_element_solution_gradient(
-                    element, _solution_field, tensor::make_integer_sequence<n_element_nodes>());
-                // 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());
-                // compute the elementary contributions to the H1 norm
-                norm += blocks::l2_norm_block<element_type, quadratureRuleT>(u_error.function())
-                            .compute(element)
-                      + blocks::l2_norm_block<element_type, quadratureRuleT>(
-                            u_error_gradient.function())
-                            .compute(element);
-            }
-
-            return std::sqrt(norm);
-        }
-
       private:
         // a const reference to the function space
         const function_space_type & _function_space;
@@ -328,8 +230,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;
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 <fields::field_c F>
+    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 <class elementT>
+        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/FemField.h b/lib/mito/fem/FemField.h
new file mode 100644
index 000000000..0bcd500b8
--- /dev/null
+++ b/lib/mito/fem/FemField.h
@@ -0,0 +1,96 @@
+// -*- 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 fieldValueT>
+    class FemField {
+
+      private:
+        // the field value type
+        using field_value_type = fieldValueT;
+        // the nodal field type
+        using nodal_field_type = discrete::nodal_field_t<field_value_type>;
+        // 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(std::move(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 operator()(const node_type & node) const -> const field_value_type &
+        {
+            return _nodal_field(node);
+        }
+
+        // get a mutable nodal value at {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?
+        // TODO: add concept that constrains {elementT} to be a finite element type
+        // localize the function on {element}
+        template <class elementT>
+        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 = []<int... a>(
+                                           const element_type & element,
+                                           const nodal_field_type & field,
+                                           tensor::integer_sequence<a...>) {
+                // assemble the solution field from the shape functions
+                return ((element.template shape<a>() * field(element.connectivity()[a])) + ...);
+            };
+            return _assemble(
+                element, _nodal_field, tensor::make_integer_sequence<element_type::n_nodes>{});
+        }
+
+        // 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;
+    };
+
+}    // namespace mito
+
+
+// end of file
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<mesh_node_type>>;
-
+        // a finite element field type
+        template <class fieldValueT>
+        using fem_field_type = fem_field_t<fieldValueT>;
 
       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 <class fieldValueT>
+        constexpr auto fem_field(std::string name) const -> fem_field_type<fieldValueT>
+        {
+            // assemble the node type
+            using node_type = discretization_node_type;
+
+            // get the discretization nodes
+            std::unordered_set<node_type, utilities::hash_function<node_type>> nodes;
+            get_discretization_nodes(*this, nodes);
+
+            // build a nodal field on the discretization nodes collected from the function space
+            return fem_field_t<fieldValueT>(discrete::nodal_field_t<fieldValueT>(nodes, name));
+        }
+
       private:
         // a collection of finite elements
         elements_type _elements;
diff --git a/lib/mito/fem/api.h b/lib/mito/fem/api.h
index b0222bb33..2a4f7215e 100644
--- a/lib/mito/fem/api.h
+++ b/lib/mito/fem/api.h
@@ -9,9 +9,17 @@
 
 namespace mito::fem {
 
-    // nodal field factory
-    template <class Y, function_space_c functionSpaceT>
-    constexpr auto nodal_field(const functionSpaceT & function_space, std::string name);
+    // domain field alias
+    template <fields::field_c F>
+    using domain_field_t = DomainField<F>;
+
+    // domain field factory
+    template <fields::field_c F>
+    constexpr auto domain_field(const F &);
+
+    // finite element field alias
+    template <class fieldValueT>
+    using fem_field_t = FemField<fieldValueT>;
 
     // the possible discretization types: continuous Galerking (CG) vs. discontinuous Galerkin (DG)
     enum class discretization_t { CG, DG };
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/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/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
index fcd527d92..ab58d5f29 100644
--- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
+++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
@@ -64,19 +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 mito::fields::linear_combination(
-                std::array{ _x0, _x1, _x2 }, phi_0, phi_1, phi_2);
-        }
-
         // get the shape function associated with local node {a}
         template <int a>
         requires(a >= 0 && a < n_nodes)
@@ -113,8 +100,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..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 mito::fields::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 &
         {
@@ -128,8 +108,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/factories.h b/lib/mito/fem/factories.h
index 6148588a8..9ac305241 100644
--- a/lib/mito/fem/factories.h
+++ b/lib/mito/fem/factories.h
@@ -9,19 +9,12 @@
 
 namespace mito::fem {
 
-    // nodal field factory
-    template <class Y, function_space_c functionSpaceT>
-    constexpr auto nodal_field(const functionSpaceT & function_space, std::string name)
+    // domain field factory
+    template <fields::field_c F>
+    constexpr auto domain_field(const F & field)
     {
-        // assemble the node type
-        using node_type = functionSpaceT::discretization_node_type;
-
-        // get the nodes in the mesh
-        std::unordered_set<node_type, utilities::hash_function<node_type>> 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<Y>(nodes, name);
+        // build a domain field from the given function
+        return domain_field_t<F>(field);
     }
 
     // TOFIX: create a constructor that takes no constraints
diff --git a/lib/mito/fem/forward.h b/lib/mito/fem/forward.h
index d212e8f4a..a5f50008c 100644
--- a/lib/mito/fem/forward.h
+++ b/lib/mito/fem/forward.h
@@ -29,6 +29,20 @@ namespace mito::fem {
     // class discrete system
     template <function_space_c functionSpaceT, class linearSystemT>
     class DiscreteSystem;
+
+    // class domain field
+    template <fields::field_c F>
+    class DomainField;
+
+    // class finite element field
+    template <class fieldValueT>
+    class FemField;
+
+    // concept of a localizable field
+    template <class F, class E>
+    concept localizable_field_c = requires(const F & f, const E & e) {
+        { f.localize(e) };
+    };
 }
 
 
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 <class quadratureRuleT, class functionSpaceT, class F1, class F2>
+    requires(
+        localizable_field_c<F1, typename functionSpaceT::element_type>
+        && localizable_field_c<F2, typename functionSpaceT::element_type>)
+    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<element_type, quadratureRuleT>(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 <class quadratureRuleT, class functionSpaceT, class F1, class F2>
+    requires(
+        localizable_field_c<F1, typename functionSpaceT::element_type>
+        && localizable_field_c<F2, typename functionSpaceT::element_type>)
+    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<element_type, quadratureRuleT>(u1_local - u2_local)
+                        .compute(element)
+                  + blocks::l2_norm_block<element_type, quadratureRuleT>(
+                        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 1ce8000ab..81c5c810f 100644
--- a/lib/mito/fem/public.h
+++ b/lib/mito/fem/public.h
@@ -19,8 +19,13 @@
 // 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"
 
@@ -30,4 +35,7 @@
 // factories implementation
 #include "factories.h"
 
+// norms implementation
+#include "norms.h"
+
 // end of file
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
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 <functions::function_c F>
-    requires(geometry::coordinates_c<typename F::input_type>)
-    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..3c0234bcb 100644
--- a/lib/mito/fields/api.h
+++ b/lib/mito/fields/api.h
@@ -9,29 +9,6 @@
 
 namespace mito::fields {
 
-    // field alias
-    template <functions::function_c F>
-    using field_t = Field<F>;
-
-    // factory for fields
-    template <functions::function_c F>
-    constexpr auto field(const F & f) -> field_t<F>;
-
-    // factory for fields
-    template <functions::function_c F>
-    constexpr auto field(F && f) -> field_t<F>;
-
-    // factory for fields (from functors)
-    template <functions::functor_c F>
-    constexpr auto field(F && f) -> field_t<F>;
-
-    // uniform field
-    template <geometry::coordinates_c coordType, class Y>
-    constexpr auto uniform_field(const Y & constant);
-
-    // the order N identity tensor in D dimensions
-    template <geometry::coordinates_c coordType, int N>
-    constexpr auto identity_tensor_field = uniform_field<coordType>(tensor::identity<N>);
 }
 
 
diff --git a/lib/mito/fields/differential.h b/lib/mito/fields/differential.h
index e18646475..9720ce8fc 100644
--- a/lib/mito/fields/differential.h
+++ b/lib/mito/fields/differential.h
@@ -10,27 +10,22 @@
 // Differential operators on Fields
 namespace mito::fields {
 
-    // the {I...}-th first partial derivative of a field
-    template <int... I>
-    constexpr auto derivative(const field_c auto & f)
-    {
-        // the {I...}-th first partial derivative
-        return field(mito::functions::derivative<I...>(f.function()));
-    }
-
     // function to compute the gradient of a scalar field
     template <scalar_field_c F>
     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;
 
         // helper function to compute the gradient of a scalar
         constexpr auto _grad = []<size_t... I>(const F & field, std::index_sequence<I...>) {
             // the vector of the partial derivatives
-            return ((derivative<I>(field) * uniform_field<coordinate_t>(tensor::e<I, D>)) + ...);
+            return (
+                (functions::derivative<I>(field)
+                 * functions::constant<coordinate_t>(tensor::e<I, D>))
+                + ...);
         };
 
         return _grad(field, std::make_index_sequence<D>{});
@@ -41,7 +36,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
@@ -54,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 (
-                    (derivative<K>(field * uniform_field<coordinate_t>(tensor::e<I, N>))
-                     * uniform_field<coordinate_t>(tensor::unit<tensor::matrix_t<N, D>, I, K>))
+                    (functions::derivative<K>(
+                         field * functions::constant<coordinate_t>(tensor::e<I, N>))
+                     * functions::constant<coordinate_t>(
+                         tensor::unit<tensor::matrix_t<N, D>, I, K>))
                     + ...);
             };
 
@@ -72,14 +69,17 @@ 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;
 
         // helper function to compute the divergence of a vector field
         constexpr auto _div = []<size_t... I>(const F & field, std::index_sequence<I...>) {
             // the summation of the I-th partial derivative of the I-th component
-            return ((derivative<I>(field * uniform_field<coordinate_t>(tensor::e<I, D>))) + ...);
+            return (
+                (functions::derivative<I>(
+                    field * functions::constant<coordinate_t>(tensor::e<I, D>)))
+                + ...);
         };
 
         // all done
@@ -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;
 
@@ -102,10 +102,10 @@ namespace mito::fields {
                 // the vector of the partial derivatives
                 // (\partial field_KI / \partial x_I) * e_K
                 return (
-                    (derivative<I>(
-                         (field * uniform_field<coordinate_t>(tensor::e<I, D>))
-                         * uniform_field<coordinate_t>(tensor::e<K, D>))
-                     * uniform_field<coordinate_t>(tensor::e<K, D>))
+                    (functions::derivative<I>(
+                         (field * functions::constant<coordinate_t>(tensor::e<I, D>))
+                         * functions::constant<coordinate_t>(tensor::e<K, D>))
+                     * functions::constant<coordinate_t>(tensor::e<K, D>))
                     + ...);
             };
 
diff --git a/lib/mito/fields/factories.h b/lib/mito/fields/factories.h
index 4360a2dc9..3c0234bcb 100644
--- a/lib/mito/fields/factories.h
+++ b/lib/mito/fields/factories.h
@@ -9,33 +9,6 @@
 
 namespace mito::fields {
 
-    // factory for fields
-    template <functions::function_c F>
-    constexpr auto field(const F & f) -> field_t<F>
-    {
-        return field_t<F>(f);
-    }
-
-    // factory for fields
-    template <functions::function_c F>
-    constexpr auto field(F && f) -> field_t<F>
-    {
-        return field_t<F>(std::forward<F>(f));
-    }
-
-    // factory for fields (from functors)
-    template <functions::functor_c F>
-    constexpr auto field(F && f)
-    {
-        return field(mito::functions::function(std::forward<F>(f)));
-    }
-
-    // uniform field
-    template <geometry::coordinates_c coordType, class Y>
-    constexpr auto uniform_field(const Y & constant)
-    {
-        return field(mito::functions::constant<coordType>(constant));
-    }
 }
 
 
diff --git a/lib/mito/fields/fields_algebra.h b/lib/mito/fields/fields_algebra.h
index b9c1dece7..d9734c32f 100644
--- a/lib/mito/fields/fields_algebra.h
+++ b/lib/mito/fields/fields_algebra.h
@@ -10,49 +10,14 @@
 // Algebraic operations on Fields
 namespace mito::fields {
 
-    // f^(-1)
-    template <tensor_field_c F>
-    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 <tensor_field_c F>
-    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 {
-            return determinant(f(x));
-        }));
-    }
-
-    // f transpose
-    template <tensor_field_c F>
-    constexpr auto transpose(const F & f)
-    {
-        return field(functions::transpose(f.function()));
-    }
-
-    // sqrt(f)
-    template <scalar_field_c F>
-    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
     template <one_form_field_c F1, one_form_field_c F2>
     constexpr auto tensor(const F1 & fA, const F2 & fB)
     requires(compatible_fields_c<F1, F2>)
     {
-        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 +25,7 @@ namespace mito::fields {
     constexpr auto tensor(const F1 & fA, const F2 & fB, const F3 & fC)
     requires(compatible_fields_c<F1, F2> && compatible_fields_c<F1, F3>)
     {
-        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 +43,9 @@ namespace mito::fields {
     constexpr auto wedge(const F1 & fA, const F2 & fB)
     requires(compatible_fields_c<F1, F2>)
     {
-        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 +53,10 @@ namespace mito::fields {
     constexpr auto wedge(const F1 & fA, const F2 & fB, const F3 & fC)
     requires(compatible_fields_c<F1, F2> && compatible_fields_c<F1, F3>)
     {
-        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 <typename T, field_c... Fs>
-    constexpr auto linear_combination(std::array<T, sizeof...(Fs)> coeffs, const Fs &... fs)
-    requires(compatible_fields_c<Fs...>)
-    {
-        return field(functions::linear_combination(coeffs, fs.function()...));
-    }
-
-    // addition of fields fa + fb
-    template <field_c F1, field_c F2>
-    constexpr auto operator+(const F1 & fA, const F2 & fB)
-    requires(compatible_fields_c<F1, F2>)
-    {
-        return field(fA.function() + fB.function());
-    }
-
-    // scalar + field
-    template <field_c F>
-    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 <field_c F>
-    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 <field_c F1, field_c F2>
-    constexpr auto operator*(const F1 & fA, const F2 & fB)
-    requires(compatible_fields_c<F1, F2>)
-    {
-        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 <scalar_field_c F>
-    constexpr auto operator/(const tensor::tensor_or_scalar_c auto & a, const F & f)
-    {
-        return field(a / f.function());
-    }
-
-    // field1 / field2
-    template <field_c F1, scalar_field_c F2>
-    constexpr auto operator/(const F1 & f1, const F2 & f2)
-    requires(compatible_fields_c<F1, F2>)
-    {
-        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 <field_c F1, functions::function_c F2>
-    constexpr auto operator+(const F1 & fA, const F2 & fB)
-    {
-        return field(fA.function() + fB);
-    }
-
-    // addition of fields fa + fb
-    template <functions::function_c F1, field_c F2>
-    constexpr auto operator+(const F1 & fA, const F2 & fB)
-    {
-        return fB + fA;
-    }
-
-    // product of fields
-    template <field_c F1, functions::function_c F2>
-    constexpr auto operator*(const F1 & fA, const F2 & fB)
-    {
-        return field(fA.function() * fB);
-    }
-
-    // product of fields
-    template <functions::function_c F1, field_c F2>
-    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 <field_c F1, functions::scalar_valued_function_c F2>
-    constexpr auto operator/(const F1 & f1, const F2 & f2)
-    {
-        return field(f1.function() / f2);
-    }
-
-    // field1 / field2
-    template <functions::function_c F1, scalar_field_c F2>
-    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 <functions::function_c F>
-    requires(geometry::coordinates_c<typename F::input_type>)
-    class Field;
-
-    // concept of a field
-    template <class FIELD>
-    concept field_c = requires(FIELD c) {
-        // require that F only binds to {Field<F>} specializations or derived classes
-        []<functions::function_c F>(const Field<F> &) {
-        }(c);
-    };
+    // concept of a field: a field is a function of coordinates
+    template <class F>
+    concept field_c = functions::function_c<F> && geometry::coordinates_c<typename F::input_type>;
 
     // concept of a field {FIELD} being a scalar field
     template <class FIELD>
@@ -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<FIELDS...>>::coordinates_type,
-             typename FIELDS::coordinates_type>
+             typename std::tuple_element_t<0, std::tuple<FIELDS...>>::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/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/api.h b/lib/mito/functions/api.h
index bb0162b8d..024545570 100644
--- a/lib/mito/functions/api.h
+++ b/lib/mito/functions/api.h
@@ -64,13 +64,9 @@ namespace mito::functions {
     constexpr auto function(F && f);
 
     // the function extracting the N-th component from an input argument of type {T}
-    template <class T, int N>
+    template <subscriptable_c T, int N>
     [[maybe_unused]] constexpr auto component = Component<T, N>();
 
-    // the function associating to a D-dimensional vector its N-th component
-    template <int N, int D>
-    [[maybe_unused]] constexpr auto x = component<tensor::vector_t<D>, N>;
-
     // the linear combination
     template <typename T, function_c... Funcs>
     constexpr auto linear_combination(std::array<T, sizeof...(Funcs)> coeffs, Funcs... funcs);
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/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<X, Y>(c);
     }
 
+    // returns the order N identity tensor as a function of X
+    template <class X, int N>
+    constexpr auto identity()
+    {
+        // the order N identity tensor
+        constexpr auto id = tensor::identity<N>;
+        // return the identity ad a constant function
+        return constant<X>(id);
+    }
+
     // returns the function wrapping the functor {f}
     template <functor_c F>
     constexpr auto function(F && f)
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 <class F>
     concept vector_function_c = vector_domain_function_c<F> and vector_valued_function_c<F>;
 
-    // concept of a function defined on tensors
+    // concept of a function defined on tensors or coordinates
     template <class F>
-    concept tensor_domain_function_c = function_c<F> and tensor::tensor_c<typename F::input_type>;
+    concept tensor_domain_function_c = function_c<F>
+                                   and (tensor::tensor_c<typename F::input_type>
+                                        or geometry::coordinates_c<typename F::input_type>);
 
     // concept of a tensor-valued function
     template <class 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 <function_c F>
+    requires(tensor::square_matrix_c<typename F::output_type>)
+    class Determinant :
+        public Function<typename F::input_type, typename F::output_type::scalar_type> {
+
+      public:
+        // the type of the function (what I derive from)
+        using function_type =
+            Function<typename F::input_type, typename F::output_type::scalar_type>;
+        // 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 <function_c H>
+        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 <function_c F>
     requires(tensor::square_matrix_c<typename F::output_type>)
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 <class T, int I>
+    template <subscriptable_c T, int I>
     class Component : public Function<T, tensor::scalar_t> {
 
       public:
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
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 <coordinate_system_c coordinateSystemT>
-        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<n_vertices - 1>([&]<int i>() {
-                // 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 = []<int... a>(
+                                           const auto & nodes, const auto & coordinate_system,
+                                           tensor::integer_sequence<a...>) {
+                // get the origin of the coordinate system
+                constexpr auto origin = coordinate_system.origin();
+
+                // assemble the parametrization as x0 * xi<0> + ...
+                // where {xi<a>} are the barycentric coordinates on the reference simplex and the
+                // {xa} are the position vectors of the nodes
+                return (
+                    ((reference_simplex_type::template xi<a>
+                      * (coordinate_system.coordinates(nodes[a]->point()) - origin)))
+                    + ...);
+            };
+            return _assemble(
+                _nodes, coordinate_system, tensor::make_integer_sequence<n_vertices>{});
         }
 
       private:
diff --git a/lib/mito/geometry/ReferenceSimplex.h b/lib/mito/geometry/ReferenceSimplex.h
index 6b94a0bd4..253f7d092 100644
--- a/lib/mito/geometry/ReferenceSimplex.h
+++ b/lib/mito/geometry/ReferenceSimplex.h
@@ -18,14 +18,29 @@ namespace mito::geometry {
         // the area of the reference simplex
         constexpr static double area = 1.0 / mito::tensor::factorial<order>();
 
+      private:
+        // helper to compute 1 - xi0 - xi1 - ... - xi(N-1)
+        template <int... a>
+        static constexpr auto _one_minus_xis(tensor::integer_sequence<a...>)
+        {
+            return (1.0 - ... - xi<a>);
+        }
+
       public:
         // the type of coordinates in the parametric space
         using parametric_coordinates_type = coordinates_t<N, CARTESIAN>;
 
         // the function extracting the I component of a parametric point
         template <int I>
-        static constexpr auto xi =
-            fields::field(functions::component<parametric_coordinates_type, I>);
+        static constexpr auto xi = [] {
+            if constexpr (I < N) {
+                // I-th parametric coordinate
+                return functions::component<parametric_coordinates_type, I>;
+            } else {
+                // xi<N> = 1 - xi<0> - xi<1> - ... - xi<N-1>
+                return _one_minus_xis(tensor::make_integer_sequence<N>{});
+            }
+        }();
     };
 
 }    // namespace mito
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<cartesian_coordinates_t<D>>(tensor::e<I, D>);
+            return functions::constant<cartesian_coordinates_t<D>>(tensor::e<I, D>);
         }
     };
 }
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<cartesian_coordinates_t<D>, D>;
+            return functions::identity<cartesian_coordinates_t<D>, D>();
         }
     };
 }
diff --git a/lib/mito/geometry/metric_space.h b/lib/mito/geometry/metric_space.h
index 75ef7c4f7..b7ef420ce 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<coordinates_type>::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 <int I>
@@ -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<D>{});
+            functions::sqrt(functions::determinant(g)) * _wedge(tensor::make_integer_sequence<D>{});
 
         // 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]<int... K>(tensor::integer_sequence<K...>) {
@@ -80,7 +80,7 @@ namespace mito::geometry {
 
             // all done
             return _one_form_components(tensor::make_integer_sequence<D>{});
-        }));
+        });
     }
 
     template <coordinates_c coordsT>
@@ -98,10 +98,10 @@ namespace mito::geometry {
     requires(fields::compatible_fields_c<vectorFieldT, tensorFieldT>)
     {
         // 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/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<polar_coordinates_t>(mito::tensor::e_0<2>);
+                return functions::constant<polar_coordinates_t>(mito::tensor::e_0<2>);
             }
 
             if constexpr (I == 1) {
                 // e_theta
-                return _r * fields::uniform_field<polar_coordinates_t>(mito::tensor::e_1<2>);
+                return _r * functions::constant<polar_coordinates_t>(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<polar_coordinates_t, 0>;
 
             // the function returning the constant e_00 tensor in 2D
-            constexpr auto e_rr = fields::uniform_field<polar_coordinates_t>(tensor::e_00<2>);
+            constexpr auto e_rr = functions::constant<polar_coordinates_t>(tensor::e_00<2>);
             // the function returning the constant e_11 tensor in 2D
-            constexpr auto e_tt = fields::uniform_field<polar_coordinates_t>(tensor::e_11<2>);
+            constexpr auto e_tt = functions::constant<polar_coordinates_t>(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<spherical_coordinates_t>(mito::tensor::e_0<3>);
+                return functions::constant<spherical_coordinates_t>(mito::tensor::e_0<3>);
             }
 
             if constexpr (I == 1) {
                 // return e_theta
-                return _r * fields::uniform_field<spherical_coordinates_t>(mito::tensor::e_1<3>);
+                return _r * functions::constant<spherical_coordinates_t>(mito::tensor::e_1<3>);
             }
 
             if constexpr (I == 2) {
                 // return e_phi
                 return _r * functions::sin(_theta)
-                     * fields::uniform_field<spherical_coordinates_t>(mito::tensor::e_2<3>);
+                     * functions::constant<spherical_coordinates_t>(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<spherical_coordinates_t, 1>;
 
             // the function returning the constant e_00 tensor in 3D
-            constexpr auto e_rr = fields::uniform_field<spherical_coordinates_t>(tensor::e_00<3>);
+            constexpr auto e_rr = functions::constant<spherical_coordinates_t>(tensor::e_00<3>);
             // the function returning the constant e_11 tensor in 3D
-            constexpr auto e_tt = fields::uniform_field<spherical_coordinates_t>(tensor::e_11<3>);
+            constexpr auto e_tt = functions::constant<spherical_coordinates_t>(tensor::e_11<3>);
             // the function returning the constant e_11 tensor in 3D
-            constexpr auto e_pp = fields::uniform_field<spherical_coordinates_t>(tensor::e_22<3>);
+            constexpr auto e_pp = functions::constant<spherical_coordinates_t>(tensor::e_22<3>);
 
             // return the metric field in spherical coordinates
             return e_rr + functions::pow<2>(r) * e_tt
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/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
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 <memory>
 #include <type_traits>
 #include <cstdint>
+#include <string>
+#include <typeinfo>
+#include <cxxabi.h>
 
 // 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 <typename T>
     using base_type = typename std::remove_cvref_t<T>;
+
+    // function returning the type name of {T} as a string
+    template <typename T>
+    std::string type_name()
+    {
+        std::string result;
+
+        // Add qualifiers manually
+        if (std::is_const_v<T>)
+            result += "const ";
+        if (std::is_volatile_v<T>)
+            result += "volatile ";
+        if (std::is_lvalue_reference_v<T>)
+            result += "& ";
+        if (std::is_rvalue_reference_v<T>)
+            result += "&& ";
+
+        using Base = std::remove_cv_t<std::remove_reference_t<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;
+    }
 }
 
 
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<cell_t>(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<mito::tensor::vector_t<3>>(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/discrete/point_field.cc b/tests/mito.lib/discrete/point_field.cc
new file mode 100644
index 000000000..ba3de1e8d
--- /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 <gtest/gtest.h>
+#include <mito/discrete.h>
+
+
+// pick a dimension in space
+constexpr auto D = 3;
+// a type for the field
+using vector_t = mito::tensor::vector_t<2>;
+// assemble the quadrature field
+using point_field_t = mito::discrete::point_field_t<D, vector_t>;
+
+// fetch the point cloud
+auto & cloud = mito::geometry::point_cloud<D>();
+// build a point
+auto point = cloud.point();
+
+
+void
+constFunction(const point_field_t & 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, PointFields)
+{
+    // a point field
+    auto field = point_field_t("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 <gtest/gtest.h>
-#include <mito/discrete.h>
-
-
-// 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<cell_type, Q, vector_type>;
-
-// 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
diff --git a/tests/mito.lib/fem/domain_field.cc b/tests/mito.lib/fem/domain_field.cc
new file mode 100644
index 000000000..ddd2564b3
--- /dev/null
+++ b/tests/mito.lib/fem/domain_field.cc
@@ -0,0 +1,68 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito/fem.h>
+
+
+// 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<coordinates_t>;
+// the x scalar field in 2D
+constexpr auto x = mito::functions::component<coordinates_t, 0>;
+// the y scalar field in 2D
+constexpr auto y = mito::functions::component<coordinates_t, 1>;
+
+
+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, { 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 });
+
+    // 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 }));
+
+    // 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;
+}
+
+
+// end of file
\ No newline at end of file
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 <gtest/gtest.h>
+#include <mito.h>
+
+
+// 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<degree, cell_t>;
+// the x scalar field in 2D
+constexpr auto x = mito::functions::component<coordinates_t, 0>;
+// the y scalar field in 2D
+constexpr auto y = mito::functions::component<coordinates_t, 1>;
+
+
+TEST(Fem, FemField)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system<coordinates_t>();
+
+    // read the mesh of a square in 2D
+    std::ifstream fileStream("square.summit");
+    auto mesh = mito::io::summit::reader<cell_t>(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<coordinates_t>;
+
+    // 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<finite_element_t>(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<mito::tensor::scalar_t>("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
diff --git a/tests/mito.lib/fields/calculus_identities.cc b/tests/mito.lib/fields/calculus_identities.cc
index ba7458e20..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<coordinates_t>(e_0);
-
-    // the field returning the constant {e_1} unit vector in 2D
-    constexpr auto e1 = mito::fields::uniform_field<coordinates_t>(e_1);
-
     // the sine function
     constexpr auto sin = mito::functions::sin;
 
@@ -42,14 +36,14 @@ TEST(Identities, DivGrad)
     constexpr auto x1 = mito::functions::component<coordinates_t, 1>;
 
     // 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<coordinates_t>({ pi_sixth, pi_fourth });
 
     // 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<coordinates_t, 1>;
 
     // 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/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<coordinates_t>(e_0);
-
-    // the field returning the constant {e_1} unit vector in 2D
-    constexpr auto e1 = mito::fields::uniform_field<coordinates_t>(e_1);
-
     // the sine function
     constexpr auto sin = mito::functions::sin;
 
@@ -48,7 +42,7 @@ TEST(Laplacian, VectorFields)
     constexpr auto x1 = mito::functions::component<coordinates_t, 1>;
 
     // 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<coordinates_t>({ pi_sixth, pi_fourth });
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<coordinates_t, 1>;
 
     // 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<coordinates_t>({ 0.0, 0.0 });
@@ -36,7 +36,7 @@ TEST(Fields, VectorFields)
     constexpr auto e1 = mito::functions::constant<coordinates_t>(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..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<coordinates_t, 0>;
 static constexpr auto x1 = mito::functions::component<coordinates_t, 1>;
@@ -17,30 +21,19 @@ static constexpr auto x1 = mito::functions::component<coordinates_t, 1>;
 
 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<decltype(f)>);
 
-    // the function returning the constant e0 unit vector in 2D
-    [[maybe_unused]] constexpr auto e0 =
-        mito::fields::uniform_field<coordinates_t>(mito::tensor::e_0<2>);
-
-    // the function returning the constant e1 unit vector in 2D
-    [[maybe_unused]] constexpr auto e1 =
-        mito::fields::uniform_field<coordinates_t>(mito::tensor::e_1<2>);
-
     // a vector field
-    constexpr auto g = mito::fields::field(cos * (e0 + e1));
+    constexpr auto g = f * (e_0 + e_1);
     // assert that {g} is a vector field
     static_assert(mito::fields::vector_field_c<decltype(g)>);
 
     // a second-order tensor field (2x2 identity tensor field in 3 dimensional space)
-    constexpr auto i = mito::fields::identity_tensor_field<coordinates_t, 3>;
+    constexpr auto i = mito::functions::identity<coordinates_t, 3>();
     // assert that {i} is a tensor field
     static_assert(mito::fields::tensor_field_c<decltype(i)>);
     // assert that {i} is a symmetric tensor field
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<coordinates_t>(e_0);
-
-    // the field returning the constant {e_1} unit vector in 2D
-    constexpr auto e1 = mito::fields::uniform_field<coordinates_t>(e_1);
-
     // the sine function
     constexpr auto sin = mito::functions::sin;
 
@@ -53,7 +47,7 @@ TEST(Gradient, NonSquare)
     constexpr auto x2 = mito::functions::component<coordinates_t, 2>;
 
     // 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<coordinates_t>({ 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 1bf087a30..2d5992f19 100644
--- a/tests/mito.lib/fields/polar_metric_field.cc
+++ b/tests/mito.lib/fields/polar_metric_field.cc
@@ -20,24 +20,24 @@ static constexpr auto r = mito::functions::component<coordinates_t, 0>;
 TEST(Fields, PolarCoordinates)
 {
     // the basis for vector fields (e_r and e_theta)
-    constexpr auto e_r = mito::fields::uniform_field<coordinates_t>(mito::tensor::e_0<2>);
-    constexpr auto e_t = r * mito::fields::uniform_field<coordinates_t>(mito::tensor::e_1<2>);
+    constexpr auto e_r = mito::functions::constant<coordinates_t>(mito::tensor::e_0<2>);
+    constexpr auto e_t = r * mito::functions::constant<coordinates_t>(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<coordinates_t>(mito::tensor::e_00<2>);
-    constexpr auto e_tt = mito::fields::uniform_field<coordinates_t>(mito::tensor::e_11<2>);
+    constexpr auto e_rr = mito::functions::constant<coordinates_t>(mito::tensor::e_00<2>);
+    constexpr auto e_tt = mito::functions::constant<coordinates_t>(mito::tensor::e_11<2>);
 
     // the metric field
     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::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 841175d71..9bb762d83 100644
--- a/tests/mito.lib/fields/spherical_metric_field.cc
+++ b/tests/mito.lib/fields/spherical_metric_field.cc
@@ -21,28 +21,28 @@ static constexpr auto t = mito::functions::component<coordinates_t, 1>;
 TEST(Manifolds, SphericalCoordinates)
 {
     // the basis for vector fields (e_r, e_theta, e_phi)
-    constexpr auto e_r = mito::fields::uniform_field<coordinates_t>(mito::tensor::e_0<3>);
-    constexpr auto e_t = r * mito::fields::uniform_field<coordinates_t>(mito::tensor::e_1<3>);
+    constexpr auto e_r = mito::functions::constant<coordinates_t>(mito::tensor::e_0<3>);
+    constexpr auto e_t = r * mito::functions::constant<coordinates_t>(mito::tensor::e_1<3>);
     constexpr auto e_p = r * mito::functions::sin(t)
-                       * mito::fields::uniform_field<coordinates_t>(mito::tensor::e_2<3>);
+                       * mito::functions::constant<coordinates_t>(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<coordinates_t>(mito::tensor::e_00<3>);
-    constexpr auto e_tt = mito::fields::uniform_field<coordinates_t>(mito::tensor::e_11<3>);
-    constexpr auto e_pp = mito::fields::uniform_field<coordinates_t>(mito::tensor::e_22<3>);
+    constexpr auto e_rr = mito::functions::constant<coordinates_t>(mito::tensor::e_00<3>);
+    constexpr auto e_tt = mito::functions::constant<coordinates_t>(mito::tensor::e_11<3>);
+    constexpr auto e_pp = mito::functions::constant<coordinates_t>(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;
 
     // 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::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 };
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<vector_t, 0>;
 
     // the function extracting the x_1 component
-    constexpr auto x1 = mito::functions::x<1, 2>;
+    constexpr auto x1 = mito::functions::component<vector_t, 1>;
 
     // a scalar-valued function
     constexpr auto function1 = mito::functions::cos(x0 * x1);
 
     // a scalar-valued function
-    constexpr auto function2 = 5.0 * mito::functions::one<mito::tensor::vector_t<2>>;
+    constexpr auto function2 = 5.0 * mito::functions::one<vector_t>;
 
     // 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::vector_t<2>>(mito::tensor::e_0<3>);
-    constexpr auto e1 = mito::functions::constant<mito::tensor::vector_t<2>>(mito::tensor::e_1<3>);
-    constexpr auto e2 = mito::functions::constant<mito::tensor::vector_t<2>>(mito::tensor::e_2<3>);
+    constexpr auto e0 = mito::functions::constant<vector_t>(mito::tensor::e_0<3>);
+    constexpr auto e1 = mito::functions::constant<vector_t>(mito::tensor::e_1<3>);
+    constexpr auto e2 = mito::functions::constant<vector_t>(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<vector_t, 0>;
     // 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<vector_t, 1>;
     // 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<vector_t, 2>;
+
     // 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::vector_t<2>>(mito::tensor::e_0<2>);
+    constexpr auto e0 = mito::functions::constant<vector_t>(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<vector_t, 0>;
+
     // the function extracting the x_1 component
-    constexpr auto x1 = mito::functions::x<1, 2>;
+    constexpr auto x1 = mito::functions::component<vector_t, 1>;
 
     // 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::vector_t<2>>(mito::tensor::e_0<2>);
+    constexpr auto e0 = mito::functions::constant<vector_t>(mito::tensor::e_0<2>);
 
     // the function returning the constant e1 unit vector in 2D
-    constexpr auto e1 = mito::functions::constant<mito::tensor::vector_t<2>>(mito::tensor::e_1<2>);
+    constexpr auto e1 = mito::functions::constant<vector_t>(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<vector_t, 0>;
 
     // 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<vector_t, 1>;
 
     // 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<vector_t, 0>;
 
     // 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<vector_t, 1>;
 
     // 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<vector_t, 0>;
     // the function extracting the x_1 component
-    constexpr auto x1 = mito::functions::x<1, 2>;
+    constexpr auto x1 = mito::functions::component<vector_t, 1>;
 
     // 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<vector_t, 0>;
     // 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<vector_t, 1>;
     // 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<vector_t, 0>;
     // the function extracting the x_1 component
-    constexpr auto x1 = mito::functions::x<1, 2>;
+    constexpr auto x1 = mito::functions::component<vector_t, 1>;
 
     // 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 <mito/functions.h>
 
 
-// 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<vector_t, 0>;
+static constexpr auto x1 = mito::functions::component<vector_t, 1>;
 
 // unit vectors in 3D
 static constexpr auto e0 = mito::tensor::e_0<3>;
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..6de6bdd25 100644
--- a/tests/mito.lib/integration/quadrature_flip_segment_3D.cc
+++ b/tests/mito.lib/integration/quadrature_flip_segment_3D.cc
@@ -39,11 +39,12 @@ 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<coordinates_t>(v_2);
-    constexpr auto normal_field_2 = mito::fields::uniform_field<coordinates_t>(v_3);
+    constexpr auto normal_field_1 = mito::functions::constant<coordinates_t>(v_2);
+    constexpr auto normal_field_2 = mito::functions::constant<coordinates_t>(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<mito::geometry::segment_t<3>>();
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<mito::quadrature::GAUSS, 2>(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<mito::quadrature::GAUSS, 2>(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<mito::tensor::vector_t<3>>(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<mito::tensor::vector_t<3>>(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..e2cda7bcc 100644
--- a/tests/mito.lib/manifolds/euclidean_gradient.cc
+++ b/tests/mito.lib/manifolds/euclidean_gradient.cc
@@ -29,13 +29,13 @@ 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);
+    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 4e009804f..1301af678 100644
--- a/tests/mito.lib/manifolds/polar_gradient.cc
+++ b/tests/mito.lib/manifolds/polar_gradient.cc
@@ -33,13 +33,13 @@ 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);
+    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 84e4e8995..dc39edfec 100644
--- a/tests/mito.lib/manifolds/spherical_gradient.cc
+++ b/tests/mito.lib/manifolds/spherical_gradient.cc
@@ -38,17 +38,16 @@ 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);
+    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 });
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<mito::geometry::triangle_t<3>>(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/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<spherical_coordinates_t>(mito::tensor::e_0<3>);
+    mito::functions::constant<spherical_coordinates_t>(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<coordinates_t>(normal_vector);
+    constexpr auto normal_field = mito::functions::constant<coordinates_t>(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_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<mito::quadrature::GAUSS, 2>(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..3702d2bcf 100644
--- a/tests/mito.lib/quadrature/quadrature_segment_3D.cc
+++ b/tests/mito.lib/quadrature/quadrature_segment_3D.cc
@@ -38,11 +38,11 @@ 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<coordinates_t>(v_2);
-    constexpr auto normal_field_2 = mito::fields::uniform_field<coordinates_t>(v_3);
+    constexpr auto normal_field_1 = mito::functions::constant<coordinates_t>(v_2);
+    constexpr auto normal_field_2 = mito::functions::constant<coordinates_t>(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<mito::geometry::segment_t<3>>();
diff --git a/tests/mito.lib/quadrature/quadrature_square_3D.cc b/tests/mito.lib/quadrature/quadrature_square_3D.cc
index 8c4967007..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<coordinates_t>(normal_vector);
+    constexpr auto normal_field = mito::functions::constant<coordinates_t>(normal_vector);
 
     // create a submanifold on {mesh} with the appropriate normal fields
     auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field);
@@ -77,7 +77,8 @@ TEST(Quadrature, Square)
     auto integrator = mito::quadrature::integrator<GAUSS, 2 /* degree of exactness */>(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<GAUSS, 2 /* degree of exactness */>(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<coordinates_t>(1.0));
+    auto f_one = mito::functions::constant<coordinates_t>(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