Template over mesh geometry type (#2603)

* Template over mesh geometry type * Improve type deduction * Doc fix * Simplification * Demo updates * More changes * Remove file * Updates * Test update * Updates * Get type deduction working * Updates * Test update * Update demos * Simplify ADIOS2 mesh extraction * ADIOS2 IO re-factoring * Template some IO functions * Move IO code * Doc fix * Fix linkage * Demo update * Get IO working for different mesh types * Uncomment code * Doc fix * Work on interpolation types * Remove file * Add missing implementation * Wrapper update * More transition * Updates * Updates * Tidy up * More changes * Tidy up * Refactoring * Template fix * Undo some changes * BB update * Fix for gcc * Fixes * gcc updates * Updates * More transition * More generic code * Clean up * Simplify * Fix demo * Fix doc error * Small improvements * Tidy up * Tidy up * Small improvements * Fix cpp/python error * Simplifications * Demo update * Small demo update * Updates * Template updates * Assembler updates * Tidy up * Generalise testing * Updates * Formatting fixes * Misc fixes * Small doc update * More concepts * concepts updates --------- Co-authored-by: Garth Wells <garth@gnw20mac.eng.cam.ac.uk>
FEniCS · Mar 28, 2023 · 71e3e3f · 71e3e3f
1 parent 8fc8e1c
commit 71e3e3f
Show file tree

Hide file tree

Showing 89 changed files with 6,829 additions and 6,939 deletions.
diff --git a/cpp/cmake/templates/DOLFINXConfig.cmake.in b/cpp/cmake/templates/DOLFINXConfig.cmake.in
@@ -75,6 +75,10 @@ if(@SLEPC_FOUND@)
   endif()
 endif()
 
+if(@ADIOS2_FOUND@)
+  find_dependency(ADIOS2 2.8.1)
+endif()
+
 if(NOT TARGET dolfinx)
   include("${CMAKE_CURRENT_LIST_DIR}/DOLFINXTargets.cmake")
 endif()

diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp
@@ -123,6 +123,7 @@
 #include "biharmonic.h"
 #include <cmath>
 #include <dolfinx.h>
+#include <dolfinx/common/types.h>
 #include <dolfinx/fem/Constant.h>
 #include <dolfinx/fem/petsc.h>
 #include <numbers>
@@ -149,7 +150,7 @@ int main(int argc, char* argv[])
   {
     // Create mesh
     auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
-    auto mesh = std::make_shared<mesh::Mesh>(
+    auto mesh = std::make_shared<mesh::Mesh<double>>(
         mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}},
                                {32, 32}, mesh::CellType::triangle, part));
 
@@ -159,7 +160,7 @@ int main(int argc, char* argv[])
     // .. code-block:: cpp
 
     // Create function space
-    auto V = std::make_shared<fem::FunctionSpace>(
+    auto V = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(functionspace_form_biharmonic_a, "u", mesh));
 
     // The source function ::math:`f` and the penalty term
@@ -219,7 +220,8 @@ int main(int argc, char* argv[])
           }
           return marker;
         });
-    const auto bdofs = fem::locate_dofs_topological({*V}, 1, facets);
+    const auto bdofs = fem::locate_dofs_topological(
+        V->mesh()->topology_mutable(), *V->dofmap(), 1, facets);
     auto bc = std::make_shared<const fem::DirichletBC<T>>(0.0, bdofs, V);
 
     // Now, we have specified the variational forms and can consider the
@@ -249,9 +251,9 @@ int main(int argc, char* argv[])
 
     b.set(0.0);
     fem::assemble_vector(b.mutable_array(), *L);
-    fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, T(1.0));
+    fem::apply_lifting<T, double>(b.mutable_array(), {a}, {{bc}}, {}, T(1.0));
     b.scatter_rev(std::plus<T>());
-    fem::set_bc(b.mutable_array(), {bc});
+    fem::set_bc<T, double>(b.mutable_array(), {bc});
 
     la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
     la::petsc::options::set("ksp_type", "preonly");

diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
@@ -127,13 +127,14 @@ int main(int argc, char* argv[])
     // .. code-block:: cpp
 
     // Create mesh and define function space
-    auto mesh = std::make_shared<mesh::Mesh>(
+    auto mesh = std::make_shared<mesh::Mesh<double>>(
         mesh::create_box(MPI_COMM_WORLD, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}},
                          {10, 10, 10}, mesh::CellType::tetrahedron,
                          mesh::create_cell_partitioner(mesh::GhostMode::none)));
 
