Add C++ demo that shows mixed elements and consistency fixes for elem…

…ent spaces (#3500) * Start on C++ mixed method example * Tidy up demo * Undo change * Lint fix * Tidy up * Simplification * Fix forwarding * Remove redundant function argument * Formatting * Wrapper update * Consistency fixes * Fix * Undo changes * Consistency fixes * Undo some changes * Updates * Update * Small fix * More updates * Enable demo * Updates * Improve readability * Small updates * Add CMake file * Solver fix * Change solver * Try mumps * Verbose test * Action update * Print matrix norm * Test image path update * Update solve * Docker update * CI update * Demo update * Split create FunctionSpace function * Updates * Update * Revert * Revert * Updates * Tidy up * Small fixes * Small edits * Simplifications * Logic fix * Attempt fix * Revert * Change solver * Updates * Update * Simplify * Simplify * Update test * Test update * Enable demo * Tidy up * Simplification * Simplification * Demo update * Doc fix * Updates * Wrapper fix * Test updates * Quad element update * Update test * Test fix * Updates * Test update * Remove value_shape from function space * Lint * Constructor fix * Tidy * Tidy up * Tidy * Work on docs * Doc updates * Doc updates * Simplification * Revert some lines * Add comment * Small fix * Tidy up * Simplify * Work on doc * Type fix * Doc updates * Update docs * Doc improvements * Type fix * Switch solver for PETSc 64-bit int * Small fix * Small change
FEniCS · Nov 15, 2024 · 594f190 · 594f190
1 parent 5323846
commit 594f190
Show file tree

Hide file tree

Showing 53 changed files with 1,312 additions and 682 deletions.
diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml
@@ -207,7 +207,7 @@ jobs:
           cmake --build build/demo --parallel 3
           cd build/demo
           ctest -V -R demo -R serial
-          ctest -V -R demo -R mpi_2
+          ctest --output-on-failure -V -R demo -R mpi_2
 
       - name: Install Python build dependencies
         run: pip install -r python/build-requirements.txt

diff --git a/cpp/demo/CMakeLists.txt b/cpp/demo/CMakeLists.txt
@@ -19,8 +19,9 @@ endmacro(add_demo_subdirectory)
 add_demo_subdirectory(biharmonic)
 add_demo_subdirectory(codim_0_assembly)
 add_demo_subdirectory(custom_kernel)
-add_demo_subdirectory(poisson)
-add_demo_subdirectory(poisson_matrix_free)
 add_demo_subdirectory(hyperelasticity)
 add_demo_subdirectory(interpolation-io)
 add_demo_subdirectory(interpolation_different_meshes)
+add_demo_subdirectory(mixed_poisson)
+add_demo_subdirectory(poisson)
+add_demo_subdirectory(poisson_matrix_free)
diff --git a/cpp/demo/biharmonic/main.cpp b/cpp/demo/biharmonic/main.cpp
@@ -134,6 +134,7 @@
 #include <dolfinx/common/types.h>
 #include <dolfinx/fem/Constant.h>
 #include <dolfinx/fem/petsc.h>
+#include <memory>
 #include <numbers>
 #include <utility>
 #include <vector>
@@ -169,8 +170,9 @@ int main(int argc, char* argv[])
         basix::element::dpc_variant::unset, false);
 
     //  Create function space
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(element)));
 
     // The source function $f$ and the penalty term $\alpha$ are
     // declared:
@@ -190,10 +192,10 @@ int main(int argc, char* argv[])
     auto alpha = std::make_shared<fem::Constant<T>>(8.0);
 
     //  Define variational forms
-    auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_biharmonic_a, {V, V}, {}, {{"alpha", alpha}}, {}, {}));
-    auto L = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}, {}));
+    fem::Form<T> a = fem::create_form<T>(*form_biharmonic_a, {V, V}, {},
+                                         {{"alpha", alpha}}, {}, {});
+    fem::Form<T> L
+        = fem::create_form<T>(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}, {});
 
     //  Now, the Dirichlet boundary condition ($u = 0$) can be
     //  created using the class {cpp:class}`DirichletBC`. A
@@ -210,7 +212,7 @@ int main(int argc, char* argv[])
     auto facets = mesh::exterior_facet_indices(*mesh->topology());
     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);
+    fem::DirichletBC<T> bc(0.0, bdofs, V);
 
     //  Now, we have specified the variational forms and can consider
     //  the solution of the variational problem. First, we need to
@@ -221,13 +223,13 @@ int main(int argc, char* argv[])
 
     //  Compute solution
     fem::Function<T> u(V);
