From 594f190eb1543a359e4b5f58e9b3c0ed1d024d2d Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Fri, 15 Nov 2024 15:31:12 +0000 Subject: [PATCH] Add C++ demo that shows mixed elements and consistency fixes for element 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 --- .github/workflows/ccpp.yml | 2 +- cpp/demo/CMakeLists.txt | 5 +- cpp/demo/biharmonic/main.cpp | 28 +- cpp/demo/codim_0_assembly/main.cpp | 33 +- cpp/demo/custom_kernel/main.cpp | 4 +- cpp/demo/hyperelasticity/main.cpp | 68 ++-- cpp/demo/interpolation-io/main.cpp | 18 +- .../interpolation_different_meshes/main.cpp | 10 +- cpp/demo/mixed_poisson/CMakeLists.txt | 66 ++++ cpp/demo/mixed_poisson/main.cpp | 322 ++++++++++++++++++ cpp/demo/mixed_poisson/poisson.py | 38 +++ cpp/demo/poisson/main.cpp | 31 +- cpp/demo/poisson_matrix_free/main.cpp | 30 +- cpp/doc/source/demo.rst | 1 + cpp/dolfinx/fem/DirichletBC.h | 12 +- cpp/dolfinx/fem/FiniteElement.cpp | 239 +++++++++---- cpp/dolfinx/fem/FiniteElement.h | 221 ++++++++---- cpp/dolfinx/fem/Function.h | 11 +- cpp/dolfinx/fem/FunctionSpace.h | 115 ++----- cpp/dolfinx/fem/assemble_matrix_impl.h | 4 +- cpp/dolfinx/fem/assemble_scalar_impl.h | 1 - cpp/dolfinx/fem/assemble_vector_impl.h | 31 +- cpp/dolfinx/fem/assembler.h | 44 +-- cpp/dolfinx/fem/discreteoperators.h | 18 +- cpp/dolfinx/fem/dofmapbuilder.cpp | 9 +- cpp/dolfinx/fem/dofmapbuilder.h | 4 +- cpp/dolfinx/fem/interpolate.h | 47 +-- cpp/dolfinx/fem/petsc.h | 22 +- cpp/dolfinx/fem/utils.cpp | 13 +- cpp/dolfinx/fem/utils.h | 119 +++---- cpp/dolfinx/io/ADIOS2Writers.h | 3 +- cpp/dolfinx/io/VTKFile.cpp | 6 +- cpp/dolfinx/io/xdmf_function.cpp | 3 +- cpp/dolfinx/mesh/generation.h | 3 +- cpp/dolfinx/mesh/utils.cpp | 4 +- cpp/dolfinx/mesh/utils.h | 17 +- cpp/dolfinx/refinement/interval.h | 4 +- cpp/dolfinx/refinement/plaza.h | 2 +- cpp/test/common/index_map.cpp | 10 +- cpp/test/fem/functionspace.cpp | 3 +- cpp/test/io.cpp | 4 +- cpp/test/matrix.cpp | 6 +- python/demo/demo_stokes.py | 6 +- python/dolfinx/fem/__init__.py | 34 +- python/dolfinx/fem/forms.py | 6 +- python/dolfinx/fem/function.py | 44 +-- python/dolfinx/wrappers/assemble.cpp | 97 ++++-- python/dolfinx/wrappers/fem.cpp | 53 ++- python/dolfinx/wrappers/petsc.cpp | 98 ++++-- .../test_assemble_mesh_independent_form.py | 5 +- python/test/unit/fem/test_assemble_submesh.py | 3 +- .../test/unit/fem/test_mixed_mesh_dofmap.py | 2 +- python/test/unit/fem/test_symmetric_tensor.py | 15 +- 53 files changed, 1312 insertions(+), 682 deletions(-) create mode 100644 cpp/demo/mixed_poisson/CMakeLists.txt create mode 100644 cpp/demo/mixed_poisson/main.cpp create mode 100644 cpp/demo/mixed_poisson/poisson.py diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml index b35ecf84e69..656ae2f362b 100644 --- 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 index 25c5200e6eb..9111ed7fd78 100644 --- 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 index 4dec92293bc..b75d9a51d7f 100644 --- a/cpp/demo/biharmonic/main.cpp +++ b/cpp/demo/biharmonic/main.cpp @@ -134,6 +134,7 @@ #include #include #include +#include #include #include #include @@ -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::create_functionspace(mesh, element)); + auto V + = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(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>(8.0); // Define variational forms - auto a = std::make_shared>(fem::create_form( - *form_biharmonic_a, {V, V}, {}, {{"alpha", alpha}}, {}, {})); - auto L = std::make_shared>( - fem::create_form(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}, {})); + fem::Form a = fem::create_form(*form_biharmonic_a, {V, V}, {}, + {{"alpha", alpha}}, {}, {}); + fem::Form L + = fem::create_form(*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>(0.0, bdofs, V); + fem::DirichletBC 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 u(V); - auto A = la::petsc::Matrix(fem::petsc::create_matrix(*a), false); - la::Vector 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 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(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(b.mutable_array(), {a}, {{bc}}, {}, T(1.0)); b.scatter_rev(std::plus()); - 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 index 018e26756c4..89858c3cf19 100644 --- 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::create_functionspace(mesh, element)); + auto V + = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(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::create_functionspace(submesh, element, {})); + auto W + = std::make_shared>(fem::create_functionspace( + submesh, std::make_shared>(element))); // A mixed-domain form has functions defined over different meshes. // The mesh associated with the measure (dx, ds, etc.) is called the @@ -108,9 +110,8 @@ int main(int argc, char* argv[]) // Next we compute the integration entities on the integration // domain `mesh` std::vector 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>>> @@ -118,23 +119,23 @@ int main(int argc, char* argv[]) = {{fem::IntegralType::cell, {{3, integration_entities}}}}; // We can now create the bilinear form - auto a_mixed = std::make_shared>( - fem::create_form(*form_mixed_codim0_a_mixed, {V, W}, {}, {}, - subdomain_data, entity_maps, V->mesh())); + fem::Form a_mixed + = fem::create_form(*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 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::create_form(*form_mixed_codim0_a, {W, W}, {}, {}, {}, {})); + fem::Form a + = fem::create_form(*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 A(sp); - fem::assemble_matrix(A.mat_add_values(), *a, {}); + fem::assemble_matrix(A.mat_add_values(), a, {}); A.scatter_rev(); std::vector A_mixed_flattened = A_mixed.to_dense(); diff --git a/cpp/demo/custom_kernel/main.cpp b/cpp/demo/custom_kernel/main.cpp index 95f5bb9bedc..d350ef4f68d 100644 --- 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 X(X_b.data(), weights.size(), 2); // Create a scalar function space - auto V = std::make_shared>( - fem::create_functionspace(mesh, e)); + auto V = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(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 index 01790d7408e..ebb7dc91311 100644 --- a/cpp/demo/hyperelasticity/main.cpp +++ b/cpp/demo/hyperelasticity/main.cpp @@ -43,16 +43,15 @@ class HyperElasticProblem { public: /// Constructor - HyperElasticProblem( - std::shared_ptr> L, std::shared_ptr> J, - std::vector>> 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& L, fem::Form& J, + const std::vector>& 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 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(b, *_l); + fem::assemble_vector(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> _l, _j; - std::vector>> _bcs; + fem::Form& _l; + fem::Form& _j; + std::vector>> _bcs; la::Vector _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::create_functionspace(mesh, element, {3})); + auto V + = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>( + element, std::vector{3}))); auto B = std::make_shared>(std::vector{0, 0, 0}); auto traction = std::make_shared>(std::vector{0, 0, 0}); // Define solution function auto u = std::make_shared>(V); - auto a = std::make_shared>( - fem::create_form(*form_hyperelasticity_J_form, {V, V}, {{"u", u}}, - {{"B", B}, {"T", traction}}, {}, {})); - auto L = std::make_shared>( - fem::create_form(*form_hyperelasticity_F_form, {V}, {{"u", u}}, - {{"B", B}, {"T", traction}}, {}, {})); + fem::Form a + = fem::create_form(*form_hyperelasticity_J_form, {V, V}, {{"u", u}}, + {{"B", B}, {"T", traction}}, {}, {}); + fem::Form L + = fem::create_form(*form_hyperelasticity_F_form, {V}, {{"u", u}}, + {{"B", B}, {"T", traction}}, {}, {}); auto u_rotation = std::make_shared>(V); u_rotation->interpolate( @@ -243,10 +245,9 @@ int main(int argc, char* argv[]) } return marker; }); - std::vector bcs = { - std::make_shared>(std::vector{0, 0, 0}, - bdofs_left, V), - std::make_shared>(u_rotation, bdofs_right)}; + std::vector bcs + = {fem::DirichletBC(std::vector{0, 0, 0}, bdofs_left, V), + fem::DirichletBC(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( + *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( family, cell_type, k, basix::element::lagrange_variant::unset, basix::element::dpc_variant::unset, discontinuous); - auto S = std::make_shared>(fem::create_functionspace( - mesh, S_element, std::vector{3, 3})); - auto sigma_expression = fem::create_expression( - *expression_hyperelasticity_sigma, {{"u", u}}, {}); + auto S + = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>( + S_element, std::vector{3, 3}))); - auto sigma = fem::Function(S); + fem::Function 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 index a4904cac9e7..18a54549f75 100644 --- 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, basix::element::dpc_variant::unset, false); // Create a scalar function space - auto V = std::make_shared>( - fem::create_functionspace(mesh, e)); + auto V = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(e))); // Create a finite element Function auto u = std::make_shared>(V); @@ -98,8 +98,8 @@ void interpolate_nedelec(std::shared_ptr> mesh, basix::element::dpc_variant::unset, false); // Create a Nedelec function space - auto V = std::make_shared>( - fem::create_functionspace(mesh, e)); + auto V = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(e))); // Create a Nedelec finite element Function auto u = std::make_shared>(V); @@ -114,7 +114,7 @@ void interpolate_nedelec(std::shared_ptr> 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, 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, basix::element::dpc_variant::unset, true); // Create a function space - auto V_l = std::make_shared>( - fem::create_functionspace(mesh, e_l, std::vector{2})); + auto V_l + = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>( + e_l, std::vector{2}))); auto u_l = std::make_shared>(V_l); diff --git a/cpp/demo/interpolation_different_meshes/main.cpp b/cpp/demo/interpolation_different_meshes/main.cpp index c83cd0852fc..ac7bb577494 100644 --- 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::create_functionspace(mesh_tet, element_tet, - std::vector{3})); + fem::create_functionspace( + mesh_tet, std::make_shared>( + element_tet, std::vector{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::create_functionspace(mesh_hex, element_hex, - std::vector{3})); + fem::create_functionspace( + mesh_hex, std::make_shared>( + element_hex, std::vector{3}))); auto u_tet = std::make_shared>(V_tet); auto u_hex = std::make_shared>(V_hex); diff --git a/cpp/demo/mixed_poisson/CMakeLists.txt b/cpp/demo/mixed_poisson/CMakeLists.txt new file mode 100644 index 00000000000..c9266288a2b --- /dev/null +++ b/cpp/demo/mixed_poisson/CMakeLists.txt @@ -0,0 +1,66 @@ +# This file was generated by running +# +# python cmake/scripts/generate-cmakefiles from dolfinx/cpp +# +cmake_minimum_required(VERSION 3.21) + +set(PROJECT_NAME demo_mixed_poisson) +project(${PROJECT_NAME} LANGUAGES C CXX) + +if(NOT TARGET dolfinx) + find_package(DOLFINX REQUIRED) +endif() + +include(CheckSymbolExists) +set(CMAKE_REQUIRED_INCLUDES ${PETSC_INCLUDE_DIRS}) +check_symbol_exists(PETSC_USE_COMPLEX petscsystypes.h PETSC_SCALAR_COMPLEX) +check_symbol_exists(PETSC_USE_REAL_DOUBLE petscsystypes.h PETSC_REAL_DOUBLE) + +# Add target to compile UFL files +if(PETSC_SCALAR_COMPLEX EQUAL 1) + if(PETSC_REAL_DOUBLE EQUAL 1) + set(SCALAR_TYPE "--scalar_type=complex128") + else() + set(SCALAR_TYPE "--scalar_type=complex64") + endif() +else() + if(PETSC_REAL_DOUBLE EQUAL 1) + set(SCALAR_TYPE "--scalar_type=float64") + else() + set(SCALAR_TYPE "--scalar_type=float32") + endif() +endif() +add_custom_command( + OUTPUT poisson.c + COMMAND ffcx ${CMAKE_CURRENT_SOURCE_DIR}/poisson.py ${SCALAR_TYPE} + VERBATIM + DEPENDS poisson.py + COMMENT "Compile poisson.py using FFCx" +) + +set(CMAKE_INCLUDE_CURRENT_DIR ON) + +add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/poisson.c) +target_link_libraries(${PROJECT_NAME} dolfinx) + +# Set C++20 standard +set(CMAKE_CXX_EXTENSIONS OFF) +target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20) + +# Do not throw error for 'multi-line comments' (these are typical in rst which +# includes LaTeX) +include(CheckCXXCompilerFlag) +check_cxx_compiler_flag("-Wno-comment" HAVE_NO_MULTLINE) +set_source_files_properties( + main.cpp + PROPERTIES + COMPILE_FLAGS + "$<$:-Wno-comment -Wall -Wextra -pedantic -Werror>" +) + +# Test targets (used by DOLFINx testing system) +set(TEST_PARAMETERS2 -np 2 ${MPIEXEC_PARAMS} "./${PROJECT_NAME}") +set(TEST_PARAMETERS3 -np 3 ${MPIEXEC_PARAMS} "./${PROJECT_NAME}") +add_test(NAME ${PROJECT_NAME}_mpi_2 COMMAND "mpirun" ${TEST_PARAMETERS2}) +add_test(NAME ${PROJECT_NAME}_mpi_3 COMMAND "mpirun" ${TEST_PARAMETERS3}) +add_test(NAME ${PROJECT_NAME}_serial COMMAND ${PROJECT_NAME}) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp new file mode 100644 index 00000000000..be4dd3cc334 --- /dev/null +++ b/cpp/demo/mixed_poisson/main.cpp @@ -0,0 +1,322 @@ +// # Mixed Poisson equation +// +// This demo illustrates how to solve Poisson equation using a mixed +// (two-field) formulation. In particular, it illustrates how to +// +// * Create a mixed finite element problem. +// * Extract subspaces. +// * Apply boundary conditions to different fields in a mixed problem. +// * Create integration domain data to execute finite element kernels. +// over subsets of the boundary. +// +// The full implementation is in +// {download}`demo_hyperelasticity/main.cpp`. +// +// +// # Mixed formulation for the Poisson equation +// +// ## Equation and problem definition +// +// A mixed formulation of Poisson equation can be formulated by +// introducing an additional (vector) variable, namely the (negative) +// flux: $\sigma = \nabla u$. The partial differential equations +// then read +// +// $$ +// \begin{align} +// \sigma - \nabla u &= 0 \quad {\rm in} \ \Omega, \\ +// \nabla \cdot \sigma &= - f \quad {\rm in} \ \Omega, +// \end{align} +// $$ +// with boundary conditions +// +// $$ +// u = u_0 \quad {\rm on} \ \Gamma_{D}, \\ +// \sigma \cdot n = g \quad {\rm on} \ \Gamma_{N}. +// $$ +// +// where $n$ denotes the outward unit normal vector on the boundary. We +// see that the boundary condition for the flux ($\sigma \cdot n = g$) +// is an essential boundary condition (which should be enforced in +// the function space), while the other boundary condition ($u = u_0$) +// is a natural boundary condition (which should be applied to the +// variational form). Inserting the boundary conditions, this +// variational problem can be phrased in the general form: find +// $(\sigma, u) \in \Sigma_g \times V$ such that +// +// $$ +// a((\sigma, u), (\tau, v)) = L((\tau, v)) +// \quad \forall \ (\tau, v) \in \Sigma_0 \times V, +// $$ +// +// where the forms $a$ and $L$ are defined as +// +// $$ +// \begin{align} +// a((\sigma, u), (\tau, v)) &:= +// \int_{\Omega} \sigma \cdot \tau + \nabla \cdot \tau \ u +// + \nabla \cdot \sigma \ v \ {\rm d} x, \\ +// L((\tau, v)) &:= - \int_{\Omega} f v \ {\rm d} x +// + \int_{\Gamma_D} u_0 \tau \cdot n \ {\rm d} s, +// \end{align} +// $$ +// and $\Sigma_g := \{ \tau \in H({\rm div})$ such that $\tau \cdot +// n|_{\Gamma_N} = g \}$ and $V := L^2(\Omega)$. +// +// To discretize the above formulation, discrete function spaces +// $\Sigma_h \subset \Sigma$ and $V_h \subset V$ are needed to form a +// mixed function space $\Sigma_h \times V_h$. A stable choice of finite +// element spaces is to let $\Sigma_h$ be a Raviart-Thomas elements of +// polynomial order $k$ and $V_h$ be discontinuous elements of +// polynomial order $k-1$. +// +// We will use the same definitions of functions and boundaries as in the +// demo for {doc}`the Poisson equation `. These are: +// +// * $\Omega = [0,1] \times [0,1]$ (a unit square) +// * $\Gamma_{D} = \{(0, y) \cup (1, y) \in \partial \Omega\}$ +// * $\Gamma_{N} = \{(x, 0) \cup (x, 1) \in \partial \Omega\}$ +// * $u_0 = 20 y + 1$ on $\Gamma_{D}$ +// * $g = 10$ (flux) on $\Gamma_{N}$ +// * $f = \sin(5x - 0.5) + 1 (source term) + +// ## UFL form file +// +// The UFL file is implemented in +// {download}`demo_mixed_poisson/poisson.py`. +// ````{admonition} UFL form implemented in python +// :class: dropdown +// ![ufl-code] +// ```` +// + +#include "poisson.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dolfinx; +using T = PetscScalar; +using U = typename dolfinx::scalar_value_type_t; + +int main(int argc, char* argv[]) +{ + dolfinx::init_logging(argc, argv); + PetscInitialize(&argc, &argv, nullptr, nullptr); + + { + // Create mesh + auto mesh = std::make_shared>( + mesh::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, + {32, 32}, mesh::CellType::triangle)); + + // Create Basix elements + auto RT = basix::create_element( + basix::element::family::RT, basix::cell::type::triangle, 1, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + auto P0 = basix::create_element( + basix::element::family::P, basix::cell::type::triangle, 0, + basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, true); + + // Create DOLFINx mixed element + auto ME = std::make_shared>( + std::vector>{{RT}, {P0}}); + + // Create FunctionSpace + auto V = std::make_shared>( + fem::create_functionspace(mesh, ME)); + + // Get subspaces (views into V) + auto V0 = std::make_shared>(V->sub({0})); + auto V1 = std::make_shared>(V->sub({1})); + + // Collapse spaces + auto W0 = std::make_shared>(V0->collapse().first); + auto W1 = std::make_shared>(V1->collapse().first); + + // Create source function and interpolate sin(5x) + 1 + auto f = std::make_shared>(W1); + f->interpolate( + [](auto x) -> std::pair, std::vector> + { + std::vector f; + for (std::size_t p = 0; p < x.extent(1); ++p) + { + auto x0 = x(0, p); + f.push_back(std::sin(5 * x0) + 1); + } + return {f, {f.size()}}; + }); + + // Boundary condition value for u and interpolate 20y + 1 + auto u0 = std::make_shared>(W1); + u0->interpolate( + [](auto x) -> std::pair, std::vector> + { + std::vector f; + for (std::size_t p = 0; p < x.extent(1); ++p) + f.push_back(20 * x(1, p) + 1); + return {f, {f.size()}}; + }); + + // Create boundary condition for \sigma and interpolate such that + // flux = 10 (for top and bottom boundaries) + auto g = std::make_shared>(W0); + g->interpolate( + [](auto x) -> std::pair, std::vector> + { + using mspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 2, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>; + + std::vector fdata(2 * x.extent(1), 0); + mspan_t f(fdata.data(), 2, x.extent(1)); + for (std::size_t p = 0; p < x.extent(1); ++p) + f(1, p) = x(1, p) < 0.5 ? -10 : 10; + return {std::move(fdata), {2, x.extent(1)}}; + }); + + // Get list of all boundary facets + mesh->topology()->create_connectivity(1, 2); + std::vector bfacets = mesh::exterior_facet_indices(*mesh->topology()); + + // Get facets with boundary condition on u + std::vector dfacets = mesh::locate_entities_boundary( + *mesh, 1, + [](auto x) + { + using U = typename decltype(x)::value_type; + constexpr U eps = 1.0e-8; + std::vector marker(x.extent(1), false); + for (std::size_t p = 0; p < x.extent(1); ++p) + { + auto x0 = x(0, p); + if (std::abs(x0) < eps or std::abs(x0 - 1) < eps) + marker[p] = true; + } + return marker; + }); + + // Compute facets with \sigma (flux) boundary condition facets, which is + // {all boundary facet} - {u0 boundary facets } + std::vector nfacets; + std::ranges::set_difference(bfacets, dfacets, std::back_inserter(nfacets)); + + // Get dofs that are constrained by \sigma + std::array, 2> ndofs + = fem::locate_dofs_topological( + *mesh->topology(), {*V0->dofmap(), *W0->dofmap()}, 1, nfacets); + + // Create boundary condition for \sigma. \sigma \cdot n will be + // constrained to to be equal to the normal component of g. The + // boundary conditions are applied to degrees-of-freedom ndofs, and + // V0 is the subspace that is constrained. + fem::DirichletBC bc(g, ndofs, V0); + + // Create integration domain data for u0 boundary condition (applied + // on the ds(1) in the UFL file). First we get facet data + // integration data for facets in dfacets. + std::vector domains = fem::compute_integration_domains( + fem::IntegralType::exterior_facet, *mesh->topology(), dfacets); + + // Create data structure for the ds(1) integration domain in form + // (see the UFL file). It is for en exterior facet integral (the key + // in the map), and exterior facet domain marked as '1' in the UFL + // file, and 'domains' holds the necessary data to perform + // integration of selected facets. + std::map< + fem::IntegralType, + std::vector>>> + subdomain_data{{fem::IntegralType::exterior_facet, {{1, domains}}}}; + + // Define variational forms and attach he required data + fem::Form a = fem::create_form(*form_poisson_a, {V, V}, {}, {}, + subdomain_data, {}); + fem::Form L = fem::create_form( + *form_poisson_L, {V}, {{"f", f}, {"u0", u0}}, {}, subdomain_data, {}); + + // Create solution finite element Function + auto u = std::make_shared>(V); + + // Create matrix and RHS vector data structures + auto A = la::petsc::Matrix(fem::petsc::create_matrix(a), false); + la::Vector b(L.function_spaces()[0]->dofmap()->index_map, + L.function_spaces()[0]->dofmap()->index_map_bs()); + + // Assemble the bilinear form into a matrix. The PETSc matrix is + // 'flushed' so we can set values in it in the subsequent step. + MatZeroEntries(A.mat()); + fem::assemble_matrix(la::petsc::Matrix::set_fn(A.mat(), ADD_VALUES), a, + {bc}); + MatAssemblyBegin(A.mat(), MAT_FLUSH_ASSEMBLY); + MatAssemblyEnd(A.mat(), MAT_FLUSH_ASSEMBLY); + + // Set '1' on diagonal for Dirichlet dofs + fem::set_diagonal(la::petsc::Matrix::set_fn(A.mat(), INSERT_VALUES), *V, + {bc}); + MatAssemblyBegin(A.mat(), MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY); + + // Assemble the linear form L into RHS vector + b.set(0); + fem::assemble_vector(b.mutable_array(), L); + + // Modify unconstrained dofs on RHS to account for Dirichlet BC dofs + // (constrained dofs), and perform parallel update on the vector. + fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, T(1)); + b.scatter_rev(std::plus()); + + // Set value for constrained dofs + bc.set(b.mutable_array(), std::nullopt); + + // Create PETSc linear solver + la::petsc::KrylovSolver lu(MPI_COMM_WORLD); + la::petsc::options::set("ksp_type", "preonly"); + la::petsc::options::set("pc_type", "lu"); + if (sizeof(PetscInt) == 4) + la::petsc::options::set("pc_factor_mat_solver_type", "mumps"); + else + la::petsc::options::set("pc_factor_mat_solver_type", "superlu_dist"); + lu.set_from_options(); + + // Solve linear system Ax = b + lu.set_operator(A.mat()); + la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false); + la::petsc::Vector _b(la::petsc::create_vector_wrap(b), false); + lu.solve(_u.vec(), _b.vec()); + + // Update ghost values before output + u->x()->scatter_fwd(); + + // Save solution in VTK format + auto u_soln = std::make_shared>(u->sub(1).collapse()); + io::VTKFile file(MPI_COMM_WORLD, "u.pvd", "w"); + file.write({*u_soln}, 0); + +#ifdef HAS_ADIOS2 + // Save solution in VTX format + io::VTXWriter vtx(MPI_COMM_WORLD, "u.bp", {u_soln}, "bp4"); + vtx.write(0); +#endif + } + + PetscFinalize(); + + return 0; +} diff --git a/cpp/demo/mixed_poisson/poisson.py b/cpp/demo/mixed_poisson/poisson.py new file mode 100644 index 00000000000..69e93ec6fcf --- /dev/null +++ b/cpp/demo/mixed_poisson/poisson.py @@ -0,0 +1,38 @@ +# The first step is to define the variational problem at hand. We define +# the variational problem in UFL terms in a separate form file +# {download}`mixed-poisson/poisson.py`. We begin by defining the finite +# element: + +from basix.ufl import element, mixed_element +from ufl import ( + Coefficient, + FacetNormal, + FunctionSpace, + Measure, + Mesh, + TestFunctions, + TrialFunctions, + div, + inner, +) + +shape = "triangle" +RT = element("RT", shape, 1) +P = element("DP", shape, 0) +ME = mixed_element([RT, P]) + +msh = Mesh(element("Lagrange", shape, 1, shape=(2,))) +n = FacetNormal(msh) +V = FunctionSpace(msh, ME) + +(sigma, u) = TrialFunctions(V) +(tau, v) = TestFunctions(V) + +V0 = FunctionSpace(msh, P) +f = Coefficient(V0) +u0 = Coefficient(V0) + +dx = Measure("dx", msh) +ds = Measure("ds", msh) +a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx +L = -inner(f, v) * dx + inner(u0 * n, tau) * ds(1) diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp index 6d5d6624d4e..d14d73153d7 100644 --- a/cpp/demo/poisson/main.cpp +++ b/cpp/demo/poisson/main.cpp @@ -124,8 +124,9 @@ int main(int argc, char* argv[]) basix::element::lagrange_variant::unset, basix::element::dpc_variant::unset, false); - auto V = std::make_shared>( - fem::create_functionspace(mesh, element, {})); + auto V + = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(element))); // Next, we define the variational formulation by initializing the // bilinear and linear forms ($a$, $L$) using the previously @@ -139,10 +140,10 @@ int main(int argc, char* argv[]) auto g = std::make_shared>(V); // Define variational forms - auto a = std::make_shared>(fem::create_form( - *form_poisson_a, {V, V}, {}, {{"kappa", kappa}}, {}, {})); - auto L = std::make_shared>(fem::create_form( - *form_poisson_L, {V}, {{"f", f}, {"g", g}}, {}, {}, {})); + fem::Form a = fem::create_form(*form_poisson_a, {V, V}, {}, + {{"kappa", kappa}}, {}, {}); + fem::Form L = fem::create_form(*form_poisson_L, {V}, + {{"f", f}, {"g", g}}, {}, {}, {}); // Now, the Dirichlet boundary condition ($u = 0$) can be created // using the class {cpp:class}`DirichletBC`. A @@ -157,7 +158,7 @@ int main(int argc, char* argv[]) // Define boundary condition - auto facets = mesh::locate_entities_boundary( + std::vector facets = mesh::locate_entities_boundary( *mesh, 1, [](auto x) { @@ -172,9 +173,9 @@ int main(int argc, char* argv[]) } return marker; }); - const auto bdofs = fem::locate_dofs_topological( + std::vector bdofs = fem::locate_dofs_topological( *V->mesh()->topology_mutable(), *V->dofmap(), 1, facets); - auto bc = std::make_shared>(0.0, bdofs, V); + fem::DirichletBC bc(0.0, bdofs, V); f->interpolate( [](auto x) -> std::pair, std::vector> @@ -207,13 +208,13 @@ int main(int argc, char* argv[]) // and `bc` as follows: auto u = std::make_shared>(V); - auto A = la::petsc::Matrix(fem::petsc::create_matrix(*a), false); - la::Vector b(L->function_spaces()[0]->dofmap()->index_map, - L->function_spaces()[0]->dofmap()->index_map_bs()); + la::petsc::Matrix A(fem::petsc::create_matrix(a), false); + la::Vector 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(la::petsc::Matrix::set_fn(A.mat(), INSERT_VALUES), *V, @@ -222,10 +223,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(b.mutable_array(), {a}, {{bc}}, {}, T(1)); b.scatter_rev(std::plus()); - 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/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp index df7bb19aa41..270ef7223da 100644 --- a/cpp/demo/poisson_matrix_free/main.cpp +++ b/cpp/demo/poisson_matrix_free/main.cpp @@ -143,20 +143,20 @@ void solver(MPI_Comm comm) basix::element::family::P, basix::cell::type::triangle, 2, basix::element::lagrange_variant::unset, basix::element::dpc_variant::unset, false); - auto V = std::make_shared>( - fem::create_functionspace(mesh, element, {})); + auto V = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(element))); // Prepare and set Constants for the bilinear form auto f = std::make_shared>(-6.0); // Define variational forms - auto L = std::make_shared>( - fem::create_form(*form_poisson_L, {V}, {}, {{"f", f}}, {}, {})); + fem::Form L + = fem::create_form(*form_poisson_L, {V}, {}, {{"f", f}}, {}, {}); // Action of the bilinear form "a" on a function ui auto ui = std::make_shared>(V); - auto M = std::make_shared>( - fem::create_form(*form_poisson_M, {V}, {{"ui", ui}}, {{}}, {}, {})); + fem::Form M + = fem::create_form(*form_poisson_M, {V}, {{"ui", ui}}, {{}}, {}, {}); // Define boundary condition auto u_D = std::make_shared>(V); @@ -178,12 +178,12 @@ void solver(MPI_Comm comm) // Assemble RHS vector la::Vector b(V->dofmap()->index_map, V->dofmap()->index_map_bs()); - fem::assemble_vector(b.mutable_array(), *L); + fem::assemble_vector(b.mutable_array(), L); // Apply lifting to account for Dirichlet boundary condition // b <- b - A * x_bc bc->set(ui->x()->mutable_array(), std::nullopt, T(-1)); - fem::assemble_vector(b.mutable_array(), *M); + fem::assemble_vector(b.mutable_array(), M); // Communicate ghost values b.scatter_rev(std::plus()); @@ -194,8 +194,8 @@ void solver(MPI_Comm comm) b.scatter_fwd(); // Pack coefficients and constants - auto coeff = fem::allocate_coefficient_storage(*M); - std::vector constants = fem::pack_constants(*M); + auto coeff = fem::allocate_coefficient_storage(M); + std::vector constants = fem::pack_constants(M); // Create function for computing the action of A on x (y = Ax) auto action = [&M, &ui, &bc, &coeff, &constants](auto& x, auto& y) @@ -207,8 +207,8 @@ void solver(MPI_Comm comm) std::ranges::copy(x.array(), ui->x()->mutable_array().begin()); // Compute action of A on x - fem::pack_coefficients(*M, coeff); - fem::assemble_vector(y.mutable_array(), *M, std::span(constants), + fem::pack_coefficients(M, coeff); + fem::assemble_vector(y.mutable_array(), M, std::span(constants), fem::make_coefficients_span(coeff)); // Set BC dofs to zero (effectively zeroes rows of A) @@ -230,9 +230,9 @@ void solver(MPI_Comm comm) // Compute L2 error (squared) of the solution vector e = (u - u_d, u // - u_d)*dx - auto E = std::make_shared>(fem::create_form( - *form_poisson_E, {}, {{"uexact", u_D}, {"usol", u}}, {}, {}, {}, mesh)); - T error = fem::assemble_scalar(*E); + fem::Form E = fem::create_form( + *form_poisson_E, {}, {{"uexact", u_D}, {"usol", u}}, {}, {}, {}, mesh); + T error = fem::assemble_scalar(E); if (dolfinx::MPI::rank(comm) == 0) { std::cout << "Number of CG iterations " << num_it << std::endl; diff --git a/cpp/doc/source/demo.rst b/cpp/doc/source/demo.rst index 7ce9d8d80db..8b6c8a35349 100644 --- a/cpp/doc/source/demo.rst +++ b/cpp/doc/source/demo.rst @@ -23,6 +23,7 @@ Intermediate :maxdepth: 1 demos/demo_hyperelasticity.md + demos/demo_mixed_poisson.md demos/demo_poisson_matrix_free.md demos/demo_interpolation_different_meshes.md demos/demo_interpolation-io.md diff --git a/cpp/dolfinx/fem/DirichletBC.h b/cpp/dolfinx/fem/DirichletBC.h index b43066ff253..e0aa9f29c4a 100644 --- a/cpp/dolfinx/fem/DirichletBC.h +++ b/cpp/dolfinx/fem/DirichletBC.h @@ -336,7 +336,7 @@ class DirichletBC { assert(g); assert(V); - if (g->shape.size() != V->value_shape().size()) + if (g->shape.size() != V->element()->value_shape().size()) { throw std::runtime_error( "Rank mis-match between Constant and function space in DirichletBC"); @@ -419,8 +419,10 @@ class DirichletBC DirichletBC(std::shared_ptr> g, X&& V_g_dofs, std::shared_ptr> V) : _function_space(V), _g(g), - _dofs0(std::forward(V_g_dofs[0])), - _dofs1_g(std::forward(V_g_dofs[1])), + _dofs0(std::forward::value_type>( + V_g_dofs[0])), + _dofs1_g(std::forward::value_type>( + V_g_dofs[1])), _owned_indices0(num_owned(*_function_space->dofmap(), _dofs0)) { } @@ -513,7 +515,7 @@ class DirichletBC auto dofs1_g = _dofs1_g.empty() ? std::span(_dofs0) : std::span(_dofs1_g); std::span values = g->x()->array(); - if (x0.has_value()) + if (x0) { std::span _x0 = x0.value(); assert(x.size() <= _x0.size()); @@ -543,7 +545,7 @@ class DirichletBC auto g = std::get>>(_g); const std::vector& value = g->value; std::int32_t bs = _function_space->dofmap()->bs(); - if (x0.has_value()) + if (x0) { assert(x.size() <= x0.value().size()); std::ranges::for_each( diff --git a/cpp/dolfinx/fem/FiniteElement.cpp b/cpp/dolfinx/fem/FiniteElement.cpp index 28b7f95f8d3..464d288ea63 100644 --- a/cpp/dolfinx/fem/FiniteElement.cpp +++ b/cpp/dolfinx/fem/FiniteElement.cpp @@ -21,6 +21,26 @@ using namespace dolfinx::fem; namespace { +/// @brief Create a list of fem::FiniteElement from Basix elements and +/// other data. +/// @tparam T +/// @param elements +/// @return List of DOLFINx elements +template +std::vector>> +_build_element_list(std::vector> elements) +{ + std::vector>> _e; + std::ranges::transform(elements, std::back_inserter(_e), + [](auto data) + { + auto& [e, bs, symm] = data; + return std::make_shared>(e, bs, + symm); + }); + return _e; +} + /// Recursively extract sub finite element template std::shared_ptr> @@ -61,62 +81,64 @@ _extract_sub_element(const FiniteElement& finite_element, return _extract_sub_element(*sub_element, sub_component); } -} // namespace -//----------------------------------------------------------------------------- -template -FiniteElement::FiniteElement(mesh::CellType cell_type, - std::span points, - std::array pshape, - std::size_t block_size, bool symmetric) - : _signature("Quadrature element " + std::to_string(pshape[0]) + " " - + std::to_string(block_size)), - _space_dim(pshape[0] * block_size), - _reference_value_shape(std::vector{}), _bs(block_size), - _symmetric(symmetric), _needs_dof_permutations(false), - _needs_dof_transformations(false), - _entity_dofs(mesh::cell_dim(cell_type) + 1), - _entity_closure_dofs(mesh::cell_dim(cell_type) + 1), - _points(std::vector(points.begin(), points.end()), pshape) -{ - const int tdim = mesh::cell_dim(cell_type); - for (int d = 0; d <= tdim; ++d) - { - int num_entities = mesh::cell_num_entities(cell_type, d); - _entity_dofs[d].resize(num_entities); - _entity_closure_dofs[d].resize(num_entities); - } +} // namespace - for (std::size_t i = 0; i < pshape[0]; ++i) - { - _entity_dofs[tdim][0].push_back(i); - _entity_closure_dofs[tdim][0].push_back(i); - } -} //----------------------------------------------------------------------------- template -FiniteElement::FiniteElement(const basix::FiniteElement& element, - std::size_t block_size, bool symmetric) - : _space_dim(block_size * element.dim()), - _reference_value_shape(element.value_shape()), _bs(block_size), +FiniteElement::FiniteElement( + const basix::FiniteElement& element, + std::optional> value_shape, bool symmetric) + : _value_shape(value_shape ? *value_shape : element.value_shape()), + _bs(value_shape + ? std::accumulate(value_shape->begin(), value_shape->end(), 1, + std::multiplies{}) + : 1), + _cell_type(mesh::cell_type_from_basix_type(element.cell_type())), + _space_dim(_bs * element.dim()), + _sub_elements( + value_shape + ? std::vector< + std::shared_ptr>>( + _bs, std::make_shared>(element)) + : std::vector< + std::shared_ptr>>()), + _reference_value_shape(element.value_shape()), _element(std::make_unique>(element)), _symmetric(symmetric), _needs_dof_permutations( - !_element->dof_transformations_are_identity() - and _element->dof_transformations_are_permutations()), + !element.dof_transformations_are_identity() + and element.dof_transformations_are_permutations()), _needs_dof_transformations( - !_element->dof_transformations_are_identity() - and !_element->dof_transformations_are_permutations()), + !element.dof_transformations_are_identity() + and !element.dof_transformations_are_permutations()), _entity_dofs(element.entity_dofs()), _entity_closure_dofs(element.entity_closure_dofs()) { - // Create all sub-elements - if (_bs > 1) + // TODO: symmetric rank-2 symmetric tensors are presently constructed + // as rank-1 tensors, e.g. a rank-2 symmetric tensor in 3D is + // constructed as rank-1 with shape (6,). It should be really be + // shape=(3, 3) with block size 6. + + // If element is blocked, check that base element is scalar + if (value_shape and !element.value_shape().empty()) { - _sub_elements - = std::vector>>( - _bs, std::make_shared>(element, 1)); - _reference_value_shape = {block_size}; + throw std::runtime_error("Blocked finite elements can be constructed only " + "from scalar base elements."); + } + + if (symmetric) // Consistency check for symmetric elements + { + if (!value_shape) + { + throw std::runtime_error( + "Symmetric elements required value shape to be supplied."); + } + // else if (block_shape->size() + // != 2) // See below TODO on symmetric rank-2 tensors + // { + // throw std::runtime_error("Symmetric elements must be rank-2."); + // } } std::string family; @@ -137,11 +159,19 @@ FiniteElement::FiniteElement(const basix::FiniteElement& element, } //----------------------------------------------------------------------------- template +FiniteElement::FiniteElement(std::vector> elements) + : FiniteElement(_build_element_list(elements)) +{ +} +//----------------------------------------------------------------------------- +template FiniteElement::FiniteElement( const std::vector>>& elements) - : _space_dim(0), _sub_elements(elements), - _reference_value_shape(std::nullopt), _bs(1), _symmetric(false), - _needs_dof_permutations(false), _needs_dof_transformations(false) + : _value_shape(std::nullopt), _bs(1), + _cell_type(elements.front()->cell_type()), _space_dim(-1), + _sub_elements(elements), _reference_value_shape(std::nullopt), + _symmetric(false), _needs_dof_permutations(false), + _needs_dof_transformations(false) { if (elements.size() < 2) { @@ -152,8 +182,9 @@ FiniteElement::FiniteElement( _signature = "Mixed element ("; const std::vector>>& ed - = elements[0]->entity_dofs(); + = elements.front()->entity_dofs(); _entity_dofs.resize(ed.size()); + _entity_closure_dofs.resize(ed.size()); for (std::size_t i = 0; i < ed.size(); ++i) { @@ -176,8 +207,8 @@ FiniteElement::FiniteElement( { for (std::size_t j = 0; j < _entity_dofs[i].size(); ++j) { - const std::vector sub_ed = e->entity_dofs()[i][j]; - const std::vector sub_ecd = e->entity_closure_dofs()[i][j]; + std::vector sub_ed = e->entity_dofs()[i][j]; + std::vector sub_ecd = e->entity_closure_dofs()[i][j]; for (auto k : sub_ed) { for (std::size_t b = 0; b < sub_bs; ++b) @@ -199,6 +230,46 @@ FiniteElement::FiniteElement( } //----------------------------------------------------------------------------- template +FiniteElement::FiniteElement(mesh::CellType cell_type, + std::span points, + std::array pshape, + std::vector value_shape, + bool symmetric) + : _value_shape(value_shape), + _bs(std::accumulate(value_shape.begin(), value_shape.end(), 1, + std::multiplies{})), + _cell_type(cell_type), + _signature("Quadrature element " + std::to_string(pshape[0])), + _space_dim(_bs * pshape[0]), _sub_elements({}), + _reference_value_shape(std::vector()), _element(nullptr), + _symmetric(symmetric), _needs_dof_permutations(false), + _needs_dof_transformations(false), + _entity_dofs(mesh::cell_dim(cell_type) + 1), + _entity_closure_dofs(mesh::cell_dim(cell_type) + 1), + _points(std::vector(points.begin(), points.end()), pshape) +{ + assert(value_shape.size() <= 1); + _bs = std::accumulate(value_shape.begin(), value_shape.end(), 1, + std::multiplies{}); + _space_dim = pshape[0] * _bs; + _signature += " " + std::to_string(_bs); + + const int tdim = mesh::cell_dim(cell_type); + for (int d = 0; d <= tdim; ++d) + { + int num_entities = mesh::cell_num_entities(cell_type, d); + _entity_dofs[d].resize(num_entities); + _entity_closure_dofs[d].resize(num_entities); + } + + for (std::size_t i = 0; i < pshape[0]; ++i) + { + _entity_dofs[tdim][0].push_back(i); + _entity_closure_dofs[tdim][0].push_back(i); + } +} +//----------------------------------------------------------------------------- +template bool FiniteElement::operator==(const FiniteElement& e) const { if (!_element or !e._element) @@ -217,6 +288,12 @@ bool FiniteElement::operator!=(const FiniteElement& e) const } //----------------------------------------------------------------------------- template +mesh::CellType FiniteElement::cell_type() const noexcept +{ + return _cell_type; +} +//----------------------------------------------------------------------------- +template std::string FiniteElement::signature() const noexcept { return _signature; @@ -229,6 +306,56 @@ int FiniteElement::space_dimension() const noexcept } //----------------------------------------------------------------------------- template +int FiniteElement::value_size() const +{ + if (_value_shape) + { + int vs = std::accumulate(_value_shape->begin(), _value_shape->end(), 1, + std::multiplies{}); + + // See comments in constructor on why special handling for the + // symmetric case is required. + if (_symmetric) + { + if (vs == 3) + return 4; + else if (vs == 6) + return 9; + else + { + throw std::runtime_error( + "Inconsistent size for symmetric rank-2 tensor."); + } + } + else + return vs; + } + else + throw std::runtime_error("Element does not have a value_shape."); +} +//----------------------------------------------------------------------------- +template +std::span FiniteElement::value_shape() const +{ + if (_value_shape) + return *_value_shape; + else + throw std::runtime_error("Element does not have a value_shape."); +} +//----------------------------------------------------------------------------- +template +int FiniteElement::reference_value_size() const +{ + if (_reference_value_shape) + { + return std::accumulate(_reference_value_shape->begin(), + _reference_value_shape->end(), 1, std::multiplies{}); + } + else + throw std::runtime_error("Element does not have a reference_value_shape."); +} +//----------------------------------------------------------------------------- +template std::span FiniteElement::reference_value_shape() const { if (_reference_value_shape) @@ -258,18 +385,6 @@ bool FiniteElement::symmetric() const } //----------------------------------------------------------------------------- template -int FiniteElement::reference_value_size() const -{ - if (_reference_value_shape) - { - return std::accumulate(_reference_value_shape->begin(), - _reference_value_shape->end(), 1, std::multiplies{}); - } - else - throw std::runtime_error("Element does not have a reference_value_shape."); -} -//----------------------------------------------------------------------------- -template int FiniteElement::block_size() const noexcept { return _bs; diff --git a/cpp/dolfinx/fem/FiniteElement.h b/cpp/dolfinx/fem/FiniteElement.h index c30cf18e9d1..b16056a0b50 100644 --- a/cpp/dolfinx/fem/FiniteElement.h +++ b/cpp/dolfinx/fem/FiniteElement.h @@ -1,4 +1,4 @@ -// Copyright (C) 2020-2021 Garth N. Wells and Matthew W. Scroggs +// Copyright (C) 2020-2024 Garth N. Wells and Matthew W. Scroggs // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -30,6 +30,24 @@ enum class doftransform inverse_transpose = 3, ///< Transpose inverse }; +/// @brief Basix element holder +/// @tparam T Scalar type +template +struct BasixElementData +{ + std::reference_wrapper> + element; ///< Finite element. + std::optional> value_shape + = std::nullopt; ///< Value shape. Can only be set for scalar `element`. + bool symmetry = false; ///< Symmetry. Should ony set set for 2nd-order tensor + ///< blocked elements. +}; + +/// Type deduction helper +template +BasixElementData(U element, V bs, W symmetry) + -> BasixElementData::type::scalar_type>; + /// @brief Model of a finite element. /// /// Provides the dof layout on a reference element, and various methods @@ -42,26 +60,39 @@ class FiniteElement using geometry_type = T; /// @brief Create a finite element from a Basix finite element. - /// @param[in] element Basix finite element - /// @param[in] block_size The block size for the element - /// @param[in] symmetric Is the element a symmetric tensor? + /// @param[in] element Basix finite element. + /// @param[in] value_shape Value shape for blocked element, e.g. `{3}` + /// for a vector in 3D or `{2, 2}` for a rank-2 tensor in 2D. Can only + /// be set for blocked scalar `element`. For other elements and scalar + /// elements it should be `std::nullopt`. + /// @param[in] symmetric Is the element a symmetric tensor? Should ony + /// set for 2nd-order tensor blocked elements. FiniteElement(const basix::FiniteElement& element, - std::size_t block_size, bool symmetric = false); + std::optional> value_shape + = std::nullopt, + bool symmetric = false); + + /// @brief Create a mixed finite element from Basix finite elements. + /// @param[in] elements List of (Basix finite element, block size, + /// symmetric) tuples, one for each element in the mixed element. + FiniteElement(std::vector> elements); - /// @brief Create mixed finite element from a list of finite elements. - /// @param[in] elements Basix finite elements + /// @brief Create a mixed finite element from a list of finite + /// elements. + /// @param[in] elements Finite elements to compose the mixed element. FiniteElement( const std::vector>>& elements); - /// @brief Create a quadrature element - /// @param[in] cell_type The cell type - /// @param[in] points Quadrature points - /// @param[in] pshape Shape of points array - /// @param[in] block_size The block size for the element + /// @brief Create a quadrature element. + /// @param[in] cell_type Cell type. + /// @param[in] points Quadrature points. + /// @param[in] pshape Shape of `points` array. + /// @param[in] value_shape Value shape for the element. /// @param[in] symmetric Is the element a symmetric tensor? FiniteElement(mesh::CellType cell_type, std::span points, - std::array pshape, std::size_t block_size, + std::array pshape, + std::vector value_shape = {}, bool symmetric = false); /// Copy constructor @@ -79,19 +110,22 @@ class FiniteElement /// Move assignment FiniteElement& operator=(FiniteElement&& element) = default; - /// Check if two elements are equivalent - /// @return True is the two elements are the same + /// @brief Check if two elements are equivalent. + /// @return True is the two elements are the same. /// @note Equality can be checked only for non-mixed elements. For a - /// mixed element, this function will raise an exception. + /// mixed element, this function will throw an exception. bool operator==(const FiniteElement& e) const; - /// Check if two elements are not equivalent - /// @return True is the two elements are not the same + /// @brief Check if two elements are not equivalent. + /// @return True is the two elements are not the same. /// @note Equality can be checked only for non-mixed elements. For a /// mixed element, this function will raise an exception. bool operator!=(const FiniteElement& e) const; - /// String identifying the finite element + /// @brief Cell shape that the element is defined on. + mesh::CellType cell_type() const noexcept; + + /// @brief String identifying the finite element. /// @return Element signature /// @warning The function is provided for convenience, but it should /// not be relied upon for determining the element type. Use other @@ -99,43 +133,95 @@ class FiniteElement /// properties. std::string signature() const noexcept; - /// Dimension of the finite element function space (the number of - /// degrees-of-freedom for the element) - /// @return Dimension of the finite element space + /// @brief Dimension of the finite element function space (the number + /// of degrees-of-freedom for the element). + /// + /// For 'blocked' elements, this function returns the dimension of the + /// full element rather than the dimension of the base element. + /// + /// @return Dimension of the finite element space. int space_dimension() const noexcept; /// @brief Block size of the finite element function space. /// - /// For BlockedElements, this is the number of DOFs colocated at each - /// DOF point. For other elements, this is always 1. - /// @return Block size of the finite element space + /// For non-blocked elements, this is always 1. For blocked elements, + /// this is the number of DOFs collocated at each DOF point. For + /// blocked elements the block size is equal to the value size, except + /// for symmetric rank-2 tensor blocked elements. For a symmetric + /// rank-2 tensor blocked element the block size is 3 in 2D and 6 in + /// 3D. + /// + /// @return Block size of the finite element space. int block_size() const noexcept; - /// @brief Value size. + /// @brief Value size of the finite element field. + /// + /// The value size is the number of components in the finite element + /// field. It is the product of the value shape, e.g. is is 1 for a + /// scalar function, 2 for a 2D vector, 9 for a second-order tensor in + /// 3D, etc. For blocked elements, this function returns the value + /// size for the full 'blocked' element. + /// + /// @note The return value of this function is inconsistent with + /// value_shape() for rank-2 'symmetric' elements. Due to issues + /// elsewhere in the code base, rank-2 symmetric fields have value + /// shape `{3}` (2D) or `{6}` rather than `{2, 2}` and `{3, 3}`, + /// respectively. For symmetric rank-2 tensors this function returns 4 + /// for 2D cases and 9 for 3D cases. This inconsistency will be fixed + /// in the future. + /// + /// @throws Exception is thrown for a mixed element as mixed elements + /// do not have a value shape. + /// @return The value size. + int value_size() const; + + /// @brief Value shape of the finite element field. + /// + /// The value shape describes the shape of the finite element field, + /// e.g. `{}` for a scalar, `{2}` for a vector in 2D, `{3, 3}` for a + /// rank-2 tensor in 3D, etc. + /// + /// @throws Exception is thrown for a mixed element as mixed elements + /// do not have a value shape. + /// @return The value shape. + std::span value_shape() const; + + /// @brief Value size of the base (non-blocked) finite element field. + /// + /// The reference value size is the product of the reference value + /// shape, e.g. it is 1 for a scalar element, 2 for a 2D + /// (non-blocked) vector, 9 for a (non-blocked) second-order tensor in + /// 3D, etc. + /// + /// For blocked elements, this function returns the value shape for + /// the 'base' element from which the blocked element is composed. For + /// other elements, the return value is the same as + /// FiniteElement::value_shape. /// - /// The value size is the product of the value shape, e.g. is is 1 - /// for a scalar function, 2 for a 2D vector, 9 for a second-order - /// tensor in 3D. /// @throws Exception is thrown for a mixed element as mixed elements /// do not have a value shape. /// @return The value size. int reference_value_size() const; - /// @brief Value shape. + /// @brief Value shape of the base (non-blocked) finite element field. + /// + /// For non-blocked elements, this function returns the same as + /// FiniteElement::value_shape. For blocked and quadrature elements + /// the returned shape will be `{}`. + /// + /// Mixed elements do not have a reference value shape. /// - /// The value shape described the shape of the finite element field, - /// e.g. {} for a scalar, {3, 3} for a tensor in 3D. Mixed elements do - /// not have a value shape. /// @throws Exception is thrown for a mixed element as mixed elements /// do not have a value shape. /// @return The value shape. std::span reference_value_shape() const; - /// The local DOFs associated with each subentity of the cell + /// @brief Local DOFs associated with each sub-entity of the cell. const std::vector>>& entity_dofs() const noexcept; - /// The local DOFs associated with the closure of each subentity of the cell + /// @brief Local DOFs associated with the closure of each sub-entity + /// of the cell. const std::vector>>& entity_closure_dofs() const noexcept; @@ -144,41 +230,43 @@ class FiniteElement /// @brief Evaluate derivatives of the basis functions up to given order /// at points in the reference cell. + /// /// @param[in,out] values Array that will be filled with the tabulated /// basis values. Must have shape `(num_derivatives, num_points, /// num_dofs, reference_value_size)` (row-major storage) /// @param[in] X The reference coordinates at which to evaluate the /// basis functions. Shape is `(num_points, topological dimension)` - /// (row-major storage) - /// @param[in] shape The shape of `X` - /// @param[in] order The number of derivatives (up to and including - /// this order) to tabulate for + /// (row-major storage). + /// @param[in] shape Shape of `X`. + /// @param[in] order Number of derivatives (up to and including + /// this order) to tabulate for. void tabulate(std::span values, std::span X, std::array shape, int order) const; - /// Evaluate all derivatives of the basis functions up to given order - /// at given points in reference cell + /// @brief Evaluate all derivatives of the basis functions up to given + /// order at given points in reference cell. + /// /// @param[in] X The reference coordinates at which to evaluate the /// basis functions. Shape is `(num_points, topological dimension)` - /// (row-major storage) - /// @param[in] shape The shape of `X` - /// @param[in] order The number of derivatives (up to and including - /// this order) to tabulate for - /// @return Basis function values and array shape (row-major storage) + /// (row-major storage). + /// @param[in] shape Shape of `X`. + /// @param[in] order Number of derivatives (up to and including this + /// order) to tabulate for. + /// @return Basis function values and array shape (row-major storage). std::pair, std::array> tabulate(std::span X, std::array shape, int order) const; /// @brief Number of sub elements (for a mixed or blocked element). - /// @return The number of sub elements + /// @return Number of sub elements. int num_sub_elements() const noexcept; /// @brief Check if element is a mixed element. /// - /// A mixed element i composed of two or more elements of different - /// types (a block element, e.g. a Lagrange element with block size >= - /// 1 is not considered mixed). + /// A mixed element is composed of two or more elements of different + /// types. A blocked element, e.g. a Lagrange element with block size + /// >= 1 is not considered mixed. /// /// @return True if element is mixed. bool is_mixed() const noexcept; @@ -240,15 +328,15 @@ class FiniteElement /// @brief Create a matrix that maps degrees of freedom from one /// element to this element (interpolation). /// - /// @param[in] from The element to interpolate from + /// @param[in] from The element to interpolate from. /// @return Matrix operator that maps the `from` degrees-of-freedom to /// the degrees-of-freedom of this element. The (0) matrix data /// (row-major storage) and (1) the shape (num_dofs of `this` element, /// num_dofs of `from`) are returned. /// /// @pre The two elements must use the same mapping between the - /// reference and physical cells - /// @note Does not support mixed elements + /// reference and physical cells. + /// @note Does not support mixed elements. std::pair, std::array> create_interpolation_operator(const FiniteElement& from) const; @@ -264,7 +352,7 @@ class FiniteElement /// orientation of a basis function, and this orientation cannot be /// corrected for by permuting the DOF numbers on each cell. /// - /// @return True if DOF transformations are required + /// @return True if DOF transformations are required. bool needs_dof_transformations() const noexcept; /// @brief Check if DOF permutations are needed for this element. @@ -280,7 +368,7 @@ class FiniteElement /// this can be corrected for by permuting the DOF numbers on each /// cell. /// - /// @return True if DOF transformations are required + /// @return True if DOF transformations are required. bool needs_dof_permutations() const noexcept; /// @brief Return a function that applies a DOF transformation @@ -321,7 +409,7 @@ class FiniteElement /// data* from the reference element to the conforming (physical) /// ordering, e.g. \f$u = T^{-t} \tilde{u}\f$. /// @param[in] scalar_element Indicates whether the scalar - /// transformations should be returned for a vector element + /// transformations should be returned for a vector element. template std::function, std::span, std::int32_t, int)> @@ -728,8 +816,24 @@ class FiniteElement dof_permutation_fn(bool inverse = false, bool scalar_element = false) const; private: + // Value shape. For blocked elements this is larger than + // _reference_value_shape. For non-blocked 'primal' elements it is + // equal to _reference_value_shape. For mixed elements, it is + // std::nullopt. + std::optional> _value_shape; + + // Block size for BlockedElements. This gives the number of DOFs + // co-located at each dof 'point'. + // TODO: add docs for symmetric elements + int _bs; + + // Element cell shape + mesh::CellType _cell_type; + + // Element signature std::string _signature; + // Dimension of the finite element space (accounting for any blocking) int _space_dim; // List of sub-elements (if any) @@ -740,10 +844,6 @@ class FiniteElement // For a mixed element it is std::nullopt. std::optional> _reference_value_shape; - // Block size for BlockedElements. This gives the number of DOFs - // co-located at each dof 'point'. - int _bs; - // Basix Element (nullptr for mixed elements) std::unique_ptr> _element; @@ -761,4 +861,5 @@ class FiniteElement // all elements except quadrature elements) std::pair, std::array> _points; }; + } // namespace dolfinx::fem diff --git a/cpp/dolfinx/fem/Function.h b/cpp/dolfinx/fem/Function.h index 2ec8298e71c..c2ca69f95b4 100644 --- a/cpp/dolfinx/fem/Function.h +++ b/cpp/dolfinx/fem/Function.h @@ -200,7 +200,8 @@ class Function const auto [fx, fshape] = f(_x); assert(fshape.size() <= 2); - if (int vs = _function_space->value_size(); vs == 1 and fshape.size() == 1) + if (int vs = _function_space->element()->value_size(); + vs == 1 and fshape.size() == 1) { // Check for scalar-valued functions if (fshape.front() != x.size() / 3) @@ -353,7 +354,7 @@ class Function std::size_t value_size = e0.value_size(); if (e0.argument_function_space()) throw std::runtime_error("Cannot interpolate Expression with Argument."); - if (value_size != _function_space->value_size()) + if (value_size != _function_space->element()->value_size()) { throw std::runtime_error( "Function value size not equal to Expression value size."); @@ -485,9 +486,9 @@ class Function auto element = _function_space->element(); assert(element); const int bs_element = element->block_size(); - const std::size_t reference_value_size - = element->reference_value_size() / bs_element; - const std::size_t value_size = _function_space->value_size() / bs_element; + const std::size_t reference_value_size = element->reference_value_size(); + const std::size_t value_size + = _function_space->element()->reference_value_size(); const std::size_t space_dimension = element->space_dimension() / bs_element; // If the space has sub elements, concatenate the evaluations on the diff --git a/cpp/dolfinx/fem/FunctionSpace.h b/cpp/dolfinx/fem/FunctionSpace.h index ea00a6add02..e7b82f9c8a3 100644 --- a/cpp/dolfinx/fem/FunctionSpace.h +++ b/cpp/dolfinx/fem/FunctionSpace.h @@ -36,18 +36,16 @@ class FunctionSpace /// Geometry type of the Mesh that the FunctionSpace is defined on. using geometry_type = T; - /// @brief Create function space for given mesh, element and dofmap. + /// @brief Create function space for given mesh, element and + /// degree-of-freedom map. /// @param[in] mesh Mesh that the space is defined on. /// @param[in] element Finite element for the space. /// @param[in] dofmap Degree-of-freedom map for the space. - /// @param[in] value_shape The shape of the value space on the physical cell FunctionSpace(std::shared_ptr> mesh, std::shared_ptr> element, - std::shared_ptr dofmap, - std::vector value_shape) + std::shared_ptr dofmap) : _mesh(mesh), _element(element), _dofmap(dofmap), - _id(boost::uuids::random_generator()()), _root_space_id(_id), - _value_shape(value_shape) + _id(boost::uuids::random_generator()()), _root_space_id(_id) { // Do nothing } @@ -67,7 +65,7 @@ class FunctionSpace /// Move assignment operator FunctionSpace& operator=(FunctionSpace&& V) = default; - /// @brief Create a subspace (view) for a specific component. + /// @brief Create a subspace (view) for a specific component. /// /// @note If the subspace is re-used, for performance reasons the /// returned subspace should be stored by the caller to avoid repeated @@ -93,10 +91,7 @@ class FunctionSpace = std::make_shared(_dofmap->extract_sub_dofmap(component)); // Create new sub space - FunctionSpace sub_space(_mesh, element, dofmap, - compute_value_shape(element, - _mesh->topology()->dim(), - _mesh->geometry().dim())); + FunctionSpace sub_space(_mesh, element, dofmap); // Set root space id and component w.r.t. root sub_space._root_space_id = _root_space_id; @@ -112,33 +107,24 @@ class FunctionSpace /// FunctionSpace bool contains(const FunctionSpace& V) const { - if (this == std::addressof(V)) - { - // Spaces are the same (same memory address) + if (this == std::addressof(V)) // Spaces are the same (same memory address) return true; - } - else if (_root_space_id != V._root_space_id) - { - // Different root spaces + else if (_root_space_id != V._root_space_id) // Different root spaces return false; - } - else if (_component.size() > V._component.size()) + else if (_component.size() + > V._component.size()) // V is a superspace of *this { - // V is a superspace of *this return false; } else if (!std::equal(_component.begin(), _component.end(), - V._component.begin())) + V._component.begin())) // Components of 'this' are not + // the same as the leading + // components of V { - // Components of 'this' are not the same as the leading components - // of V return false; } - else - { - // Ok, V is really our subspace + else // Ok, V is really our subspace return true; - } } /// Collapse a subspace and return a new function space and a map from @@ -155,11 +141,8 @@ class FunctionSpace auto collapsed_dofmap = std::make_shared(std::move(_collapsed_dofmap)); - return { - FunctionSpace(_mesh, _element, collapsed_dofmap, - compute_value_shape(_element, _mesh->topology()->dim(), - _mesh->geometry().dim())), - std::move(collapsed_dofs)}; + return {FunctionSpace(_mesh, _element, collapsed_dofmap), + std::move(collapsed_dofs)}; } /// @brief Get the component with respect to the root superspace. @@ -168,7 +151,7 @@ class FunctionSpace std::vector component() const { return _component; } /// @brief Indicate whether this function space represents a symmetric - /// 2-tensor + /// 2-tensor. bool symmetric() const { if (_element) @@ -334,23 +317,6 @@ class FunctionSpace /// The dofmap std::shared_ptr dofmap() const { return _dofmap; } - /// The shape of the value space - std::span value_shape() const noexcept - { - return _value_shape; - } - - /// The value size, e.g. 1 for a scalar-valued function, 2 for a 2D vector, 9 - /// for a second-order tensor in 3D. - /// @note The return value of this function is equivalent to - /// `std::accumulate(value_shape().begin(), value_shape().end(), 1, - /// std::multiplies{})`. - int value_size() const - { - return std::accumulate(_value_shape.begin(), _value_shape.end(), 1, - std::multiplies{}); - } - private: // The mesh std::shared_ptr> _mesh; @@ -367,14 +333,13 @@ class FunctionSpace // Unique identifier for the space and for its root space boost::uuids::uuid _id; boost::uuids::uuid _root_space_id; - - std::vector _value_shape; }; -/// Extract FunctionSpaces for (0) rows blocks and (1) columns blocks -/// from a rectangular array of (test, trial) space pairs. The test -/// space must be the same for each row and the trial spaces must be the -/// same for each column. Raises an exception if there is an +/// @brief Extract FunctionSpaces for (0) rows blocks and (1) columns +/// blocks from a rectangular array of (test, trial) space pairs. +/// +/// The test space must be the same for each row and the trial spaces +/// must be the same for each column. Raises an exception if there is an /// inconsistency. e.g. if each form in row i does not have the same /// test space then an exception is raised. /// @@ -430,41 +395,9 @@ common_function_spaces( return {spaces0, spaces1}; } -/// @brief Compute the physical value shape of an element for a mesh -/// @param[in] element The element -/// @param[in] tdim Topological dimension -/// @param[in] gdim Geometric dimension -/// @return Physical valus shape -template -std::vector compute_value_shape( - std::shared_ptr> element, - std::size_t tdim, std::size_t gdim) -{ - auto rvs = element->reference_value_shape(); - std::vector value_shape(rvs.size()); - if (element->block_size() > 1) - { - for (std::size_t i = 0; i < rvs.size(); ++i) - { - value_shape[i] = rvs[i]; - } - } - else - { - for (std::size_t i = 0; i < rvs.size(); ++i) - { - if (rvs[i] == tdim) - value_shape[i] = gdim; - else - value_shape[i] = rvs[i]; - } - } - return value_shape; -} - /// Type deduction -template -FunctionSpace(U mesh, V element, W dofmap, X value_shape) +template +FunctionSpace(U mesh, V element, W dofmap) -> FunctionSpace::type::geometry_type::value_type>; diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index c3ff3838cb8..a1a0a042c50 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -24,7 +24,7 @@ namespace dolfinx::fem::impl { - +/// @brief Typedef using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; @@ -501,9 +501,11 @@ void assemble_matrix( // Integration domain mesh std::shared_ptr> mesh = a.mesh(); assert(mesh); + // Test function mesh auto mesh0 = a.function_spaces().at(0)->mesh(); assert(mesh0); + // Trial function mesh auto mesh1 = a.function_spaces().at(1)->mesh(); assert(mesh1); diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 2f957d0d9f7..59202ffe626 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -20,7 +20,6 @@ namespace dolfinx::fem::impl { - /// Assemble functional over cells template T assemble_cells(mdspan2_t x_dofmap, std::span> x, diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index fb1980f009e..6e7a440fd75 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -28,7 +28,6 @@ namespace dolfinx::fem::impl { - /// @cond using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const std::int32_t, @@ -938,9 +937,11 @@ void lift_bc(std::span b, const Form& a, mdspan2_t x_dofmap, // Integration domain mesh std::shared_ptr> mesh = a.mesh(); assert(mesh); + // Test function mesh auto mesh0 = a.function_spaces().at(0)->mesh(); assert(mesh0); + // Trial function mesh auto mesh1 = a.function_spaces().at(1)->mesh(); assert(mesh1); @@ -1072,12 +1073,13 @@ void lift_bc(std::span b, const Form& a, mdspan2_t x_dofmap, /// @param[in] alpha Scaling to apply template void apply_lifting( - std::span b, const std::vector>> a, + std::span b, + std::vector>>> a, const std::vector>& constants, const std::vector, std::pair, int>>>& coeffs, - const std::vector>>>& - bcs1, + const std::vector< + std::vector>>>& bcs1, const std::vector>& x0, T alpha) { if (!x0.empty() and x0.size() != a.size()) @@ -1099,15 +1101,14 @@ void apply_lifting( if (a[j] and !bcs1[j].empty()) { // Extract data from mesh - std::shared_ptr> mesh = a[j]->mesh(); + std::shared_ptr> mesh = a[j]->get().mesh(); if (!mesh) throw std::runtime_error("Unable to extract a mesh."); mdspan2_t x_dofmap = mesh->geometry().dofmap(); auto x = mesh->geometry().x(); - assert(a[j]->function_spaces().at(0)); - - auto V1 = a[j]->function_spaces()[1]; + assert(a[j]->get().function_spaces().at(0)); + auto V1 = a[j]->get().function_spaces()[1]; assert(V1); auto map1 = V1->dofmap()->index_map; const int bs1 = V1->dofmap()->index_map_bs(); @@ -1115,21 +1116,21 @@ void apply_lifting( const int crange = bs1 * (map1->size_local() + map1->num_ghosts()); bc_markers1.assign(crange, false); bc_values1.assign(crange, 0); - for (const std::shared_ptr>& bc : bcs1[j]) + for (auto& bc : bcs1[j]) { - bc->mark_dofs(bc_markers1); - bc->set(bc_values1, std::nullopt, 1); + bc.get().mark_dofs(bc_markers1); + bc.get().set(bc_values1, std::nullopt, 1); } if (!x0.empty()) { - lift_bc(b, *a[j], x_dofmap, x, constants[j], coeffs[j], bc_values1, - bc_markers1, x0[j], alpha); + lift_bc(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j], + bc_values1, bc_markers1, x0[j], alpha); } else { - lift_bc(b, *a[j], x_dofmap, x, constants[j], coeffs[j], bc_values1, - bc_markers1, std::span(), alpha); + lift_bc(b, a[j]->get(), x_dofmap, x, constants[j], coeffs[j], + bc_values1, bc_markers1, std::span(), alpha); } } } diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index decd5d35411..ab41da82792 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -150,16 +150,17 @@ void assemble_vector(std::span b, const Form& L) /// is responsible for calling VecGhostUpdateBegin/End. template void apply_lifting( - std::span b, const std::vector>>& a, + std::span b, + std::vector>>> a, const std::vector>& constants, const std::vector, std::pair, int>>>& coeffs, - const std::vector>>>& - bcs1, + const std::vector< + std::vector>>>& bcs1, const std::vector>& x0, T alpha) { // If all forms are null, there is nothing to do - if (std::ranges::all_of(a, [](auto ptr) { return ptr == nullptr; })) + if (std::ranges::all_of(a, [](auto ai) { return !ai; })) return; impl::apply_lifting(b, a, constants, coeffs, bcs1, x0, alpha); @@ -179,9 +180,10 @@ void apply_lifting( /// is responsible for calling VecGhostUpdateBegin/End. template void apply_lifting( - std::span b, const std::vector>>& a, - const std::vector>>>& - bcs1, + std::span b, + std::vector>>> a, + const std::vector< + std::vector>>>& bcs1, const std::vector>& x0, T alpha) { std::vector< @@ -192,10 +194,10 @@ void apply_lifting( { if (_a) { - auto coefficients = allocate_coefficient_storage(*_a); - pack_coefficients(*_a, coefficients); + auto coefficients = allocate_coefficient_storage(_a->get()); + pack_coefficients(_a->get(), coefficients); coeffs.push_back(coefficients); - constants.push_back(pack_constants(*_a)); + constants.push_back(pack_constants(_a->get())); } else { @@ -268,7 +270,7 @@ void assemble_matrix( auto mat_add, const Form& a, std::span constants, const std::map, std::pair, int>>& coefficients, - const std::vector>>& bcs) + const std::vector>>& bcs) { // Index maps for dof ranges auto map0 = a.function_spaces().at(0)->dofmap()->index_map; @@ -284,18 +286,17 @@ void assemble_matrix( std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts()); for (std::size_t k = 0; k < bcs.size(); ++k) { - assert(bcs[k]); - assert(bcs[k]->function_space()); - if (a.function_spaces().at(0)->contains(*bcs[k]->function_space())) + assert(bcs[k].get().function_space()); + if (a.function_spaces().at(0)->contains(*bcs[k].get().function_space())) { dof_marker0.resize(dim0, false); - bcs[k]->mark_dofs(dof_marker0); + bcs[k].get().mark_dofs(dof_marker0); } - if (a.function_spaces().at(1)->contains(*bcs[k]->function_space())) + if (a.function_spaces().at(1)->contains(*bcs[k].get().function_space())) { dof_marker1.resize(dim1, false); - bcs[k]->mark_dofs(dof_marker1); + bcs[k].get().mark_dofs(dof_marker1); } } @@ -312,7 +313,7 @@ void assemble_matrix( template void assemble_matrix( auto mat_add, const Form& a, - const std::vector>>& bcs) + const std::vector>>& bcs) { // Prepare constants and coefficients const std::vector constants = pack_constants(a); @@ -392,15 +393,14 @@ void set_diagonal(auto set_fn, std::span rows, template void set_diagonal( auto set_fn, const FunctionSpace& V, - const std::vector>>& bcs, + const std::vector>>& bcs, T diagonal = 1.0) { for (auto& bc : bcs) { - assert(bc); - if (V.contains(*bc->function_space())) + if (V.contains(*bc.get().function_space())) { - const auto [dofs, range] = bc->dof_indices(); + const auto [dofs, range] = bc.get().dof_indices(); set_diagonal(set_fn, dofs.first(range), diagonal); } } diff --git a/cpp/dolfinx/fem/discreteoperators.h b/cpp/dolfinx/fem/discreteoperators.h index 3989e77e9e0..cfe37525a9f 100644 --- a/cpp/dolfinx/fem/discreteoperators.h +++ b/cpp/dolfinx/fem/discreteoperators.h @@ -194,9 +194,9 @@ void interpolation_matrix(const FunctionSpace& V0, const std::size_t space_dim0 = e0->space_dimension(); const std::size_t space_dim1 = e1->space_dimension(); const std::size_t dim0 = space_dim0 / bs0; - const std::size_t value_size_ref0 = e0->reference_value_size() / bs0; - const std::size_t value_size0 = V0.value_size() / bs0; - const std::size_t value_size1 = V1.value_size() / bs1; + const std::size_t value_size_ref0 = e0->reference_value_size(); + const std::size_t value_size0 = V0.element()->reference_value_size(); + const std::size_t value_size1 = V1.element()->reference_value_size(); // Get geometry data const CoordinateElement& cmap = mesh->geometry().cmap(); @@ -267,12 +267,14 @@ void interpolation_matrix(const FunctionSpace& V0, // Basis values of Lagrange space unrolled for block size // (num_quadrature_points, Lagrange dof, value_size) - std::vector basis_values_b(Xshape[0] * bs0 * dim0 * V1.value_size()); + std::vector basis_values_b(Xshape[0] * bs0 * dim0 + * V1.element()->value_size()); mdspan3_t basis_values(basis_values_b.data(), Xshape[0], bs0 * dim0, - V1.value_size()); - std::vector mapped_values_b(Xshape[0] * bs0 * dim0 * V1.value_size()); + V1.element()->value_size()); + std::vector mapped_values_b(Xshape[0] * bs0 * dim0 + * V1.element()->value_size()); mdspan3_t mapped_values(mapped_values_b.data(), Xshape[0], bs0 * dim0, - V1.value_size()); + V1.element()->value_size()); auto pull_back_fn1 = e1->basix_element().template map_fn(); @@ -382,7 +384,7 @@ void interpolation_matrix(const FunctionSpace& V0, { MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> - A(Ab.data(), Xshape[0], V1.value_size(), space_dim0); + A(Ab.data(), Xshape[0], V1.element()->value_size(), space_dim0); for (std::size_t i = 0; i < mapped_values.extent(0); ++i) for (std::size_t j = 0; j < mapped_values.extent(1); ++j) for (std::size_t k = 0; k < mapped_values.extent(2); ++k) diff --git a/cpp/dolfinx/fem/dofmapbuilder.cpp b/cpp/dolfinx/fem/dofmapbuilder.cpp index f7cfcf25d7b..55c6ce3ca21 100644 --- a/cpp/dolfinx/fem/dofmapbuilder.cpp +++ b/cpp/dolfinx/fem/dofmapbuilder.cpp @@ -42,7 +42,7 @@ struct dofmap_t /// [0:owned_size). /// @param[in] dofmaps The local dofmaps (cell -> dofs) /// @param[in] owned_size Number of dofs owned by this process -/// @param[in] original_to_contiguous Map from dof indices in @p dofmap +/// @param[in] original_to_contiguous Map from dof indices in `dofmap`. /// to new indices that are ordered such that owned indices are [0, /// owned_size) /// @param[in] reorder_fn The graph reordering function to apply @@ -370,9 +370,10 @@ build_basic_dofmaps( /// same range /// @param [in] dof_entity Map from dof index to (index_map, entity_index), /// where entity_index is the local mesh entity index in the given index_map -/// @param [in] index_maps The set of IndexMaps, one for each topological -/// entity type used in the dofmap. The location in this array is referred to by -/// the first item in each entry of @p dof_entity +/// @param [in] index_maps The set of IndexMaps, one for each +/// topological entity type used in the dofmap. The location in this +/// array is referred to by the first item in each entry of +/// `dof_entity`. /// @param [in] reorder_fn Graph reordering function that is applied for /// dof re-ordering /// @return The pair (old-to-new local index map, M), where M is the diff --git a/cpp/dolfinx/fem/dofmapbuilder.h b/cpp/dolfinx/fem/dofmapbuilder.h index 28460d8267b..08baf3193c1 100644 --- a/cpp/dolfinx/fem/dofmapbuilder.h +++ b/cpp/dolfinx/fem/dofmapbuilder.h @@ -35,8 +35,8 @@ class ElementDofLayout; /// Build dofmap data for elements on a mesh topology /// @param[in] comm MPI communicator /// @param[in] topology The mesh topology -/// @param[in] element_dof_layouts The element dof layouts for each cell type in -/// @p topology +/// @param[in] element_dof_layouts The element dof layouts for each cell +/// type in `topology`. /// @param[in] reorder_fn Graph reordering function that is applied to /// the dofmaps /// @return The index map, block size, and dofmaps for each element type diff --git a/cpp/dolfinx/fem/interpolate.h b/cpp/dolfinx/fem/interpolate.h index e25fdcbfa2b..d2a098d8c8c 100644 --- a/cpp/dolfinx/fem/interpolate.h +++ b/cpp/dolfinx/fem/interpolate.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -506,8 +507,8 @@ void interpolate_nonmatching_maps(Function& u1, // Get sizes of elements const std::size_t dim0 = element0->space_dimension() / bs0; - const std::size_t value_size_ref0 = element0->reference_value_size() / bs0; - const std::size_t value_size0 = V0->value_size() / bs0; + const std::size_t value_size_ref0 = element0->reference_value_size(); + const std::size_t value_size0 = V0->element()->reference_value_size(); const CoordinateElement& cmap = mesh0->geometry().cmap(); auto x_dofmap = mesh0->geometry().dofmap(); @@ -539,13 +540,13 @@ void interpolate_nonmatching_maps(Function& u1, impl::mdspan_t basis_reference0(basis_reference0_b.data(), Xshape[0], dim0, value_size_ref0); - std::vector values0_b(Xshape[0] * 1 * V1->value_size()); + std::vector values0_b(Xshape[0] * 1 * V1->element()->value_size()); impl::mdspan_t values0(values0_b.data(), Xshape[0], 1, - V1->value_size()); + V1->element()->value_size()); - std::vector mapped_values_b(Xshape[0] * 1 * V1->value_size()); + std::vector mapped_values_b(Xshape[0] * 1 * V1->element()->value_size()); impl::mdspan_t mapped_values0(mapped_values_b.data(), Xshape[0], 1, - V1->value_size()); + V1->element()->value_size()); std::vector coord_dofs_b(num_dofs_g * gdim); impl::mdspan_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim); @@ -729,9 +730,8 @@ void interpolate(Function& u, std::span f, const int gdim = mesh->geometry().dim(); const int tdim = mesh->topology()->dim(); - const bool symmetric = u.function_space()->symmetric(); - if (fshape[0] != (std::size_t)u.function_space()->value_size()) + if (fshape[0] != (std::size_t)u.function_space()->element()->value_size()) throw std::runtime_error("Interpolation data has the wrong shape/size."); std::span cell_info; @@ -741,7 +741,8 @@ void interpolate(Function& u, std::span f, cell_info = std::span(mesh->topology()->get_cell_permutation_info()); } - const std::size_t f_shape1 = f.size() / u.function_space()->value_size(); + const std::size_t f_shape1 + = f.size() / u.function_space()->element()->value_size(); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> _f(f.data(), fshape); @@ -753,13 +754,14 @@ void interpolate(Function& u, std::span f, // Loop over cells and compute interpolation dofs const int num_scalar_dofs = element->space_dimension() / element_bs; - const int value_size = u.function_space()->value_size() / element_bs; + const int value_size = u.function_space()->element()->reference_value_size(); std::span coeffs = u.x()->mutable_array(); std::vector _coeffs(num_scalar_dofs); // This assumes that any element with an identity interpolation matrix // is a point evaluation + const bool symmetric = u.function_space()->symmetric(); if (element->map_ident() && element->interpolation_ident()) { // Point evaluation element *and* the geometric map is the identity, @@ -774,17 +776,18 @@ void interpolate(Function& u, std::span f, std::size_t matrix_size = 0; while (matrix_size * matrix_size < fshape[0]) ++matrix_size; + // Loop over cells for (std::size_t c = 0; c < cells.size(); ++c) { - // The entries of a symmetric matrix are numbered (for an example 4x4 - // element): + // The entries of a symmetric matrix are numbered (for an + // example 4x4 element): // 0 * * * // 1 2 * * // 3 4 5 * // 6 7 8 9 - // The loop extracts these elements. In this loop, row is the row of - // this matrix, and (k - rowstart) is the column + // The loop extracts these elements. In this loop, row is the + // row of this matrix, and (k - rowstart) is the column std::size_t row = 0; std::size_t rowstart = 0; const std::int32_t cell = cells[c]; @@ -796,6 +799,7 @@ void interpolate(Function& u, std::span f, ++row; rowstart = k; } + // num_scalar_dofs is the number of interpolation points per // cell in this case (interpolation matrix is identity) std::copy_n( @@ -847,9 +851,10 @@ void interpolate(Function& u, std::span f, "Interpolation into this element not supported."); } - const int element_vs = u.function_space()->value_size() / element_bs; + const int element_vs + = u.function_space()->element()->reference_value_size(); - if (element_vs > 1 && element_bs > 1) + if (element_vs > 1 and element_bs > 1) { throw std::runtime_error( "Interpolation into this element not supported."); @@ -933,8 +938,7 @@ void interpolate(Function& u, std::span f, std::vector coord_dofs_b(num_dofs_g * gdim); mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim); - const std::size_t value_size_ref - = element->reference_value_size() / element_bs; + const std::size_t value_size_ref = element->reference_value_size(); std::vector ref_data_b(Xshape[0] * 1 * value_size_ref); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> @@ -1132,7 +1136,7 @@ void interpolate(Function& u, const Function& v, assert(cell_map); auto element_u = u.function_space()->element(); assert(element_u); - const std::size_t value_size = u.function_space()->value_size(); + const std::size_t value_size = u.function_space()->element()->value_size(); const std::vector& dest_ranks = interpolation_data.src_owner; const std::vector& src_ranks = interpolation_data.dest_owners; @@ -1217,9 +1221,8 @@ void interpolate(Function& u1, std::span cells1, auto fs1 = u1.function_space(); auto element1 = fs1->element(); assert(element1); - if (fs0->value_shape().size() != fs1->value_shape().size() - or !std::equal(fs0->value_shape().begin(), fs0->value_shape().end(), - fs1->value_shape().begin())) + if (!std::ranges::equal(fs0->element()->value_shape(), + fs1->element()->value_shape())) { throw std::runtime_error( "Interpolation: elements have different value dimensions"); diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h index a94af0b0513..5fbcdcef9ad 100644 --- a/cpp/dolfinx/fem/petsc.h +++ b/cpp/dolfinx/fem/petsc.h @@ -13,6 +13,7 @@ #include "utils.h" #include #include +#include #include #include #include @@ -329,13 +330,17 @@ void assemble_vector(Vec b, const Form& L) /// is responsible for calling VecGhostUpdateBegin/End. template void apply_lifting( - Vec b, const std::vector>>& a, + Vec b, + std::vector< + std::optional>>> + a, const std::vector>& constants, const std::vector, std::pair, int>>>& coeffs, const std::vector< - std::vector>>>& bcs1, + std::vector>>>& + bcs1, const std::vector& x0, PetscScalar alpha) { Vec b_local; @@ -398,10 +403,11 @@ void apply_lifting( template void apply_lifting( Vec b, - const std::vector>>& a, - const std::vector< - std::vector>>>& - bcs1, + std::vector< + std::optional>>> + a, + const std::vector>>>& bcs1, const std::vector& x0, PetscScalar alpha) { Vec b_local; @@ -456,7 +462,7 @@ void apply_lifting( template void set_bc( Vec b, - const std::vector>>& bcs, + const std::vector>> bcs, const Vec x0, PetscScalar alpha = 1) { PetscInt n = 0; @@ -474,7 +480,7 @@ void set_bc( VecGetArrayRead(x0_local, &array); std::span _x0(array, n); for (auto& bc : bcs) - bc->set(_b, _x0, alpha); + bc.set(_b, _x0, alpha); VecRestoreArrayRead(x0_local, &array); VecGhostRestoreLocalForm(x0, &x0_local); } diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 8f098aa686d..d0f7c60a059 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -128,16 +128,13 @@ std::vector fem::get_constant_names(const ufcx_form& ufcx_form) + ufcx_form.num_constants); } //----------------------------------------------------------------------------- -std::vector fem::compute_integration_domains( - fem::IntegralType integral_type, const mesh::Topology& topology, - std::span entities, int dim) +std::vector +fem::compute_integration_domains(fem::IntegralType integral_type, + const mesh::Topology& topology, + std::span entities) { const int tdim = topology.dim(); - if ((integral_type == IntegralType::cell ? tdim : tdim - 1) != dim) - { - throw std::runtime_error("Invalid mesh entity dimension: " - + std::to_string(dim)); - } + const int dim = integral_type == IntegralType::cell ? tdim : tdim - 1; { // Create span of the owned entities (leaves off any ghosts) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 7bd2c89eb0e..d0d5698744a 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -25,10 +25,12 @@ #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -49,7 +51,6 @@ class IndexMap; namespace dolfinx::fem { - namespace impl { /// Helper function to get an array of of (cell, local_facet) pairs @@ -83,37 +84,41 @@ get_cell_facet_pairs(std::int32_t f, std::span cells, } // namespace impl -/// @brief Given an integral type and a set of entities, compute the -/// entities that should be integrated over. +/// @brief Given an integral type and a set of entities, computes and +/// return data for the entities that should be integrated over. /// -/// This function returns a list `[(id, entities)]`. For cell -/// integrals `entities` are the cell indices. For exterior facet -/// integrals, `entities` is a list of `(cell_index, local_facet_index)` -/// pairs. For interior facet integrals, `entities` is a list of -/// `(cell_index0, local_facet_index0, cell_index1, local_facet_index1)`. -/// `id` refers to the subdomain id used in the definition of the integration -/// measures of the variational form. +/// This function returns a list data, for each entity in `entities`, +/// that is used in assembly. For cell integrals it is simply the cell +/// cell indices. For exterior facet integrals, a list of `(cell_index, +/// local_facet_index)` pairs is returned. For interior facet integrals, +/// a list of `(cell_index0, local_facet_index0, cell_index1, +/// local_facet_index1)` tuples is returned. +/// The data computed by this function is typically used as input to +/// fem::create_form. /// /// @note Owned mesh entities only are returned. Ghost entities are not /// included. /// -/// @param[in] integral_type Integral type -/// @param[in] topology Mesh topology -/// @param[in] entities List of mesh entities -/// @param[in] dim Topological dimension of entities -/// @return List of integration entities /// @pre For facet integrals, the topology facet-to-cell and /// cell-to-facet connectivity must be computed before calling this /// function. +/// +/// @param[in] integral_type Integral type. +/// @param[in] topology Mesh topology. +/// @param[in] entities List of mesh entities. For +/// `integral_type==IntegralType::cell`, `entities` should be cell +/// indices. For other `IntegralType`, `entities` should be facet +/// indices. +/// @return List of integration entity data. std::vector compute_integration_domains(IntegralType integral_type, const mesh::Topology& topology, - std::span entities, int dim); + std::span entities); /// @brief Extract test (0) and trial (1) function spaces pairs for each /// bilinear form for a rectangular array of forms. /// -/// @param[in] a A rectangular block on bilinear forms +/// @param[in] a A rectangular block on bilinear forms. /// @return Rectangular array of the same shape as `a` with a pair of /// function spaces in each array entry. If a form is null, then the /// returned function space pair is (null, null). @@ -323,7 +328,8 @@ std::vector get_constant_names(const ufcx_form& ufcx_form); /// equal to the rank of the form. /// @param[in] coefficients Coefficient fields in the form. /// @param[in] constants Spatial constants in the form. -/// @param[in] subdomains Subdomain markers. +/// @param[in] subdomains Subdomain markers. The data can be computed +/// using fem::compute_integration_domains. /// @param[in] entity_maps The entity maps for the form. Empty for /// single domain problems. /// @param[in] mesh The mesh of the domain. @@ -692,7 +698,8 @@ Form create_form_factory( /// @param[in] spaces Function spaces for the Form arguments. /// @param[in] coefficients Coefficient fields in the form (by name). /// @param[in] constants Spatial constants in the form (by name). -/// @param[in] subdomains Subdomain markers. +/// @param[in] subdomains Subdomain markers. The data can be computed +/// using fem::compute_integration_domains. /// @pre Each value in `subdomains` must be sorted by domain id. /// @param[in] entity_maps The entity maps for the form. Empty for /// single domain problems. @@ -728,7 +735,6 @@ Form create_form( } // Place constants in appropriate order - std::vector>> const_map; for (const std::string& name : get_constant_names(ufcx_form)) { @@ -752,7 +758,8 @@ Form create_form( /// @param[in] spaces Function spaces for the Form arguments. /// @param[in] coefficients Coefficient fields in the form (by name), /// @param[in] constants Spatial constants in the form (by name), -/// @param[in] subdomains Subdomain markers. +/// @param[in] subdomains Subdomain markers. The data can be computed +/// using fem::compute_integration_domains. /// @pre Each value in `subdomains` must be sorted by domain id. /// @param[in] entity_maps The entity maps for the form. Empty for /// single domain problems. @@ -781,76 +788,34 @@ Form create_form( return L; } -/// @brief Create a function space from a Basix element. -/// @param[in] mesh Mesh -/// @param[in] e Basix finite element. -/// @param[in] value_shape Value shape for 'blocked' elements, e.g. -/// vector-valued Lagrange elements where each component for the vector -/// field is a Lagrange element. For example, a vector-valued element in -/// 3D will have `value_shape` equal to `{3}`, and for a second-order -/// tensor element in 2D `value_shape` equal to `{2, 2}`. -/// @param[in] reorder_fn The graph reordering function to call on the -/// dofmap. If `nullptr`, the default re-ordering is used. -/// @return The created function space +/// @brief NEW Create a function space from a fem::FiniteElement. template FunctionSpace create_functionspace( - std::shared_ptr> mesh, const basix::FiniteElement& e, - const std::vector& value_shape = {}, + std::shared_ptr> mesh, + std::shared_ptr> e, std::function(const graph::AdjacencyList&)> reorder_fn = nullptr) { - if (!e.value_shape().empty() and !value_shape.empty()) - { - throw std::runtime_error( - "Cannot specify value shape for non-scalar base element."); - } + assert(e); - if (mesh::cell_type_from_basix_type(e.cell_type()) - != mesh->topology()->cell_type()) + // TODO: check cell type of e (need to add method to fem::FiniteElement) + assert(mesh); + assert(mesh->topology()); + if (e->cell_type() != mesh->topology()->cell_type()) throw std::runtime_error("Cell type of element and mesh must match."); - std::size_t bs = value_shape.empty() - ? 1 - : std::accumulate(value_shape.begin(), value_shape.end(), - 1, std::multiplies{}); - - // Create a DOLFINx element - auto _e = std::make_shared>(e, bs); - assert(_e); - - const std::vector _value_shape - = (value_shape.empty() and !e.value_shape().empty()) - ? fem::compute_value_shape(_e, mesh->topology()->dim(), - mesh->geometry().dim()) - : value_shape; - - // Create UFC subdofmaps and compute offset - const int num_sub_elements = _e->num_sub_elements(); - std::vector sub_doflayout; - sub_doflayout.reserve(num_sub_elements); - for (int i = 0; i < num_sub_elements; ++i) - { - auto sub_element = _e->extract_sub_element({i}); - std::vector parent_map_sub(sub_element->space_dimension()); - for (std::size_t j = 0; j < parent_map_sub.size(); ++j) - parent_map_sub[j] = i + _e->block_size() * j; - sub_doflayout.emplace_back(1, e.entity_dofs(), e.entity_closure_dofs(), - parent_map_sub, std::vector()); - } + // Create element dof layout + fem::ElementDofLayout layout = fem::create_element_dof_layout(*e); // Create a dofmap - ElementDofLayout layout(_e->block_size(), e.entity_dofs(), - e.entity_closure_dofs(), {}, sub_doflayout); std::function, std::uint32_t)> permute_inv - = nullptr; - if (_e->needs_dof_permutations()) - permute_inv = _e->dof_permutation_fn(true, true); - assert(mesh); - assert(mesh->topology()); + = e->needs_dof_permutations() ? e->dof_permutation_fn(true, true) + : nullptr; auto dofmap = std::make_shared(create_dofmap( mesh->comm(), layout, *mesh->topology(), permute_inv, reorder_fn)); - return FunctionSpace(mesh, _e, dofmap, _value_shape); + + return FunctionSpace(mesh, e, dofmap); } /// @private diff --git a/cpp/dolfinx/io/ADIOS2Writers.h b/cpp/dolfinx/io/ADIOS2Writers.h index 4ea391ef022..73897339807 100644 --- a/cpp/dolfinx/io/ADIOS2Writers.h +++ b/cpp/dolfinx/io/ADIOS2Writers.h @@ -212,7 +212,8 @@ void vtx_write_data(adios2::IO& io, adios2::Engine& engine, // Pad to 3D if vector/tensor is product of dimensions is smaller than // 3**rank to ensure that we can visualize them correctly in Paraview - std::span value_shape = u.function_space()->value_shape(); + std::span value_shape + = u.function_space()->element()->value_shape(); int rank = value_shape.size(); int num_comp = std::reduce(value_shape.begin(), value_shape.end(), 1, std::multiplies{}); diff --git a/cpp/dolfinx/io/VTKFile.cpp b/cpp/dolfinx/io/VTKFile.cpp index 0f9b7d59b14..927f220c39b 100644 --- a/cpp/dolfinx/io/VTKFile.cpp +++ b/cpp/dolfinx/io/VTKFile.cpp @@ -459,7 +459,7 @@ void write_function( if (piece_node.child(data_type).empty()) piece_node.append_child(data_type); - const int rank = _u.get().function_space()->value_shape().size(); + const int rank = _u.get().function_space()->element()->value_shape().size(); pugi::xml_node data_node = piece_node.child(data_type); if (data_node.attribute(tensor_str[rank]).empty()) data_node.append_attribute(tensor_str[rank]); @@ -476,7 +476,7 @@ void write_function( // Pad to 3D if vector/tensor is product of dimensions is smaller than // 3**rank to ensure that we can visualize them correctly in Paraview - std::span value_shape = V->value_shape(); + std::span value_shape = V->element()->value_shape(); int rank = value_shape.size(); int num_components = std::reduce(value_shape.begin(), value_shape.end(), 1, std::multiplies{}); @@ -660,7 +660,7 @@ void write_function( // Pad to 3D if vector/tensor is product of dimensions is smaller than // 3**rank to ensure that we can visualize them correctly in Paraview - std::span value_shape = V->value_shape(); + std::span value_shape = V->element()->value_shape(); int rank = value_shape.size(); int num_components = std::reduce(value_shape.begin(), value_shape.end(), 1, std::multiplies{}); diff --git a/cpp/dolfinx/io/xdmf_function.cpp b/cpp/dolfinx/io/xdmf_function.cpp index 47ce57ddd6e..c67ede85af7 100644 --- a/cpp/dolfinx/io/xdmf_function.cpp +++ b/cpp/dolfinx/io/xdmf_function.cpp @@ -77,7 +77,8 @@ void xdmf_function::add_function(MPI_Comm comm, const fem::Function& u, // Pad to 3D if vector/tensor is product of dimensions is smaller than 3**rank // to ensure that we can visualize them correctly in Paraview - std::span value_shape = u.function_space()->value_shape(); + std::span value_shape + = u.function_space()->element()->value_shape(); int rank = value_shape.size(); int num_components = std::reduce(value_shape.begin(), value_shape.end(), 1, std::multiplies{}); diff --git a/cpp/dolfinx/mesh/generation.h b/cpp/dolfinx/mesh/generation.h index ba2fb2fd187..d8f004b507d 100644 --- a/cpp/dolfinx/mesh/generation.h +++ b/cpp/dolfinx/mesh/generation.h @@ -299,8 +299,7 @@ std::vector create_geom(MPI_Comm comm, std::array, 2> p, const auto [nx, ny, nz] = n; assert(std::ranges::all_of(n, [](auto e) { return e >= 1; })); - for (std::int64_t i = 0; i < 3; i++) - assert(p0[i] < p1[i]); + assert(p0 < p1); // Structured grid cuboid extents const std::array extents = { diff --git a/cpp/dolfinx/mesh/utils.cpp b/cpp/dolfinx/mesh/utils.cpp index b067a321b60..6a8b9b2fcbc 100644 --- a/cpp/dolfinx/mesh/utils.cpp +++ b/cpp/dolfinx/mesh/utils.cpp @@ -64,11 +64,11 @@ std::vector mesh::exterior_facet_indices(const Topology& topology) throw std::runtime_error( "Facet to cell connectivity has not been computed."); } + // Find all owned facets (not ghost) with only one attached cell auto facet_map = topology.index_map(tdim - 1); - const int num_facets = facet_map->size_local(); std::vector facets; - for (std::int32_t f = 0; f < num_facets; ++f) + for (std::int32_t f = 0; f < facet_map->size_local(); ++f) { if (f_to_c->num_links(f) == 1) facets.push_back(f); diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index 1fa5724e39d..e386381f562 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -59,7 +59,7 @@ void reorder_list(std::span list, std::span nodemap) } } -/// @brief The coordinates of 'vertices' for for entities of a give +/// @brief The coordinates of 'vertices' for for entities of a given /// dimension that are attached to specified facets. /// /// @pre The provided facets must be on the boundary of the mesh. @@ -67,7 +67,7 @@ void reorder_list(std::span list, std::span nodemap) /// @param[in] mesh Mesh to compute the vertex coordinates for /// @param[in] dim Topological dimension of the entities /// @param[in] facets List of facets on the meh boundary -/// @return (0) Entities attached to the boundary facets, (1) vertex +/// @return (0) Entities attached to the boundary facets (sorted), (1) vertex /// coordinates (shape is `(3, num_vertices)`) and (2) map from vertex /// in the full mesh to the position (column) in the vertex coordinates /// array (set to -1 if vertex in full mesh is not in the coordinate @@ -161,7 +161,7 @@ compute_vertex_coords_boundary(const mesh::Mesh& mesh, int dim, /// /// @note Collective /// -/// @param[in] topology Mesh topology +/// @param[in] topology Mesh topology. /// @return Sorted list of owned facet indices that are exterior facets /// of the mesh. std::vector exterior_facet_indices(const Topology& topology); @@ -550,14 +550,14 @@ std::vector locate_entities(const Mesh& mesh, int dim, /// @param[in] marker Marking function, returns `true` for a point that /// is 'marked', and `false` otherwise. /// @returns List of marked entity indices (indices local to the -/// process) +/// process). template U> std::vector locate_entities_boundary(const Mesh& mesh, int dim, U marker) { auto topology = mesh.topology(); assert(topology); - const int tdim = topology->dim(); + int tdim = topology->dim(); if (dim == tdim) { throw std::runtime_error( @@ -567,8 +567,7 @@ std::vector locate_entities_boundary(const Mesh& mesh, int dim, // Compute list of boundary facets mesh.topology_mutable()->create_entities(tdim - 1); mesh.topology_mutable()->create_connectivity(tdim - 1, tdim); - const std::vector boundary_facets - = exterior_facet_indices(*topology); + std::vector boundary_facets = exterior_facet_indices(*topology); using cmdspan3x_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const T, @@ -576,10 +575,10 @@ std::vector locate_entities_boundary(const Mesh& mesh, int dim, std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>; // Run marker function on the vertex coordinates - const auto [facet_entities, xdata, vertex_to_pos] + auto [facet_entities, xdata, vertex_to_pos] = impl::compute_vertex_coords_boundary(mesh, dim, boundary_facets); cmdspan3x_t x(xdata.data(), 3, xdata.size() / 3); - const std::vector marked = marker(x); + std::vector marked = marker(x); if (marked.size() != x.extent(1)) throw std::runtime_error("Length of array of markers is wrong."); diff --git a/cpp/dolfinx/refinement/interval.h b/cpp/dolfinx/refinement/interval.h index a8c633b183d..0b1505acd43 100644 --- a/cpp/dolfinx/refinement/interval.h +++ b/cpp/dolfinx/refinement/interval.h @@ -80,7 +80,7 @@ compute_refinement_data(const mesh::Mesh& mesh, // Mark cells for refinement std::vector> marked_for_update(ranks.size()); - if (cells.has_value()) + if (cells) { std::ranges::for_each(cells.value(), [&](auto cell) @@ -102,7 +102,7 @@ compute_refinement_data(const mesh::Mesh& mesh, // Communicate ghost cells that might have been marked. This is not // necessary for a uniform refinement. - if (cells.has_value()) + if (cells) { update_logical_edgefunction(neighbor_comm, marked_for_update, refinement_marker, *map_c); diff --git a/cpp/dolfinx/refinement/plaza.h b/cpp/dolfinx/refinement/plaza.h index 65cdec06645..f68973e7dd3 100644 --- a/cpp/dolfinx/refinement/plaza.h +++ b/cpp/dolfinx/refinement/plaza.h @@ -502,7 +502,7 @@ compute_refinement_data(const mesh::Mesh& mesh, map_e->size_local() + map_e->num_ghosts(), !edges.has_value()); std::vector> marked_for_update(ranks.size()); - if (edges.has_value()) + if (edges) { for (auto edge : edges.value()) { diff --git a/cpp/test/common/index_map.cpp b/cpp/test/common/index_map.cpp index 3a0fca1c496..128ec7c3f5e 100644 --- a/cpp/test/common/index_map.cpp +++ b/cpp/test/common/index_map.cpp @@ -45,8 +45,9 @@ void test_scatter_fwd(int n) // Scatter values to ghost and check value is correctly received sct.scatter_fwd(data_local, data_ghost); CHECK((int)data_ghost.size() == n * num_ghosts); - CHECK(std::all_of(data_ghost.begin(), data_ghost.end(), [=](auto i) - { return i == val * ((mpi_rank + 1) % mpi_size); })); + CHECK( + std::ranges::all_of(data_ghost, [=](auto i) + { return i == val * ((mpi_rank + 1) % mpi_size); })); std::vector requests = sct.create_request_vector(decltype(sct)::type::p2p); @@ -56,8 +57,9 @@ void test_scatter_fwd(int n) decltype(sct)::type::p2p); sct.scatter_fwd_end(requests); - CHECK(std::all_of(data_ghost.begin(), data_ghost.end(), [=](auto i) - { return i == val * ((mpi_rank + 1) % mpi_size); })); + CHECK( + std::ranges::all_of(data_ghost, [=](auto i) + { return i == val * ((mpi_rank + 1) % mpi_size); })); } void test_scatter_rev() diff --git a/cpp/test/fem/functionspace.cpp b/cpp/test/fem/functionspace.cpp index 77a41481ee6..68c2271932c 100644 --- a/cpp/test/fem/functionspace.cpp +++ b/cpp/test/fem/functionspace.cpp @@ -26,5 +26,6 @@ TEST_CASE("Create Function Space (mismatch of elements)", "[functionspace]") basix::element::lagrange_variant::unset, basix::element::dpc_variant::unset, false); - CHECK_THROWS(fem::create_functionspace(mesh, element, {})); + CHECK_THROWS(fem::create_functionspace( + mesh, std::make_shared>(element))); } diff --git a/cpp/test/io.cpp b/cpp/test/io.cpp index a4c9f8f7bbe..00b20013853 100644 --- a/cpp/test/io.cpp +++ b/cpp/test/io.cpp @@ -36,8 +36,8 @@ void test_vtx_reuse_mesh() basix::element::dpc_variant::unset, false); // Create a scalar function space - auto V = std::make_shared>( - fem::create_functionspace(mesh, e)); + auto V = std::make_shared>(fem::create_functionspace( + mesh, std::make_shared>(e))); // Create a finite element Function auto u = std::make_shared>(V); diff --git a/cpp/test/matrix.cpp b/cpp/test/matrix.cpp index 4cf27c2ef75..9741111964c 100644 --- a/cpp/test/matrix.cpp +++ b/cpp/test/matrix.cpp @@ -120,7 +120,8 @@ la::MatrixCSR create_operator(MPI_Comm comm) basix::element::dpc_variant::unset, false); auto V = std::make_shared>( - fem::create_functionspace(mesh, element, {})); + fem::create_functionspace( + mesh, std::make_shared>(element))); // Prepare and set Constants for the bilinear form auto kappa = std::make_shared>(2.0); @@ -158,7 +159,8 @@ la::MatrixCSR create_operator(MPI_Comm comm) basix::element::dpc_variant::unset, false); auto V = std::make_shared>( - fem::create_functionspace(mesh, element, {})); + fem::create_functionspace( + mesh, std::make_shared>(element))); // Prepare and set Constants for the bilinear form auto kappa = std::make_shared>(2.0); diff --git a/python/demo/demo_stokes.py b/python/demo/demo_stokes.py index d6e366cfde2..a7c1f23d916 100644 --- a/python/demo/demo_stokes.py +++ b/python/demo/demo_stokes.py @@ -152,8 +152,10 @@ def lid_velocity_expression(x): # piecewise linear basis (scalar). -P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type) -P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type) +P2 = element( + "Lagrange", msh.basix_cell(), degree=2, shape=(msh.geometry.dim,), dtype=default_real_type +) +P1 = element("Lagrange", msh.basix_cell(), degree=1, dtype=default_real_type) V, Q = functionspace(msh, P2), functionspace(msh, P1) # Boundary conditions for the velocity field are defined: diff --git a/python/dolfinx/fem/__init__.py b/python/dolfinx/fem/__init__.py index be186a29050..cb5cce7bf20 100644 --- a/python/dolfinx/fem/__init__.py +++ b/python/dolfinx/fem/__init__.py @@ -127,35 +127,39 @@ def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _Matri def compute_integration_domains( - integral_type: IntegralType, topology: Topology, entities: np.ndarray, dim: int + integral_type: IntegralType, topology: Topology, entities: np.ndarray ): """Given an integral type and a set of entities compute integration entities. - This function returns a list `[(id, entities)]`. - For cell integrals `entities` are the cell indices. For exterior facet - integrals, `entities` is a list of `(cell_index, local_facet_index)` - pairs. For interior facet integrals, `entities` is a list of - `(cell_index0, local_facet_index0, cell_index1, local_facet_index1)`. - `id` refers to the subdomain id used in the definition of the integration + This function returns a list `[(id, entities)]`. For cell integrals + `entities` are the cell indices. For exterior facet integrals, + `entities` is a list of `(cell_index, local_facet_index)` pairs. For + interior facet integrals, `entities` is a list of `(cell_index0, + local_facet_index0, cell_index1, local_facet_index1)`. `id` refers + to the subdomain id used in the definition of the integration measures of the variational form. Note: - Owned mesh entities only are returned. Ghost entities are not included. + Owned mesh entities only are returned. Ghost entities are not + included. Note: For facet integrals, the topology facet-to-cell and - cell-to-facet connectivity must be computed before calling this function. + cell-to-facet connectivity must be computed before calling this + function. Args: - integral_type: Integral type - topology: Mesh topology - entities: List of mesh entities - dim: Topological dimension of entities + integral_type: Integral type. + topology: Mesh topology. + entities: List of mesh entities. For + ``integral_type==IntegralType.cell``, ``entities`` should be + cell indices. For other ``IntegralType``s, ``entities`` + should be facet indices. Returns: - List of integration entities + List of integration entities. """ - return _compute_integration_domains(integral_type, topology._cpp_object, entities, dim) + return _compute_integration_domains(integral_type, topology._cpp_object, entities) __all__ = [ diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 68f42bc1dd0..54836cb51a8 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -145,15 +145,13 @@ def get_integration_domains( subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) # type: ignore subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) # type: ignore # Compute integration domains only for each subdomain id in - # the integrals - # If a process has no integral entities, insert an empty - # array + # the integrals. If a process has no integral entities, + # insert an empty array. for id in subdomain_ids: integration_entities = _cpp.fem.compute_integration_domains( integral_type, subdomain._cpp_object.topology, # type: ignore subdomain.find(id), # type: ignore - subdomain.dim, # type: ignore ) domains.append((id, integration_entities)) return [(s[0], np.array(s[1])) for s in domains] diff --git a/python/dolfinx/fem/function.py b/python/dolfinx/fem/function.py index 729294c3bbb..550d994391b 100644 --- a/python/dolfinx/fem/function.py +++ b/python/dolfinx/fem/function.py @@ -561,7 +561,6 @@ class ElementMetaData(typing.NamedTuple): def _create_dolfinx_element( - comm: _MPI.Intracomm, cell_type: _cpp.mesh.CellType, ufl_e: ufl.FiniteElementBase, dtype: np.dtype, @@ -575,15 +574,18 @@ def _create_dolfinx_element( raise ValueError(f"Unsupported dtype: {dtype}") if ufl_e.is_mixed: - elements = [_create_dolfinx_element(comm, cell_type, e, dtype) for e in ufl_e.sub_elements] + elements = [_create_dolfinx_element(cell_type, e, dtype) for e in ufl_e.sub_elements] return CppElement(elements) elif ufl_e.is_quadrature: return CppElement( - cell_type, ufl_e.custom_quadrature()[0], ufl_e.block_size, ufl_e.is_symmetric + cell_type, ufl_e.custom_quadrature()[0], [ufl_e.block_size], ufl_e.is_symmetric ) else: basix_e = ufl_e.basix_element._e - return CppElement(basix_e, ufl_e.block_size, ufl_e.is_symmetric) + bs = ufl_e.reference_value_shape if ufl_e.block_size > 1 else None + # bs = [ufl_e.block_size] if ufl_e.block_size > 1 else None + # print(bs, ufl_e.reference_value_shape) + return CppElement(basix_e, bs, ufl_e.is_symmetric) def functionspace( @@ -617,16 +619,16 @@ def functionspace( if ufl_e.cell != mesh.ufl_domain().ufl_cell(): raise ValueError("Non-matching UFL cell and mesh cell shapes.") - ufl_space = ufl.FunctionSpace(mesh.ufl_domain(), ufl_e) - value_shape = ufl_space.value_shape + # ufl_space = ufl.FunctionSpace(mesh.ufl_domain(), ufl_e) + # value_shape = ufl_space.value_shape # Compile dofmap and element and create DOLFINx objects if form_compiler_options is None: form_compiler_options = dict() form_compiler_options["scalar_type"] = dtype - cpp_element = _create_dolfinx_element(mesh.comm, mesh.topology.cell_type, ufl_e, dtype) - + print("Create element") + cpp_element = _create_dolfinx_element(mesh.topology.cell_type, ufl_e, dtype) cpp_dofmap = _cpp.fem.create_dofmap(mesh.comm, mesh.topology._cpp_object, cpp_element) assert np.issubdtype( @@ -635,13 +637,9 @@ def functionspace( # Initialize the cpp.FunctionSpace try: - cppV = _cpp.fem.FunctionSpace_float64( - mesh._cpp_object, cpp_element, cpp_dofmap, value_shape - ) + cppV = _cpp.fem.FunctionSpace_float64(mesh._cpp_object, cpp_element, cpp_dofmap) except TypeError: - cppV = _cpp.fem.FunctionSpace_float32( - mesh._cpp_object, cpp_element, cpp_dofmap, value_shape - ) + cppV = _cpp.fem.FunctionSpace_float32(mesh._cpp_object, cpp_element, cpp_dofmap) return FunctionSpace(mesh, ufl_e, cppV) @@ -695,17 +693,11 @@ def clone(self) -> FunctionSpace: """ try: Vcpp = _cpp.fem.FunctionSpace_float64( - self._cpp_object.mesh, - self._cpp_object.element, - self._cpp_object.dofmap, - self._cpp_object.value_shape, + self._cpp_object.mesh, self._cpp_object.element, self._cpp_object.dofmap ) # type: ignore except TypeError: Vcpp = _cpp.fem.FunctionSpace_float32( - self._cpp_object.mesh, - self._cpp_object.element, - self._cpp_object.dofmap, - self._cpp_object.value_shape, + self._cpp_object.mesh, self._cpp_object.element, self._cpp_object.dofmap ) # type: ignore return FunctionSpace(self._mesh, self.ufl_element(), Vcpp) @@ -714,10 +706,10 @@ def num_sub_spaces(self) -> int: """Number of sub spaces.""" return self.element.num_sub_elements - @property - def value_shape(self) -> tuple[int, ...]: - """Value shape.""" - return tuple(int(i) for i in self._cpp_object.value_shape) + # @property + # def value_shape(self) -> tuple[int, ...]: + # """Value shape.""" + # return tuple(int(i) for i in self._cpp_object.value_shape) def sub(self, i: int) -> FunctionSpace: """Return the i-th sub space. diff --git a/python/dolfinx/wrappers/assemble.cpp b/python/dolfinx/wrappers/assemble.cpp index cab064b6d8e..2d85dd7501f 100644 --- a/python/dolfinx/wrappers/assemble.cpp +++ b/python/dolfinx/wrappers/assemble.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -87,7 +88,7 @@ void declare_discrete_operators(nb::module_& m) const dolfinx::fem::FunctionSpace& V1) { // Create sparsity - auto sp = create_sparsity(V0, V1); + dolfinx::la::SparsityPattern sp = create_sparsity(V0, V1); // Build operator dolfinx::la::MatrixCSR A(sp); @@ -143,7 +144,7 @@ void declare_discrete_operators(nb::module_& m) [](const dolfinx::fem::FunctionSpace& V0, const dolfinx::fem::FunctionSpace& V1) { - auto sp = create_sparsity(V0, V1); + dolfinx::la::SparsityPattern sp = create_sparsity(V0, V1); // Build operator dolfinx::la::MatrixCSR A(sp); @@ -246,9 +247,17 @@ void declare_assembly_functions(nb::module_& m) const std::map, nb::ndarray, nb::c_contig>>& coefficients, - const std::vector< - std::shared_ptr>>& bcs) + std::vector*> bcs) { + std::vector< + std::reference_wrapper>> + _bcs; + for (auto bc : bcs) + { + assert(bc); + _bcs.push_back(*bc); + } + const std::array data_bs = {a.function_spaces().at(0)->dofmap()->index_map_bs(), a.function_spaces().at(1)->dofmap()->index_map_bs()}; @@ -264,63 +273,63 @@ void declare_assembly_functions(nb::module_& m) dolfinx::fem::assemble_matrix( A.mat_add_values(), a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 2) { auto mat_add = A.template mat_add_values<2, 2>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 3) { auto mat_add = A.template mat_add_values<3, 3>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 4) { auto mat_add = A.template mat_add_values<4, 4>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 5) { auto mat_add = A.template mat_add_values<5, 5>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 6) { auto mat_add = A.template mat_add_values<6, 6>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 7) { auto mat_add = A.template mat_add_values<7, 7>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 8) { auto mat_add = A.template mat_add_values<8, 8>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else if (data_bs[0] == 9) { auto mat_add = A.template mat_add_values<9, 9>(); dolfinx::fem::assemble_matrix( mat_add, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else throw std::runtime_error("Block size not supported in Python"); @@ -330,12 +339,19 @@ void declare_assembly_functions(nb::module_& m) m.def( "insert_diagonal", [](dolfinx::la::MatrixCSR& A, const dolfinx::fem::FunctionSpace& V, - const std::vector< - std::shared_ptr>>& bcs, - T diagonal) + std::vector*> bcs, T diagonal) { + std::vector< + std::reference_wrapper>> + _bcs; + for (auto bc : bcs) + { + assert(bc); + _bcs.push_back(*bc); + } + // NB block size of data ("diagonal") is (1, 1) - dolfinx::fem::set_diagonal(A.mat_set_values(), V, bcs, diagonal); + dolfinx::fem::set_diagonal(A.mat_set_values(), V, _bcs, diagonal); }, nb::arg("A"), nb::arg("V"), nb::arg("bcs"), nb::arg("diagonal"), "Experimental."); @@ -359,9 +375,17 @@ void declare_assembly_functions(nb::module_& m) nb::ndarray, nb::c_contig, nb::numpy>)> fin, const dolfinx::fem::Form& form, - const std::vector< - std::shared_ptr>>& bcs) + std::vector*> bcs) { + std::vector< + std::reference_wrapper>> + _bcs; + for (auto bc : bcs) + { + assert(bc); + _bcs.push_back(*bc); + } + auto f = [&fin](std::span rows, std::span cols, std::span data) @@ -373,7 +397,7 @@ void declare_assembly_functions(nb::module_& m) nb::ndarray, nb::c_contig, nb::numpy>( data.data(), {rows.size(), cols.size()})); }; - dolfinx::fem::assemble_matrix(f, form, bcs); + dolfinx::fem::assemble_matrix(f, form, _bcs); }, nb::arg("fin"), nb::arg("form"), nb::arg("bcs"), "Experimental assembly with Python insertion function. This will be " @@ -383,17 +407,40 @@ void declare_assembly_functions(nb::module_& m) m.def( "apply_lifting", [](nb::ndarray, nb::c_contig> b, - const std::vector>>& a, + std::vector*> a, const std::vector, nb::c_contig>>& constants, const std::vector< std::map, nb::ndarray, nb::c_contig>>>& coeffs, - const std::vector>>>& bcs1, + std::vector*>> bcs1, const std::vector, nb::c_contig>>& x0, T alpha) { + std::vector>>> + _bcs; + for (auto& bc1 : bcs1) + { + auto& _bcs0 = _bcs.emplace_back(); + for (auto bc : bc1) + { + assert(bc); + _bcs0.push_back(*bc); + } + } + + std::vector>>> + _a; + for (auto form : a) + { + if (form) + _a.push_back(*form); + else + _a.push_back(std::nullopt); + } + std::vector> _x0; for (auto x : x0) _x0.emplace_back(x.data(), x.size()); @@ -409,8 +456,8 @@ void declare_assembly_functions(nb::module_& m) std::ranges::transform(coeffs, std::back_inserter(_coeffs), [](auto& c) { return py_to_cpp_coeffs(c); }); - dolfinx::fem::apply_lifting(std::span(b.data(), b.size()), a, - _constants, _coeffs, bcs1, _x0, alpha); + dolfinx::fem::apply_lifting(std::span(b.data(), b.size()), _a, + _constants, _coeffs, _bcs, _x0, alpha); }, nb::arg("b").noconvert(), nb::arg("a"), nb::arg("constants"), nb::arg("coeffs"), nb::arg("bcs1"), nb::arg("x0"), nb::arg("alpha"), diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index bba5e0574f1..876e85029ef 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -72,10 +72,8 @@ void declare_function_space(nb::module_& m, std::string type) "Finite element function space") .def(nb::init>, std::shared_ptr>, - std::shared_ptr, - std::vector>(), - nb::arg("mesh"), nb::arg("element"), nb::arg("dofmap"), - nb::arg("value_shape")) + std::shared_ptr>(), + nb::arg("mesh"), nb::arg("element"), nb::arg("dofmap")) .def("collapse", &dolfinx::fem::FunctionSpace::collapse) .def("component", &dolfinx::fem::FunctionSpace::component) .def("contains", &dolfinx::fem::FunctionSpace::contains, @@ -83,15 +81,6 @@ void declare_function_space(nb::module_& m, std::string type) .def_prop_ro("element", &dolfinx::fem::FunctionSpace::element) .def_prop_ro("mesh", &dolfinx::fem::FunctionSpace::mesh) .def_prop_ro("dofmap", &dolfinx::fem::FunctionSpace::dofmap) - .def_prop_ro( - "value_shape", - [](const dolfinx::fem::FunctionSpace& self) - { - std::span vshape = self.value_shape(); - return nb::ndarray(vshape.data(), - {vshape.size()}); - }, - nb::rv_policy::reference_internal) .def("sub", &dolfinx::fem::FunctionSpace::sub, nb::arg("component")) .def("tabulate_dof_coordinates", [](const dolfinx::fem::FunctionSpace& self) @@ -109,13 +98,15 @@ void declare_function_space(nb::module_& m, std::string type) .def( "__init__", [](dolfinx::fem::FiniteElement* self, - basix::FiniteElement& element, std::size_t block_size, + basix::FiniteElement& element, + std::optional> block_shape, bool symmetric) { - new (self) dolfinx::fem::FiniteElement(element, block_size, + new (self) dolfinx::fem::FiniteElement(element, block_shape, symmetric); }, - nb::arg("element"), nb::arg("block_size"), nb::arg("symmetric")) + nb::arg("element"), nb::arg("block_shape").none(), + nb::arg("symmetric"), "Single Basix element constructor.") .def( "__init__", [](dolfinx::fem::FiniteElement* self, @@ -123,20 +114,20 @@ void declare_function_space(nb::module_& m, std::string type) std::shared_ptr>> elements) { new (self) dolfinx::fem::FiniteElement(elements); }, - nb::arg("elements")) + nb::arg("elements"), "Mixed-element constructor.") .def( "__init__", [](dolfinx::fem::FiniteElement* self, mesh::CellType cell_type, nb::ndarray, nb::numpy> points, - std::size_t block_size, bool symmetry) + std::vector block_shape, bool symmetry) { std::span pdata(points.data(), points.size()); new (self) dolfinx::fem::FiniteElement( cell_type, pdata, {points.shape(0), points.shape(1)}, - block_size, symmetry); + block_shape, symmetry); }, - nb::arg("cell_type"), nb::arg("points"), nb::arg("block_size"), - nb::arg("symmetry")) + nb::arg("cell_type"), nb::arg("points"), nb::arg("block_shape"), + nb::arg("symmetry"), "Quadrature element constructor.") .def("__eq__", &dolfinx::fem::FiniteElement::operator==) .def_prop_ro("dtype", [](const dolfinx::fem::FiniteElement&) { return dolfinx_wrappers::numpy_dtype(); }) @@ -145,6 +136,15 @@ void declare_function_space(nb::module_& m, std::string type) nb::rv_policy::reference_internal) .def_prop_ro("num_sub_elements", &dolfinx::fem::FiniteElement::num_sub_elements) + .def_prop_ro( + "value_shape", + [](const dolfinx::fem::FiniteElement& self) + { + std::span vshape = self.value_shape(); + return nb::ndarray(vshape.data(), + {vshape.size()}); + }, + nb::rv_policy::reference_internal) .def("interpolation_points", [](const dolfinx::fem::FiniteElement& self) { @@ -466,7 +466,7 @@ void declare_objects(nb::module_& m, const std::string& type) const int gdim = self.function_space()->mesh()->geometry().dim(); // Compute value size - auto vshape = self.function_space()->value_shape(); + auto vshape = self.function_space()->element()->value_shape(); std::size_t value_size = std::reduce(vshape.begin(), vshape.end(), 1, std::multiplies{}); @@ -1197,16 +1197,13 @@ void fem(nb::module_& m) "compute_integration_domains", [](dolfinx::fem::IntegralType type, const dolfinx::mesh::Topology& topology, - nb::ndarray, nb::c_contig> entities, - int dim) + nb::ndarray, nb::c_contig> entities) { return dolfinx_wrappers::as_nbarray( dolfinx::fem::compute_integration_domains( - type, topology, std::span(entities.data(), entities.size()), - dim)); + type, topology, std::span(entities.data(), entities.size()))); }, - nb::arg("integral_type"), nb::arg("topology"), nb::arg("entities"), - nb::arg("dim")); + nb::arg("integral_type"), nb::arg("topology"), nb::arg("entities")); // dolfinx::fem::ElementDofLayout nb::class_( diff --git a/python/dolfinx/wrappers/petsc.cpp b/python/dolfinx/wrappers/petsc.cpp index 5152b42b7a0..f0841f54854 100644 --- a/python/dolfinx/wrappers/petsc.cpp +++ b/python/dolfinx/wrappers/petsc.cpp @@ -147,11 +147,12 @@ void petsc_la_module(nb::module_& m) "create_index_sets", [](const std::vector>& maps) { - std::vector< - std::pair, int>> - _maps; - for (auto m : maps) - _maps.push_back({*m.first, m.second}); + using X = std::vector< + std::pair, int>>; + X _maps; + std::ranges::transform(maps, std::back_inserter(_maps), + [](auto m) -> typename X::value_type + { return {*m.first, m.second}; }); std::vector index_sets = dolfinx::la::petsc::create_index_sets(_maps); @@ -175,14 +176,15 @@ void petsc_la_module(nb::module_& m) std::shared_ptr, int>>& maps) { std::vector> _x_b; - std::vector, int>> - _maps; - for (auto& array : x_b) - _x_b.emplace_back(array.data(), array.size()); - for (auto q : maps) - _maps.push_back({*q.first, q.second}); + std::ranges::transform(x_b, std::back_inserter(_x_b), [](auto x) + { return std::span(x.data(), x.size()); }); + using X = std::vector, int>>; + X _maps; + std::ranges::transform(maps, std::back_inserter(_maps), + [](auto q) -> typename X::value_type + { return {*q.first, q.second}; }); dolfinx::la::petsc::scatter_local_vectors(x, _x_b, _maps); }, nb::arg("x"), nb::arg("x_b"), nb::arg("maps"), @@ -195,16 +197,19 @@ void petsc_la_module(nb::module_& m) const std::vector, int>>& maps) { - std::vector, int>> - _maps; - for (auto m : maps) - _maps.push_back({*m.first, m.second}); + using X = std::vector, int>>; + X _maps; + std::ranges::transform(maps, std::back_inserter(_maps), + [](auto& m) -> typename X::value_type + { return {*m.first, m.second}; }); + std::vector> vecs = dolfinx::la::petsc::get_local_vectors(x, _maps); std::vector> ret; - for (std::vector& v : vecs) - ret.push_back(dolfinx_wrappers::as_nbarray(std::move(v))); + std::ranges::transform( + vecs, std::back_inserter(ret), + [](auto& v) { return dolfinx_wrappers::as_nbarray(std::move(v)); }); return ret; }, nb::arg("x"), nb::arg("maps"), @@ -219,12 +224,12 @@ void petsc_fem_module(nb::module_& m) [](const std::vector< std::pair, int>>& maps) { - std::vector< - std::pair, int>> - _maps; - for (auto q : maps) - _maps.push_back({*q.first, q.second}); - + using X = std::vector< + std::pair, int>>; + X _maps; + std::ranges::transform(maps, std::back_inserter(_maps), + [](auto q) -> typename X::value_type + { return {*q.first, q.second}; }); return dolfinx::fem::petsc::create_vector_block(_maps); }, nb::rv_policy::take_ownership, nb::arg("maps"), @@ -234,11 +239,12 @@ void petsc_fem_module(nb::module_& m) [](const std::vector< std::pair, int>>& maps) { - std::vector< - std::pair, int>> - _maps; - for (auto m : maps) - _maps.push_back({*m.first, m.second}); + using X = std::vector< + std::pair, int>>; + X _maps; + std::ranges::transform(maps, std::back_inserter(_maps), + [](auto m) -> typename X::value_type + { return {*m.first, m.second}; }); return dolfinx::fem::petsc::create_vector_nest(_maps); }, nb::rv_policy::take_ownership, nb::arg("maps"), @@ -266,10 +272,19 @@ void petsc_fem_module(nb::module_& m) const std::map, nb::ndarray, nb::c_contig>>& coefficients, - const std::vector>>& bcs, + std::vector*> + bcs, bool unrolled) { + std::vector>> + _bcs; + for (auto bc : bcs) + { + assert(bc); + _bcs.push_back(*bc); + } + if (unrolled) { auto set_fn = dolfinx::la::petsc::Matrix::set_block_expand_fn( @@ -277,14 +292,14 @@ void petsc_fem_module(nb::module_& m) a.function_spaces()[1]->dofmap()->bs(), ADD_VALUES); dolfinx::fem::assemble_matrix( set_fn, a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } else { dolfinx::fem::assemble_matrix( dolfinx::la::petsc::Matrix::set_block_fn(A, ADD_VALUES), a, std::span(constants.data(), constants.size()), - py_to_cpp_coeffs(coefficients), bcs); + py_to_cpp_coeffs(coefficients), _bcs); } }, nb::arg("A"), nb::arg("a"), nb::arg("constants"), nb::arg("coeffs"), @@ -325,12 +340,21 @@ void petsc_fem_module(nb::module_& m) m.def( "insert_diagonal", [](Mat A, const dolfinx::fem::FunctionSpace& V, - const std::vector>>& bcs, + std::vector*> + bcs, PetscScalar diagonal) { + std::vector>> + _bcs; + for (auto bc : bcs) + { + assert(bc); + _bcs.push_back(*bc); + } + dolfinx::fem::set_diagonal( - dolfinx::la::petsc::Matrix::set_fn(A, INSERT_VALUES), V, bcs, + dolfinx::la::petsc::Matrix::set_fn(A, INSERT_VALUES), V, _bcs, diagonal); }, nb::arg("A"), nb::arg("V"), nb::arg("bcs"), nb::arg("diagonal")); diff --git a/python/test/unit/fem/test_assemble_mesh_independent_form.py b/python/test/unit/fem/test_assemble_mesh_independent_form.py index aacd934bf7e..ce8895bba22 100644 --- a/python/test/unit/fem/test_assemble_mesh_independent_form.py +++ b/python/test/unit/fem/test_assemble_mesh_independent_form.py @@ -122,10 +122,7 @@ def g(x): wh.interpolate(g) facet_entities = dolfinx.fem.compute_integration_domains( - dolfinx.fem.IntegralType.exterior_facet, - mesh.topology, - sub_to_parent, - mesh.topology.dim - 1, + dolfinx.fem.IntegralType.exterior_facet, mesh.topology, sub_to_parent ) subdomains = {dolfinx.fem.IntegralType.exterior_facet: [(subdomain_id, facet_entities)]} form = dolfinx.fem.create_form( diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index e07c91178d1..a861d6c2b92 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -456,8 +456,9 @@ def compute_mapped_interior_facet_data(mesh, facet_tag, value, parent_to_sub_map Integration data for interior facets """ mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) + assert facet_tag.dim == mesh.topology.dim - 1 integration_data = compute_integration_domains( - fem.IntegralType.interior_facet, mesh.topology, facet_tag.find(value), facet_tag.dim + fem.IntegralType.interior_facet, mesh.topology, facet_tag.find(value) ) mapped_cell_0 = parent_to_sub_map[integration_data[0::4]] mapped_cell_1 = parent_to_sub_map[integration_data[2::4]] diff --git a/python/test/unit/fem/test_mixed_mesh_dofmap.py b/python/test/unit/fem/test_mixed_mesh_dofmap.py index b3745feb398..38b9e68fcf6 100644 --- a/python/test/unit/fem/test_mixed_mesh_dofmap.py +++ b/python/test/unit/fem/test_mixed_mesh_dofmap.py @@ -15,7 +15,7 @@ def create_element_dofmap(mesh, cell_types, degree): cpp_elements = [] for cell_type in cell_types: ufl_e = basix.ufl.element("P", cell_type, degree, dtype=np.float64) - cpp_elements += [_cpp.fem.FiniteElement_float64(ufl_e.basix_element._e, 1, False)] + cpp_elements += [_cpp.fem.FiniteElement_float64(ufl_e.basix_element._e, None, False)] cpp_dofmaps = _cpp.fem.create_dofmaps(mesh.comm, mesh.topology, cpp_elements) diff --git a/python/test/unit/fem/test_symmetric_tensor.py b/python/test/unit/fem/test_symmetric_tensor.py index 588890a3d58..4831d19d0aa 100644 --- a/python/test/unit/fem/test_symmetric_tensor.py +++ b/python/test/unit/fem/test_symmetric_tensor.py @@ -20,9 +20,7 @@ def test_transpose(degree, symmetry): symmetry=symmetry, dtype=dolfinx.default_real_type, ) - space = dolfinx.fem.functionspace(mesh, e) - f = dolfinx.fem.Function(space) f.interpolate(lambda x: (x[0], x[1], x[0] ** 3, x[0])) @@ -39,17 +37,17 @@ def tensor(x): return np.broadcast_to(mat, (9, x.shape[1])) element = basix.ufl.element( - "DG", mesh.basix_cell(), 0, shape=(3, 3), dtype=dolfinx.default_real_type + "DG", mesh.basix_cell(), 0, shape=(3, 3), symmetry=False, dtype=dolfinx.default_real_type ) + space = dolfinx.fem.functionspace(mesh, element) + f = dolfinx.fem.Function(space) + f.interpolate(lambda x: tensor(x)) + symm_element = basix.ufl.element( "DG", mesh.basix_cell(), 0, shape=(3, 3), symmetry=True, dtype=dolfinx.default_real_type ) - space = dolfinx.fem.functionspace(mesh, element) symm_space = dolfinx.fem.functionspace(mesh, symm_element) - f = dolfinx.fem.Function(space) symm_f = dolfinx.fem.Function(symm_space) - - f.interpolate(lambda x: tensor(x)) symm_f.interpolate(lambda x: tensor(x)) l2_error = dolfinx.fem.assemble_scalar(dolfinx.fem.form((f - symm_f) ** 2 * ufl.dx)) @@ -71,10 +69,7 @@ def tensor(x): ) space = dolfinx.fem.functionspace(mesh, element) f = dolfinx.fem.Function(space) - f.interpolate(lambda x: tensor(x)) - value = f.eval([[0, 0, 0]], [0]) - atol = 10 * np.finfo(dolfinx.default_scalar_type).resolution assert np.allclose(value, mat, atol=atol)