-    auto V = std::make_shared<fem::FunctionSpace>(fem::create_functionspace(
-        functionspace_form_hyperelasticity_F_form, "u", mesh));
+    auto V = std::make_shared<fem::FunctionSpace<double>>(
+        fem::create_functionspace(functionspace_form_hyperelasticity_F_form,
+                                  "u", mesh));
 
     // Define solution function
     auto u = std::make_shared<fem::Function<T>>(V);
@@ -178,7 +179,7 @@ int main(int argc, char* argv[])
 
     // Create Dirichlet boundary conditions
     auto bdofs_left = fem::locate_dofs_geometrical(
-        {*V},
+        *V,
         [](auto x)
         {
           constexpr double eps = 1.0e-8;
@@ -191,7 +192,7 @@ int main(int argc, char* argv[])
           return marker;
         });
     auto bdofs_right = fem::locate_dofs_geometrical(
-        {*V},
+        *V,
         [](auto x)
         {
           constexpr double eps = 1.0e-8;
@@ -228,10 +229,11 @@ int main(int argc, char* argv[])
     const basix::FiniteElement S_element = basix::create_element(
         family, cell_type, k, basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, discontinuous);
-    auto S = std::make_shared<fem::FunctionSpace>(fem::create_functionspace(
-        mesh, S_element, pow(mesh->geometry().dim(), 2)));
+    auto S = std::make_shared<fem::FunctionSpace<double>>(
+        fem::create_functionspace(mesh, S_element,
+                                  pow(mesh->geometry().dim(), 2)));
 
-    const auto sigma_expression = fem::create_expression<T>(
+    auto sigma_expression = fem::create_expression<T, double>(
         *expression_hyperelasticity_sigma, {{"u", u}}, {}, mesh);
 
     auto sigma = fem::Function<T>(S);

diff --git a/cpp/demo/interpolation-io/main.cpp b/cpp/demo/interpolation-io/main.cpp
@@ -1,11 +1,12 @@
-// Copyright (C) 2022 Garth N. Wells
+// Copyright (C) 2022-2023 Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
 // SPDX-License-Identifier:    LGPL-3.0-or-later
 
 #include <basix/finite-element.h>
 #include <cmath>
+#include <concepts>
 #include <dolfinx/common/log.h>
 #include <dolfinx/fem/FiniteElement.h>
 #include <dolfinx/fem/FunctionSpace.h>
@@ -15,6 +16,7 @@
 #include <dolfinx/mesh/Mesh.h>
 #include <dolfinx/mesh/cell_types.h>
 #include <dolfinx/mesh/generation.h>
+#include <dolfinx/mesh/utils.h>
 #include <filesystem>
 #include <mpi.h>
 #include <numbers>
@@ -24,9 +26,9 @@ using namespace dolfinx;
 // This function interpolations a function is a finite element space and
 // outputs the finite element function to a VTK file for visualisation.
 // It also shows how to create a finite element using Basix.
-template <typename T>
-void interpolate_scalar(std::shared_ptr<mesh::Mesh> mesh,
-                        std::filesystem::path filename)
+template <typename T, std::floating_point U>
+void interpolate_scalar(std::shared_ptr<mesh::Mesh<U>> mesh,
+                        [[maybe_unused]] std::filesystem::path filename)
 {
   // Create a Basix continuous Lagrange element of degree 1
   basix::FiniteElement e = basix::create_element(
@@ -36,7 +38,7 @@ void interpolate_scalar(std::shared_ptr<mesh::Mesh> mesh,
       basix::element::dpc_variant::unset, false);
 
   // Create a scalar function space
-  auto V = std::make_shared<fem::FunctionSpace>(
+  auto V = std::make_shared<fem::FunctionSpace<U>>(
       fem::create_functionspace(mesh, e, 1));
 
   // Create a finite element Function
@@ -53,18 +55,21 @@ void interpolate_scalar(std::shared_ptr<mesh::Mesh> mesh,
         return {f, {f.size()}};
       });
 
-  // Write the function to a VTK file for visualisation, e.g. using
+#ifdef HAS_ADIOS2
+  // Write the function to a VTX file for visualisation, e.g. using
   // ParaView
-  io::VTKFile file(mesh->comm(), filename.replace_extension("pvd"), "w");
-  file.write<T>({*u}, 0.0);
+  io::VTXWriter<U> outfile(mesh->comm(), filename.replace_extension("bp"), {u});
+  outfile.write(0.0);
+  outfile.close();
+#endif
 }
 
 // This function interpolations a function is a H(curl) finite element
 // space. To visualise the function, it interpolates the H(curl) finite
 // element function in a discontinuous Lagrange space and outputs the
 // Lagrange finite element function to a VTX file for visualisation.
-template <typename T>
-void interpolate_nedelec(std::shared_ptr<mesh::Mesh> mesh,
+template <typename T, std::floating_point U>
+void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
                          [[maybe_unused]] std::filesystem::path filename)
 {
   // Create a Basix Nedelec (first kind) element of degree 2 (dim=6 on triangle)
@@ -75,7 +80,7 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh> mesh,
       basix::element::dpc_variant::unset, false);
 
   // Create a Nedelec function space
-  auto V = std::make_shared<fem::FunctionSpace>(
+  auto V = std::make_shared<fem::FunctionSpace<U>>(
       fem::create_functionspace(mesh, e, 1));
 
   // Create a Nedelec finite element Function
@@ -145,7 +150,7 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh> mesh,
       basix::element::dpc_variant::unset, true);
 
   // Create a function space
-  auto V_l = std::make_shared<fem::FunctionSpace>(
+  auto V_l = std::make_shared<fem::FunctionSpace<U>>(
       fem::create_functionspace(mesh, e_l, 2));
 
   auto u_l = std::make_shared<fem::Function<T>>(V_l);
@@ -154,13 +159,14 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh> mesh,
   // space:
   u_l->interpolate(*u);
 
-  // Output the discontinuous Lagrange space in VTK format. When
+  // Output the discontinuous Lagrange space in VTX format. When
   // plotting the x0 component the field will appear discontinuous at x0
   // = 0.5 (jump in the normal component between cells) and the x1
   // component will appear continuous (continuous tangent component
   // between cells).
 #ifdef HAS_ADIOS2
-  io::VTXWriter outfile(mesh->comm(), filename.replace_extension("bp"), {u_l});
+  io::VTXWriter<U> outfile(mesh->comm(), filename.replace_extension("bp"),
+                           {u_l});
   outfile.write(0.0);
   outfile.close();
 #endif
@@ -173,27 +179,41 @@ int main(int argc, char* argv[])
   dolfinx::init_logging(argc, argv);
   MPI_Init(&argc, &argv);
 
-  // The main body of the function is scoped with the curly braces to
-  // ensure that all objects that depend on an MPI communicator are
-  // destroyed before MPI is finalised at the end of this function.
+  // The main body of the function is scoped to ensure that all objects
+  // that depend on an MPI communicator are destroyed before MPI is
+  // finalised at the end of this function.
   {
-    // Create a mesh. For what comes later in this demo we need to
+    // Create meshes. For what comes later in this demo we need to
     // ensure that a boundary between cells is located at x0=0.5
-    auto mesh = std::make_shared<mesh::Mesh>(mesh::create_rectangle(
-        MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {32, 4},
-        mesh::CellType::triangle,
-        mesh::create_cell_partitioner(mesh::GhostMode::none)));
+
+    // Create mesh using float for geometry coordinates
+    auto mesh0
+        = std::make_shared<mesh::Mesh<float>>(mesh::create_rectangle<float>(
+            MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {32, 4},
+            mesh::CellType::triangle,
+            mesh::create_cell_partitioner(mesh::GhostMode::none)));
+
+    // Create mesh using double for geometry coordinates
+    auto mesh1
+        = std::make_shared<mesh::Mesh<double>>(mesh::create_rectangle<double>(
+            MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, {32, 4},
+            mesh::CellType::triangle,
+            mesh::create_cell_partitioner(mesh::GhostMode::none)));
 
     // Interpolate a function in a scalar Lagrange space and output the
-    // result to file for visualisation
-    interpolate_scalar<double>(mesh, "u");
-    interpolate_scalar<std::complex<double>>(mesh, "u_complex");
+    // result to file for visualisation using different types
+    interpolate_scalar<float>(mesh0, "u32");
+    interpolate_scalar<double>(mesh1, "u64");
+    interpolate_scalar<std::complex<float>>(mesh0, "u_complex64");
+    interpolate_scalar<std::complex<double>>(mesh1, "u_complex128");
 
     // Interpolate a function in a H(curl) finite element space, and
     // then interpolate the H(curl) function in a discontinuous Lagrange
-    // space for visualisation
-    interpolate_nedelec<double>(mesh, "u_nedelec");
-    interpolate_nedelec<std::complex<double>>(mesh, "u_nedelec_complex");
+    // space for visualisation using different types
+    interpolate_nedelec<float>(mesh0, "u_nedelec32");
+    interpolate_nedelec<double>(mesh1, "u_nedelec64");
+    interpolate_nedelec<std::complex<float>>(mesh0, "u_nedelec_complex64");
+    interpolate_nedelec<std::complex<double>>(mesh1, "u_nedelec_complex12");
   }
 
   MPI_Finalize();

diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -23,24 +23,25 @@ int main(int argc, char* argv[])
     MPI_Comm comm{MPI_COMM_WORLD};
 
     // Create a tetrahedral mesh
-    auto mesh_tet = std::make_shared<mesh::Mesh>(
+    auto mesh_tet = std::make_shared<mesh::Mesh<double>>(
         mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {20, 20, 20},
                          mesh::CellType::tetrahedron));
+
     // Create a hexahedral mesh
-    auto mesh_hex = std::make_shared<mesh::Mesh>(
+    auto mesh_hex = std::make_shared<mesh::Mesh<double>>(
         mesh::create_box(comm, {{{0, 0, 0}, {1, 1, 1}}}, {15, 15, 15},
                          mesh::CellType::hexahedron));
 
     basix::FiniteElement element_tet = basix::element::create_lagrange(
         mesh::cell_type_to_basix_type(mesh_tet->topology().cell_type()), 1,
         basix::element::lagrange_variant::equispaced, false);
-    auto V_tet = std::make_shared<fem::FunctionSpace>(
+    auto V_tet = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(mesh_tet, element_tet, 3));
 
     basix::FiniteElement element_hex = basix::element::create_lagrange(
         mesh::cell_type_to_basix_type(mesh_hex->topology().cell_type()), 2,
         basix::element::lagrange_variant::equispaced, false);
-    auto V_hex = std::make_shared<fem::FunctionSpace>(
+    auto V_hex = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(mesh_hex, element_hex, 3));
 
     auto u_tet = std::make_shared<fem::Function<T>>(V_tet);
@@ -72,10 +73,10 @@ int main(int argc, char* argv[])
     u_hex->interpolate(*u_tet, nmm_interpolation_data);
 
 #ifdef HAS_ADIOS2
-    io::VTXWriter write_tet(mesh_tet->comm(), "u_tet.vtx", {u_tet});
+    io::VTXWriter<double> write_tet(mesh_tet->comm(), "u_tet.bp", {u_tet});
     write_tet.write(0.0);
 
-    io::VTXWriter write_hex(mesh_hex->comm(), "u_hex.vtx", {u_hex});
+    io::VTXWriter<double> write_hex(mesh_hex->comm(), "u_hex.bp", {u_hex});
     write_hex.write(0.0);
 #endif
   }

diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp
@@ -121,11 +121,11 @@ int main(int argc, char* argv[])
   {
     // Create mesh and function space
     auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
-    auto mesh = std::make_shared<mesh::Mesh>(
+    auto mesh = std::make_shared<mesh::Mesh<double>>(
         mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 1.0}}},
                                {32, 16}, mesh::CellType::triangle, part));
 
-    auto V = std::make_shared<fem::FunctionSpace>(
+    auto V = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));
 
     // Next, we define the variational formulation by initializing the
@@ -166,17 +166,19 @@ int main(int argc, char* argv[])
         *mesh, 1,
         [](auto x)
         {
-          constexpr double eps = 1.0e-8;
+          using U = typename decltype(x)::value_type;
+          constexpr U eps = 1.0e-8;
           std::vector<std::int8_t> marker(x.extent(1), false);
           for (std::size_t p = 0; p < x.extent(1); ++p)
           {
-            double x0 = x(0, p);
+            auto x0 = x(0, p);
             if (std::abs(x0) < eps or std::abs(x0 - 2) < eps)
               marker[p] = true;
           }
           return marker;
         });
-    const auto bdofs = fem::locate_dofs_topological({*V}, 1, facets);
+    const auto bdofs = fem::locate_dofs_topological(
+        V->mesh()->topology_mutable(), *V->dofmap(), 1, facets);
     auto bc = std::make_shared<const fem::DirichletBC<T>>(0.0, bdofs, V);
 
     f->interpolate(
@@ -185,8 +187,8 @@ int main(int argc, char* argv[])
           std::vector<T> f;
           for (std::size_t p = 0; p < x.extent(1); ++p)
           {
-            double dx = (x(0, p) - 0.5) * (x(0, p) - 0.5);
-            double dy = (x(1, p) - 0.5) * (x(1, p) - 0.5);
+            auto dx = (x(0, p) - 0.5) * (x(0, p) - 0.5);
+            auto dy = (x(1, p) - 0.5) * (x(1, p) - 0.5);
             f.push_back(10 * std::exp(-(dx + dy) / 0.02));
           }
 
@@ -229,9 +231,9 @@ int main(int argc, char* argv[])
 
     b.set(0.0);
     fem::assemble_vector(b.mutable_array(), *L);
-    fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, T(1));
+    fem::apply_lifting<T, double>(b.mutable_array(), {a}, {{bc}}, {}, T(1));
     b.scatter_rev(std::plus<T>());
-    fem::set_bc(b.mutable_array(), {bc});
+    fem::set_bc<T, double>(b.mutable_array(), {bc});
 
     la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
     la::petsc::options::set("ksp_type", "preonly");