-    auto A = la::petsc::Matrix(fem::petsc::create_matrix(*a), false);
-    la::Vector<T> b(L->function_spaces()[0]->dofmap()->index_map,
-                    L->function_spaces()[0]->dofmap()->index_map_bs());
+    auto A = la::petsc::Matrix(fem::petsc::create_matrix(a), false);
+    la::Vector<T> b(L.function_spaces()[0]->dofmap()->index_map,
+                    L.function_spaces()[0]->dofmap()->index_map_bs());
 
     MatZeroEntries(A.mat());
     fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A.mat(), ADD_VALUES),
-                         *a, {bc});
+                         a, {bc});
     MatAssemblyBegin(A.mat(), MAT_FLUSH_ASSEMBLY);
     MatAssemblyEnd(A.mat(), MAT_FLUSH_ASSEMBLY);
     fem::set_diagonal<T>(la::petsc::Matrix::set_fn(A.mat(), INSERT_VALUES), *V,
@@ -236,10 +238,10 @@ int main(int argc, char* argv[])
     MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY);
 
     b.set(0.0);
-    fem::assemble_vector(b.mutable_array(), *L);
+    fem::assemble_vector(b.mutable_array(), L);
     fem::apply_lifting<T, U>(b.mutable_array(), {a}, {{bc}}, {}, T(1.0));
     b.scatter_rev(std::plus<T>());
-    bc->set(b.mutable_array(), std::nullopt);
+    bc.set(b.mutable_array(), std::nullopt);
 
     la::petsc::KrylovSolver lu(MPI_COMM_WORLD);
     la::petsc::options::set("ksp_type", "preonly");

diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp
@@ -41,8 +41,9 @@ int main(int argc, char* argv[])
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(element)));
 
     // Next we find all cells of the mesh with y<0.5
     const int tdim = mesh->topology()->dim();
@@ -85,8 +86,9 @@ int main(int argc, char* argv[])
     }
 
     // We create the function space used for the trial space
-    auto W = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(submesh, element, {}));
+    auto W
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            submesh, std::make_shared<fem::FiniteElement<U>>(element)));
 
     // A mixed-domain form has functions defined over different meshes.
     // The mesh associated with the measure (dx, ds, etc.) is called the
@@ -108,33 +110,32 @@ int main(int argc, char* argv[])
     // Next we compute the integration entities on the integration
     // domain `mesh`
     std::vector<std::int32_t> integration_entities
-        = fem::compute_integration_domains(fem::IntegralType::cell,
-                                           *mesh->topology(),
-                                           cell_marker.find(2), tdim);
+        = fem::compute_integration_domains(
+            fem::IntegralType::cell, *mesh->topology(), cell_marker.find(2));
     std::map<
         fem::IntegralType,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
         subdomain_data
         = {{fem::IntegralType::cell, {{3, integration_entities}}}};
 
     // We can now create the bilinear form
-    auto a_mixed = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_mixed_codim0_a_mixed, {V, W}, {}, {},
-                            subdomain_data, entity_maps, V->mesh()));
+    fem::Form<T> a_mixed
+        = fem::create_form<T>(*form_mixed_codim0_a_mixed, {V, W}, {}, {},
+                              subdomain_data, entity_maps, V->mesh());
 
-    la::SparsityPattern sp_mixed = fem::create_sparsity_pattern(*a_mixed);
+    la::SparsityPattern sp_mixed = fem::create_sparsity_pattern(a_mixed);
     sp_mixed.finalize();
     la::MatrixCSR<PetscScalar> A_mixed(sp_mixed);
-    fem::assemble_matrix(A_mixed.mat_add_values(), *a_mixed, {});
+    fem::assemble_matrix(A_mixed.mat_add_values(), a_mixed, {});
     A_mixed.scatter_rev();
 
-    auto a = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_mixed_codim0_a, {W, W}, {}, {}, {}, {}));
+    fem::Form<T> a
+        = fem::create_form<T>(*form_mixed_codim0_a, {W, W}, {}, {}, {}, {});
 
-    la::SparsityPattern sp = fem::create_sparsity_pattern(*a);
+    la::SparsityPattern sp = fem::create_sparsity_pattern(a);
     sp.finalize();
     la::MatrixCSR<PetscScalar> A(sp);
-    fem::assemble_matrix(A.mat_add_values(), *a, {});
+    fem::assemble_matrix(A.mat_add_values(), a, {});
     A.scatter_rev();
 
     std::vector<T> A_mixed_flattened = A_mixed.to_dense();

diff --git a/cpp/demo/custom_kernel/main.cpp b/cpp/demo/custom_kernel/main.cpp
@@ -209,8 +209,8 @@ void assemble(MPI_Comm comm)
   mdspand_t<const T, 2> X(X_b.data(), weights.size(), 2);
 
   // Create a scalar function space
-  auto V = std::make_shared<fem::FunctionSpace<T>>(
-      fem::create_functionspace(mesh, e));
+  auto V = std::make_shared<fem::FunctionSpace<T>>(fem::create_functionspace<T>(
+      mesh, std::make_shared<fem::FiniteElement<T>>(e)));
 
   // Build list of cells to assembler over (all cells owned by this
   // rank)

diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp
@@ -43,16 +43,15 @@ class HyperElasticProblem
 {
 public:
   /// Constructor
-  HyperElasticProblem(
-      std::shared_ptr<fem::Form<T>> L, std::shared_ptr<fem::Form<T>> J,
-      std::vector<std::shared_ptr<const fem::DirichletBC<T>>> bcs)
-      : _l(L), _j(J), _bcs(bcs),
-        _b(L->function_spaces()[0]->dofmap()->index_map,
-           L->function_spaces()[0]->dofmap()->index_map_bs()),
-        _matA(la::petsc::Matrix(fem::petsc::create_matrix(*J, "aij"), false))
+  HyperElasticProblem(fem::Form<T>& L, fem::Form<T>& J,
+                      const std::vector<fem::DirichletBC<T>>& bcs)
+      : _l(L), _j(J), _bcs(bcs.begin(), bcs.end()),
+        _b(L.function_spaces()[0]->dofmap()->index_map,
+           L.function_spaces()[0]->dofmap()->index_map_bs()),
+        _matA(la::petsc::Matrix(fem::petsc::create_matrix(J, "aij"), false))
   {
-    auto map = L->function_spaces()[0]->dofmap()->index_map;
-    const int bs = L->function_spaces()[0]->dofmap()->index_map_bs();
+    auto map = L.function_spaces()[0]->dofmap()->index_map;
+    const int bs = L.function_spaces()[0]->dofmap()->index_map_bs();
     std::int32_t size_local = bs * map->size_local();
 
     std::vector<PetscInt> ghosts(map->ghosts().begin(), map->ghosts().end());
@@ -88,7 +87,7 @@ class HyperElasticProblem
       // Assemble b and update ghosts
       std::span b(_b.mutable_array());
       std::ranges::fill(b, 0);
-      fem::assemble_vector<T>(b, *_l);
+      fem::assemble_vector<T>(b, _l);
       VecGhostUpdateBegin(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
       VecGhostUpdateEnd(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
 
@@ -100,7 +99,7 @@ class HyperElasticProblem
       const T* _x = nullptr;
       VecGetArrayRead(x_local, &_x);
       std::ranges::for_each(_bcs, [b, x = std::span(_x, n)](auto& bc)
-                            { bc->set(b, x, -1); });
+                            { bc.get().set(b, x, -1); });
       VecRestoreArrayRead(x_local, &_x);
     };
   }
@@ -111,12 +110,12 @@ class HyperElasticProblem
     return [&](const Vec, Mat A)
     {
       MatZeroEntries(A);
-      fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), *_j,
+      fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), _j,
                            _bcs);
       MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
       MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
       fem::set_diagonal(la::petsc::Matrix::set_fn(A, INSERT_VALUES),
-                        *_j->function_spaces()[0], _bcs);
+                        *_j.function_spaces()[0], _bcs);
       MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
       MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
     };
@@ -129,8 +128,9 @@ class HyperElasticProblem
   Mat matrix() { return _matA.mat(); }
 
 private:
-  std::shared_ptr<fem::Form<T>> _l, _j;
-  std::vector<std::shared_ptr<const fem::DirichletBC<T>>> _bcs;
+  fem::Form<T>& _l;
+  fem::Form<T>& _j;
+  std::vector<std::reference_wrapper<const fem::DirichletBC<T>>> _bcs;
   la::Vector<T> _b;
   Vec _b_petsc = nullptr;
   la::petsc::Matrix _matA;
@@ -166,20 +166,22 @@ int main(int argc, char* argv[])
         basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, false);
 
-    auto V = std::make_shared<fem::FunctionSpace<U>>(
-        fem::create_functionspace(mesh, element, {3}));
+    auto V
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(
+                      element, std::vector<std::size_t>{3})));
 
     auto B = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
     auto traction = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
 
     // Define solution function
     auto u = std::make_shared<fem::Function<T>>(V);
-    auto a = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
-                            {{"B", B}, {"T", traction}}, {}, {}));
-    auto L = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
-                            {{"B", B}, {"T", traction}}, {}, {}));
+    fem::Form<T> a
+        = fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
+                              {{"B", B}, {"T", traction}}, {}, {});
+    fem::Form<T> L
+        = fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
+                              {{"B", B}, {"T", traction}}, {}, {});
 
     auto u_rotation = std::make_shared<fem::Function<T>>(V);
     u_rotation->interpolate(
@@ -243,10 +245,9 @@ int main(int argc, char* argv[])
           }
           return marker;
         });
-    std::vector bcs = {
-        std::make_shared<const fem::DirichletBC<T>>(std::vector<T>{0, 0, 0},
-                                                    bdofs_left, V),
-        std::make_shared<const fem::DirichletBC<T>>(u_rotation, bdofs_right)};
+    std::vector bcs
+        = {fem::DirichletBC<T>(std::vector<T>{0, 0, 0}, bdofs_left, V),
+           fem::DirichletBC<T>(u_rotation, bdofs_right)};
 
     HyperElasticProblem problem(L, a, bcs);
     nls::petsc::NewtonSolver newton_solver(mesh->comm());
@@ -262,6 +263,9 @@ int main(int argc, char* argv[])
 
     // Compute Cauchy stress. Construct appropriate Basix element for
     // stress.
+    fem::Expression sigma_expression = fem::create_expression<T, U>(
+        *expression_hyperelasticity_sigma, {{"u", u}}, {});
+
     constexpr auto family = basix::element::family::P;
     auto cell_type
         = mesh::cell_type_to_basix_type(mesh->topology()->cell_type());
@@ -270,12 +274,12 @@ int main(int argc, char* argv[])
     basix::FiniteElement S_element = basix::create_element<U>(
         family, cell_type, k, basix::element::lagrange_variant::unset,
         basix::element::dpc_variant::unset, discontinuous);
-    auto S = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
-        mesh, S_element, std::vector<std::size_t>{3, 3}));
-    auto sigma_expression = fem::create_expression<T, U>(
-        *expression_hyperelasticity_sigma, {{"u", u}}, {});
+    auto S
+        = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+            mesh, std::make_shared<fem::FiniteElement<U>>(
+                      S_element, std::vector<std::size_t>{3, 3})));
 
-    auto sigma = fem::Function<T>(S);
+    fem::Function<T> sigma(S);
     sigma.name = "cauchy_stress";
     sigma.interpolate(sigma_expression);
 

diff --git a/cpp/demo/interpolation-io/main.cpp b/cpp/demo/interpolation-io/main.cpp
@@ -46,8 +46,8 @@ void interpolate_scalar(std::shared_ptr<mesh::Mesh<U>> mesh,
       basix::element::dpc_variant::unset, false);
 
   // Create a scalar function space
-  auto V = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, e));
+  auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+      mesh, std::make_shared<fem::FiniteElement<U>>(e)));
 
   // Create a finite element Function
   auto u = std::make_shared<fem::Function<T>>(V);
@@ -98,8 +98,8 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
       basix::element::dpc_variant::unset, false);
 
   // Create a Nedelec function space
-  auto V = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, e));
+  auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+      mesh, std::make_shared<fem::FiniteElement<U>>(e)));
 
   // Create a Nedelec finite element Function
   auto u = std::make_shared<fem::Function<T>>(V);
@@ -114,7 +114,7 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
   // the Nedelec space when there are cell edges aligned to x0 = 0.5.
 
   // Find cells with all vertices satisfying (0) x0 <= 0.5 and (1) x0 >= 0.5
-  auto cells0
+  std::vector cells0
       = mesh::locate_entities(*mesh, 2,
                               [](auto x)
                               {
@@ -123,7 +123,7 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
                                   marked.push_back(x(0, i) <= 0.5);
                                 return marked;
                               });
-  auto cells1
+  std::vector cells1
       = mesh::locate_entities(*mesh, 2,
                               [](auto x)
                               {
@@ -167,8 +167,10 @@ void interpolate_nedelec(std::shared_ptr<mesh::Mesh<U>> mesh,
       basix::element::dpc_variant::unset, true);
 
   // Create a function space
-  auto V_l = std::make_shared<fem::FunctionSpace<U>>(
-      fem::create_functionspace(mesh, e_l, std::vector<std::size_t>{2}));
+  auto V_l
+      = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
+          mesh, std::make_shared<fem::FiniteElement<U>>(
+                    e_l, std::vector<std::size_t>{2})));
 
   auto u_l = std::make_shared<fem::Function<T>>(V_l);
 

diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp
@@ -36,15 +36,17 @@ int main(int argc, char* argv[])
         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<double>>(
-        fem::create_functionspace(mesh_tet, element_tet,
-                                  std::vector<std::size_t>{3}));
+        fem::create_functionspace<double>(
+            mesh_tet, std::make_shared<fem::FiniteElement<double>>(
+                          element_tet, std::vector<std::size_t>{3})));
 
     basix::FiniteElement element_hex = basix::element::create_lagrange<double>(
         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<double>>(
-        fem::create_functionspace(mesh_hex, element_hex,
-                                  std::vector<std::size_t>{3}));
+        fem::create_functionspace<double>(
+            mesh_hex, std::make_shared<fem::FiniteElement<double>>(
+                          element_hex, std::vector<std::size_t>{3})));
 
     auto u_tet = std::make_shared<fem::Function<T>>(V_tet);
     auto u_hex = std::make_shared<fem::Function<T>>(V_hex);