From c5953acb71044832a3a8ceaddaeaaa53a4b71b3a Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 21:18:32 +0000 Subject: [PATCH 01/30] Start reworking mpc with nanobind. Some tests pass. Currently debugging periodic conditions --- Changelog.md | 7 +- README.md | 2 +- cpp/ContactConstraint.h | 86 ++++---- cpp/MultiPointConstraint.h | 16 +- cpp/assemble_matrix.cpp | 4 +- cpp/assemble_vector.cpp | 3 +- cpp/lifting.h | 5 +- cpp/mpc_helpers.h | 36 ++-- cpp/utils.h | 4 +- python/CMakeLists.txt | 102 ++++++--- python/README.md | 8 +- python/dolfinx_mpc/dictcondition.py | 2 +- python/dolfinx_mpc/dolfinx_mpc.cpp | 15 +- python/dolfinx_mpc/mpc.cpp | 255 ++++++++++++----------- python/dolfinx_mpc/utils/test.py | 11 +- python/pyproject.toml | 53 +++++ python/setup.cfg | 10 - python/setup.py | 94 --------- python/tests/test_integration_domains.py | 7 +- 19 files changed, 371 insertions(+), 349 deletions(-) create mode 100644 python/pyproject.toml delete mode 100644 python/setup.cfg delete mode 100644 python/setup.py diff --git a/Changelog.md b/Changelog.md index e3df5e6d..7ad81edd 100644 --- a/Changelog.md +++ b/Changelog.md @@ -1,6 +1,11 @@ # Changelog ## Main + - **API** + - Various shared pointers in C++ interface is changed to const references + - Multipoint-constraint now accept `std::span` instead of vectors + - Now using [nanobind](https://github.com/wjakob/nanobind) for Python bindings + - Switch to `pyproject.toml`, **see installation notes** for updated instructions - **DOLFINx API-changes** - `dolfinx.fem.FunctionSpaceBase` replaced by `dolfinx.fem.FunctionSpace` - `ufl.FiniteElement` and `ufl.VectorElement` is replaced by `basix.ufl.element` @@ -10,7 +15,7 @@ - **API**: - Change input of `dolfinx_mpc.MultiPointConstraint.homogenize` and `dolfinx_mpc.backsubstitution` to `dolfinx.fem.Function` instead of `PETSc.Vec`. - **New feature**: Add support for more floating types (float32, float64, complex64, complex128). The floating type of a MPC is related to the mesh geometry. - - This resulted in a minor refactoring of the pybindings, meaning that tte class `dolfinx_mpc.cpp.mpc.MultiPointConstraint` is replaced by `dolfinx_mpc.cpp.mpc.MultiPointConstraint_{dtype}` + - This resulted in a minor refactoring of the pybindings, meaning that the class `dolfinx_mpc.cpp.mpc.MultiPointConstraint` is replaced by `dolfinx_mpc.cpp.mpc.MultiPointConstraint_{dtype}` - Casting scalar-type with `dolfinx.default_scalar_type` instead of `PETSc.ScalarType` - Remove usage of `VectorFunctionSpace`. Use blocked basix element instead. - **DOLFINX API-changes**: diff --git a/README.md b/README.md index 6f105085..c520bcc1 100644 --- a/README.md +++ b/README.md @@ -61,5 +61,5 @@ To install the `dolfinx_mpc`-library run the following code from this directory: ```bash cmake -G Ninja -DCMAKE_BUILD_TYPE=Release -B build-dir cpp/ ninja -j3 install -C build-dir -python3 -m pip install python/. --upgrade +python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./python -U ``` \ No newline at end of file diff --git a/cpp/ContactConstraint.h b/cpp/ContactConstraint.h index 84574f59..9d3ffff5 100644 --- a/cpp/ContactConstraint.h +++ b/cpp/ContactConstraint.h @@ -63,17 +63,16 @@ dolfinx_mpc::mpc_data compute_master_contributions( U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> normals, - std::shared_ptr> V, + const dolfinx::fem::FunctionSpace& V, MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> tabulated_basis_values) { const double tol = 1e-6; - auto mesh = V->mesh(); - const std::int32_t block_size = V->dofmap()->index_map_bs(); - std::shared_ptr imap - = V->dofmap()->index_map; - const int bs = V->dofmap()->index_map_bs(); + auto mesh = V.mesh(); + const std::int32_t block_size = V.dofmap()->index_map_bs(); + std::shared_ptr imap = V.dofmap()->index_map; + const int bs = V.dofmap()->index_map_bs(); const std::int32_t size_local = imap->size_local(); MPI_Comm comm = mesh->comm(); int rank = -1; @@ -89,7 +88,7 @@ dolfinx_mpc::mpc_data compute_master_contributions( { if (const std::int32_t cell = local_colliding_cell[i]; cell != -1) { - auto cell_blocks = V->dofmap()->cell_dofs(cell); + auto cell_blocks = V.dofmap()->cell_dofs(cell); for (std::size_t j = 0; j < cell_blocks.size(); ++j) { for (int b = 0; b < bs; b++) @@ -124,7 +123,7 @@ dolfinx_mpc::mpc_data compute_master_contributions( { if (const std::int32_t cell = local_colliding_cell[i]; cell != -1) { - auto cell_blocks = V->dofmap()->cell_dofs(cell); + auto cell_blocks = V.dofmap()->cell_dofs(cell); global_blocks.resize(cell_blocks.size()); imap->local_to_global(cell_blocks, global_blocks); // Compute coefficients for each master @@ -169,7 +168,7 @@ dolfinx_mpc::mpc_data compute_master_contributions( /// tagged with the marker template std::vector -locate_slave_dofs(std::shared_ptr> V, +locate_slave_dofs(const dolfinx::fem::FunctionSpace& V, const dolfinx::mesh::MeshTags& meshtags, std::int32_t slave_marker) { @@ -184,13 +183,12 @@ locate_slave_dofs(std::shared_ptr> V, // Find all dofs on slave facets // NOTE: Assumption that we are only working with vector spaces, which is // ordered as xyz,xyzgeometry - auto [V0, map] = V->sub({0})->collapse(); + auto [V0, map] = V.sub({0})->collapse(); std::array, 2> slave_dofs = dolfinx::fem::locate_dofs_topological( - *V->mesh()->topology_mutable(), - {*V->sub({0})->dofmap(), *V0.dofmap()}, edim, - std::span(slave_facets)); + *V.mesh()->topology_mutable(), {*V.sub({0})->dofmap(), *V0.dofmap()}, + edim, std::span(slave_facets)); return slave_dofs[0]; } @@ -347,27 +345,27 @@ namespace dolfinx_mpc /// collision template mpc_data create_contact_slip_condition( - std::shared_ptr> V, + const dolfinx::fem::FunctionSpace& V, const dolfinx::mesh::MeshTags& meshtags, std::int32_t slave_marker, std::int32_t master_marker, - std::shared_ptr> nh, const U eps2 = 1e-20) + const dolfinx::fem::Function& nh, const U eps2 = 1e-20) { dolfinx::common::Timer timer("~MPC: Create slip constraint"); - std::shared_ptr> mesh = V->mesh(); + std::shared_ptr> mesh = V.mesh(); MPI_Comm comm = mesh->comm(); int rank = -1; MPI_Comm_rank(comm, &rank); // Extract some const information from function-space const std::shared_ptr imap - = V->dofmap()->index_map; + = V.dofmap()->index_map; assert(mesh->topology() == meshtags.topology()); const int tdim = mesh->topology()->dim(); const int gdim = mesh->geometry().dim(); const int fdim = tdim - 1; - const int block_size = V->dofmap()->index_map_bs(); - std::int32_t size_local = V->dofmap()->index_map->size_local(); + const int block_size = V.dofmap()->index_map_bs(); + std::int32_t size_local = V.dofmap()->index_map->size_local(); mesh->topology_mutable()->create_connectivity(fdim, tdim); mesh->topology_mutable()->create_connectivity(tdim, tdim); @@ -408,7 +406,7 @@ mpc_data create_contact_slip_condition( // Note that this function casts the normal array from being potentially // complex to real valued std::vector dofs(block_size); - std::span normal_array = nh->x()->array(); + std::span normal_array = nh.x()->array(); const auto largest_normal_component = [&dofs, block_size, &normal_array, gdim](const std::int32_t block, std::span normal) @@ -449,19 +447,19 @@ mpc_data create_contact_slip_condition( // Create map from slave dof blocks to a cell containing them std::vector slave_cells = dolfinx_mpc::create_block_to_cell_map( - *mesh->topology(), *V->dofmap(), local_slave_blocks); + *mesh->topology(), *V.dofmap(), local_slave_blocks); std::vector slave_coordinates; std::array coord_shape; { std::tie(slave_coordinates, coord_shape) - = dolfinx_mpc::tabulate_dof_coordinates(*V, local_slave_blocks, + = dolfinx_mpc::tabulate_dof_coordinates(V, local_slave_blocks, slave_cells); std::vector local_cell_collisions = dolfinx_mpc::find_local_collisions(*mesh, bb_tree, slave_coordinates, eps2); auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( - *V, slave_coordinates, local_cell_collisions); + V, slave_coordinates, local_cell_collisions); assert(basis_shape.back() == 1); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> @@ -582,7 +580,7 @@ mpc_data create_contact_slip_condition( = dolfinx_mpc::find_local_collisions(*mesh, bb_tree, recv_coords, eps2); auto [recv_basis_values, shape] = dolfinx_mpc::evaluate_basis_functions( - *V, recv_coords, remote_cell_collisions); + V, recv_coords, remote_cell_collisions); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> basis_span(recv_basis_values.data(), shape[0], shape[1]); @@ -884,30 +882,30 @@ mpc_data create_contact_slip_condition( /// collision template mpc_data create_contact_inelastic_condition( - std::shared_ptr> V, + const dolfinx::fem::FunctionSpace& V, dolfinx::mesh::MeshTags meshtags, std::int32_t slave_marker, std::int32_t master_marker, const double eps2 = 1e-20) { dolfinx::common::Timer timer("~MPC: Inelastic condition"); - MPI_Comm comm = V->mesh()->comm(); + MPI_Comm comm = V.mesh()->comm(); int rank = -1; MPI_Comm_rank(comm, &rank); // Extract some const information from function-space const std::shared_ptr imap - = V->dofmap()->index_map; - const int tdim = V->mesh()->topology()->dim(); + = V.dofmap()->index_map; + const int tdim = V.mesh()->topology()->dim(); const int fdim = tdim - 1; - const int block_size = V->dofmap()->index_map_bs(); - std::int32_t size_local = V->dofmap()->index_map->size_local(); + const int block_size = V.dofmap()->index_map_bs(); + std::int32_t size_local = V.dofmap()->index_map->size_local(); // Create entity permutations needed in evaluate_basis_functions - V->mesh()->topology_mutable()->create_entity_permutations(); + V.mesh()->topology_mutable()->create_entity_permutations(); // Create connectivities needed for evaluate_basis_functions and // select_colliding cells - V->mesh()->topology_mutable()->create_connectivity(fdim, tdim); - V->mesh()->topology_mutable()->create_connectivity(tdim, tdim); + V.mesh()->topology_mutable()->create_connectivity(fdim, tdim); + V.mesh()->topology_mutable()->create_connectivity(tdim, tdim); std::vector slave_blocks = impl::locate_slave_dofs(V, meshtags, slave_marker); @@ -962,19 +960,19 @@ mpc_data create_contact_inelastic_condition( } // Create boundingboxtree for master surface - auto facet_to_cell = V->mesh()->topology()->connectivity(fdim, tdim); + auto facet_to_cell = V.mesh()->topology()->connectivity(fdim, tdim); assert(facet_to_cell); dolfinx::geometry::BoundingBoxTree bb_tree = impl::create_boundingbox_tree( - *V->mesh(), meshtags, master_marker, std::sqrt(eps2)); + *V.mesh(), meshtags, master_marker, std::sqrt(eps2)); // Tabulate slave block coordinates and find colliding cells std::vector slave_cells = dolfinx_mpc::create_block_to_cell_map( - *V->mesh()->topology(), *V->dofmap(), local_blocks); + *V.mesh()->topology(), *V.dofmap(), local_blocks); std::vector slave_coordinates; { std::array c_shape; std::tie(slave_coordinates, c_shape) - = dolfinx_mpc::tabulate_dof_coordinates(*V, local_blocks, + = dolfinx_mpc::tabulate_dof_coordinates(V, local_blocks, slave_cells); } @@ -987,10 +985,10 @@ mpc_data create_contact_inelastic_condition( std::vector collision_to_local; { std::vector colliding_cells - = dolfinx_mpc::find_local_collisions(*V->mesh(), bb_tree, + = dolfinx_mpc::find_local_collisions(*V.mesh(), bb_tree, slave_coordinates, eps2); auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( - *V, slave_coordinates, colliding_cells); + V, slave_coordinates, colliding_cells); assert(basis_shape.back() == 1); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> @@ -1003,7 +1001,7 @@ mpc_data create_contact_inelastic_condition( { if (const auto& cell = colliding_cells[i]; cell != -1) { - auto cell_blocks = V->dofmap()->cell_dofs(cell); + auto cell_blocks = V.dofmap()->cell_dofs(cell); l_master.reserve(cell_blocks.size()); coeff.reserve(cell_blocks.size()); assert(l_master.empty()); @@ -1172,10 +1170,10 @@ mpc_data create_contact_inelastic_condition( collision_block_offsets(indegree); { std::vector remote_cell_collisions - = dolfinx_mpc::find_local_collisions(*V->mesh(), bb_tree, - recv_coords, eps2); + = dolfinx_mpc::find_local_collisions(*V.mesh(), bb_tree, recv_coords, + eps2); auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( - *V, recv_coords, remote_cell_collisions); + V, recv_coords, remote_cell_collisions); assert(basis_shape.back() == 1); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> @@ -1196,7 +1194,7 @@ mpc_data create_contact_inelastic_condition( = std::vector(tdim, 0); if (const auto& cell = remote_cell_collisions[j]; cell != -1) { - auto cell_blocks = V->dofmap()->cell_dofs(cell); + auto cell_blocks = V.dofmap()->cell_dofs(cell); r_master.reserve(cell_blocks.size()); r_coeff.reserve(cell_blocks.size()); assert(r_master.empty()); diff --git a/cpp/MultiPointConstraint.h b/cpp/MultiPointConstraint.h index f1b5e8ec..3f54dca9 100644 --- a/cpp/MultiPointConstraint.h +++ b/cpp/MultiPointConstraint.h @@ -34,11 +34,11 @@ class MultiPointConstraint /// @param[in] offsets Offsets for masters /// @tparam The floating type of the mesh MultiPointConstraint(std::shared_ptr> V, - const std::vector& slaves, - const std::vector& masters, - const std::vector& coeffs, - const std::vector& owners, - const std::vector& offsets) + std::span slaves, + std::span masters, + std::span coeffs, + std::span owners, + std::span offsets) : _slaves(), _is_slave(), _cell_to_slaves_map(), _num_local_slaves(), _master_map(), _coeff_map(), _owner_map(), _mpc_constants(), _V() { @@ -64,7 +64,7 @@ class MultiPointConstraint _is_slave = std::move(_slave_data); // Create a map for cells owned by the process to the slaves - _cell_to_slaves_map = create_cell_to_dofs_map(V, slaves); + _cell_to_slaves_map = create_cell_to_dofs_map(*V, slaves); // Create adjacency list with all local dofs, where the slave dofs maps to // its masters @@ -117,11 +117,11 @@ class MultiPointConstraint // Create new function space with extended index map _V = std::make_shared>( - create_extended_functionspace(V, _master_data, _owner_data)); + create_extended_functionspace(*V, _master_data, _owner_data)); // Map global masters to local index in extended function space std::vector masters_local - = map_dofs_global_to_local(_V, _master_data); + = map_dofs_global_to_local(*_V, _master_data); _master_map = std::make_shared>( masters_local, masters_offsets); } diff --git a/cpp/assemble_matrix.cpp b/cpp/assemble_matrix.cpp index a336aba1..7e1d0f5d 100644 --- a/cpp/assemble_matrix.cpp +++ b/cpp/assemble_matrix.cpp @@ -561,11 +561,11 @@ void assemble_matrix_impl( std::function, const std::span, const std::int32_t, const int)> apply_dof_transformation - = element0->template get_dof_transformation_function(); + = element0->template get_pre_dof_transformation_function(); std::function, const std::span, const std::int32_t, const int)> apply_dof_transformation_to_transpose - = element1->template get_dof_transformation_to_transpose_function(); + = element1->template get_post_dof_transformation_function(); const bool needs_transformation_data = element0->needs_dof_transformations() diff --git a/cpp/assemble_vector.cpp b/cpp/assemble_vector.cpp index a762272b..5a5126ba 100644 --- a/cpp/assemble_vector.cpp +++ b/cpp/assemble_vector.cpp @@ -119,7 +119,8 @@ void _assemble_vector( const std::function&, const std::span&, std::int32_t, int)> - dof_transform = element->template get_dof_transformation_function(); + dof_transform + = element->template get_pre_dof_transformation_function(); const bool needs_transformation_data = element->needs_dof_transformations() or L.needs_facet_permutations(); std::span cell_info; diff --git a/cpp/lifting.h b/cpp/lifting.h index 767a661d..032f9911 100644 --- a/cpp/lifting.h +++ b/cpp/lifting.h @@ -203,12 +203,13 @@ void apply_lifting( const std::function&, const std::span&, std::int32_t, int)> - dof_transform = element0->template get_dof_transformation_function(); + dof_transform + = element0->template get_pre_dof_transformation_function(); const std::function&, const std::span&, std::int32_t, int)> dof_transform_to_transpose - = element1->template get_dof_transformation_to_transpose_function(); + = element1->template get_post_dof_transformation_function(); // Loop over cell integrals and lift bc if (a->num_integrals(dolfinx::fem::IntegralType::cell) > 0) diff --git a/cpp/mpc_helpers.h b/cpp/mpc_helpers.h index afbc29ba..8ee66850 100644 --- a/cpp/mpc_helpers.h +++ b/cpp/mpc_helpers.h @@ -18,11 +18,11 @@ namespace dolfinx_mpc /// process) in the cell template std::shared_ptr> -create_cell_to_dofs_map(std::shared_ptr> V, +create_cell_to_dofs_map(const dolfinx::fem::FunctionSpace& V, std::span dofs) { - const auto& mesh = *(V->mesh()); - const dolfinx::fem::DofMap& dofmap = *(V->dofmap()); + const auto& mesh = *(V.mesh()); + const dolfinx::fem::DofMap& dofmap = *(V.dofmap()); const int tdim = mesh.topology()->dim(); const int num_cells = mesh.topology()->index_map(tdim)->size_local(); @@ -99,15 +99,15 @@ create_cell_to_dofs_map(std::shared_ptr> V, /// @tparam The floating type of the mesh /// @returns List of local dofs template -std::vector map_dofs_global_to_local( - std::shared_ptr> V, - const std::vector& global_dofs) +std::vector +map_dofs_global_to_local(const dolfinx::fem::FunctionSpace& V, + const std::vector& global_dofs) { const std::size_t num_dofs = global_dofs.size(); - const std::int32_t& block_size = V->dofmap()->index_map_bs(); + const std::int32_t& block_size = V.dofmap()->index_map_bs(); const std::shared_ptr imap - = V->dofmap()->index_map; + = V.dofmap()->index_map; std::vector global_blocks; global_blocks.reserve(num_dofs); @@ -136,19 +136,19 @@ std::vector map_dofs_global_to_local( /// @param[in] owners The owners of the master degrees of freedom /// @tparam The floating type of the mesh template -dolfinx::fem::FunctionSpace create_extended_functionspace( - std::shared_ptr> V, - const std::vector& global_dofs, - const std::vector& owners) +dolfinx::fem::FunctionSpace +create_extended_functionspace(const dolfinx::fem::FunctionSpace& V, + const std::vector& global_dofs, + const std::vector& owners) { dolfinx::common::Timer timer( "~MPC: Create new index map with additional ghosts"); - MPI_Comm comm = V->mesh()->comm(); - const dolfinx::fem::DofMap& old_dofmap = *(V->dofmap()); + MPI_Comm comm = V.mesh()->comm(); + const dolfinx::fem::DofMap& old_dofmap = *(V.dofmap()); std::shared_ptr old_index_map = old_dofmap.index_map; - const std::int32_t& block_size = V->dofmap()->index_map_bs(); + const std::int32_t& block_size = V.dofmap()->index_map_bs(); // Compute local master block index. const std::size_t num_dofs = global_dofs.size(); @@ -168,7 +168,7 @@ dolfinx::fem::FunctionSpace create_extended_functionspace( else { // Map global master blocks to local blocks - V->dofmap()->index_map->global_to_local(global_blocks, local_blocks); + V.dofmap()->index_map->global_to_local(global_blocks, local_blocks); // Check which local masters that are not on the process already std::vector new_ghosts; @@ -220,14 +220,14 @@ dolfinx::fem::FunctionSpace create_extended_functionspace( for (std::size_t j = 0; j < dofmap_adj.extent(1); ++j) flattened_dofmap.push_back(dofmap_adj(i, j)); - auto element = V->element(); + auto element = V.element(); // Create the new dofmap based on the extended index map auto new_dofmap = std::make_shared( old_dofmap.element_dof_layout(), new_index_map, old_dofmap.bs(), std::move(flattened_dofmap), old_dofmap.bs()); - return dolfinx::fem::FunctionSpace(V->mesh(), element, new_dofmap); + return dolfinx::fem::FunctionSpace(V.mesh(), element, new_dofmap); } } // namespace dolfinx_mpc \ No newline at end of file diff --git a/cpp/utils.h b/cpp/utils.h index a8cb0ab7..8b6da82a 100644 --- a/cpp/utils.h +++ b/cpp/utils.h @@ -1150,7 +1150,7 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, = element->basix_element().template map_fn(); auto apply_dof_transformation - = element->template get_dof_transformation_function(); + = element->template get_pre_dof_transformation_function(); mdspan3_t full_basis(output_basis.data(), reference_shape); for (std::size_t p = 0; p < cells.size(); ++p) @@ -1276,7 +1276,7 @@ std::pair, std::array> tabulate_dof_coordinates( } const auto apply_dof_transformation - = element->template get_dof_transformation_function(); + = element->template get_pre_dof_transformation_function(); const std::array bsize = cmaps[0].tabulate_shape(0, X_shape[0]); diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 89df3265..34efc54d 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -2,25 +2,82 @@ cmake_minimum_required(VERSION 3.19) PROJECT(dolfinx_mpc_pybind11) -find_package(pybind11 REQUIRED CONFIG HINTS ${PYBIND11_DIR} ${PYBIND11_ROOT} - $ENV{PYBIND11_DIR} $ENV{PYBIND11_ROOT}) - -find_package(DOLFINX_MPC REQUIRED) -MESSAGE(STATUS "Found dolfinx_mpc ${DOLFINX_MPC_VERSION}") -MESSAGE(STATUS "Found dolfinx ${DOLFINX_VERSION}") - -# Set C++ standard before finding pybind11 +# Set C++ standard set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) -# Create the binding library -pybind11_add_module(cpp SHARED - dolfinx_mpc/mpc.cpp - dolfinx_mpc/dolfinx_mpc.cpp) +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE) + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo") +endif() + +find_package(Python COMPONENTS Interpreter Development REQUIRED) + +# Detect the installed nanobind package and import it into CMake +execute_process( + COMMAND "${Python_EXECUTABLE}" -m nanobind --cmake_dir + OUTPUT_STRIP_TRAILING_WHITESPACE OUTPUT_VARIABLE NB_DIR) +list(APPEND CMAKE_PREFIX_PATH "${NB_DIR}") +find_package(nanobind CONFIG REQUIRED) + +execute_process( + COMMAND + ${Python3_EXECUTABLE} -c "import os, sys, basix; sys.stdout.write(os.path.dirname(basix.__file__))" + OUTPUT_VARIABLE BASIX_PY_DIR + RESULT_VARIABLE BASIX_PY_COMMAND_RESULT + ERROR_VARIABLE BASIX_ERROR_OUT OUTPUT_STRIP_TRAILING_WHITESPACE +) +find_package(Basix REQUIRED CONFIG HINTS ${BASIX_PY_DIR}) +if (Basix_FOUND) + message(STATUS "Found Basix at ${Basix_DIR}") +endif() +find_package(DOLFINX REQUIRED CONFIG) +if (DOLFINX_FOUND) + message(STATUS "Found DOLFINx at ${DOLFINX_DIR}") +endif() + + +# Add DOLFINx_mpc libraries +find_package(DOLFINX_MPC REQUIRED CONFIG) +if (DOLFINX_MPC_FOUND) + message(STATUS "Found DOLFINx_MPC at ${DOLFINX_MPC_DIR}") +endif() + +# Create the binding library nanobind handles its own calls to +# target_link_libraries +nanobind_add_module( + cpp + NOMINSIZE + MODULE + dolfinx_mpc/dolfinx_mpc.cpp + dolfinx_mpc/mpc.cpp +) +# Add strict compiler flags +# include(CheckCXXCompilerFlag) +# check_cxx_compiler_flag("-Wall -Werror -pedantic" HAVE_PEDANTIC) + +# if(HAVE_PEDANTIC) +# # target_compile_options(cpp PRIVATE -Wall;-Werror;-pedantic) +# endif() +target_link_libraries(cpp PRIVATE dolfinx_mpc) + +# Check for petsc4py +execute_process( + COMMAND ${Python_EXECUTABLE} -c "import petsc4py; print(petsc4py.get_include())" + OUTPUT_VARIABLE PETSC4PY_INCLUDE_DIR + RESULT_VARIABLE PETSC4PY_COMMAND_RESULT + ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE +) + +if(NOT PETSC4PY_COMMAND_RESULT) + message(STATUS "Found petsc4py include directory at ${PETSC4PY_INCLUDE_DIR}") + target_include_directories(cpp PRIVATE ${PETSC4PY_INCLUDE_DIR}) +else() + message(FATAL_ERROR "petsc4py could not be found.") +endif() + -target_link_libraries(cpp PRIVATE pybind11::module) -target_link_libraries(cpp PUBLIC dolfinx_mpc) # Get python include-dirs execute_process( @@ -33,20 +90,7 @@ if (DOLFINX_PY_DIR) target_include_directories(cpp PRIVATE ${DOLFINX_PY_DIR}) endif() -# Add strict compiler flags -include(CheckCXXCompilerFlag) -CHECK_CXX_COMPILER_FLAG("-Wall -Werror -pedantic" HAVE_PEDANTIC) -if (HAVE_PEDANTIC) - target_compile_options(cpp PRIVATE -Wall;-Werror;-pedantic) -endif() +set_target_properties(cpp PROPERTIES INSTALL_RPATH_USE_LINK_PATH TRUE) -# Find petsc4py through python -execute_process( - COMMAND ${PYTHON_EXECUTABLE} -c "import petsc4py; print(petsc4py.get_include())" - OUTPUT_VARIABLE PETSC4PY_INCLUDE_DIR - RESULT_VARIABLE PETSC4PY_NOT_FOUND - ERROR_QUIET - OUTPUT_STRIP_TRAILING_WHITESPACE - ) -target_include_directories(cpp PRIVATE ${PETSC4PY_INCLUDE_DIR}) +install(TARGETS cpp DESTINATION dolfinx_mpc) diff --git a/python/README.md b/python/README.md index f7f747a6..6fa3bb3b 100644 --- a/python/README.md +++ b/python/README.md @@ -1 +1,7 @@ -# Python interface for DolfinX-MPC \ No newline at end of file +# Python interface for DolfinX-MPC + +Can be installed with + +```bash +python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation . +``` \ No newline at end of file diff --git a/python/dolfinx_mpc/dictcondition.py b/python/dolfinx_mpc/dictcondition.py index 8412d80e..98b777d2 100644 --- a/python/dolfinx_mpc/dictcondition.py +++ b/python/dolfinx_mpc/dictcondition.py @@ -117,7 +117,7 @@ def create_dictionary_constraint(V: fem.functionspace, slave_master_dict: if len(master_dofs) == 1: master_block = master_dofs[0] // bs master_rem = master_dofs % bs - glob_master = index_map.local_to_global([master_block])[0] + glob_master = index_map.local_to_global(np.asarray([master_block], dtype=np.int32))[0] if slave_status == -1: if i in non_local_entities.keys(): non_local_entities[i]["masters"].append(glob_master * bs + master_rem) diff --git a/python/dolfinx_mpc/dolfinx_mpc.cpp b/python/dolfinx_mpc/dolfinx_mpc.cpp index 52165e0b..4d58e68f 100644 --- a/python/dolfinx_mpc/dolfinx_mpc.cpp +++ b/python/dolfinx_mpc/dolfinx_mpc.cpp @@ -4,23 +4,26 @@ // // SPDX-License-Identifier: MIT -#include -#include +#include -namespace py = pybind11; +namespace nb = nanobind; namespace dolfinx_mpc_wrappers { -void mpc(py::module& m); +void mpc(nb::module_& m); } // namespace dolfinx_mpc_wrappers -PYBIND11_MODULE(cpp, m) +NB_MODULE(cpp, m) { // Create module for C++ wrappers m.doc() = "DOLFINX MultiPointConstraint Python interface"; m.attr("__version__") = DOLFINX_MPC_VERSION; +#ifdef NDEBUG + nanobind::set_leak_warnings(false); +#endif + // Create mpc submodule [mpc] - py::module mpc = m.def_submodule("mpc", "General module"); + nb::module_ mpc = m.def_submodule("mpc", "General module"); dolfinx_mpc_wrappers::mpc(mpc); } diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index 8ffb3494..58af9fd7 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -1,8 +1,8 @@ -// Copyright (C) 2020 Jørgen S. Dokken -// -// This file is part of DOLFINX-MPC -// -// SPDX-License-Identifier: MIT +// // Copyright (C) 2020 Jørgen S. Dokken +// // +// // This file is part of DOLFINX-MPC +// // +// // SPDX-License-Identifier: MIT #include #include @@ -23,13 +23,20 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include -#include -#include -#include - -namespace py = pybind11; +namespace nb = nanobind; +using namespace nb::literals; namespace @@ -37,23 +44,31 @@ namespace // Templating over mesh resolution template -void declare_mpc(py::module& m, std::string type) +void declare_mpc(nb::module_& m, std::string type) { - std::string pyclass_name = "MultiPointConstraint_" + type; + std::string nbclass_name = "MultiPointConstraint_" + type; // dolfinx_mpc::MultiPointConstraint - py::class_, - std::shared_ptr>> - multipointconstraint(m, pyclass_name.c_str(), - "Object for representing contact " - "(non-penetrating) conditions"); - multipointconstraint - .def(py::init>, - std::vector, std::vector, - std::vector, std::vector, - std::vector>()) - .def_property_readonly("masters", - &dolfinx_mpc::MultiPointConstraint::masters) + nb::class_>( + m, nbclass_name.c_str(), + "Object for representing contact (non-penetrating) conditions") + .def("__init__", + [](dolfinx_mpc::MultiPointConstraint* mpc, + std::shared_ptr> V, + nb::ndarray>& slaves, + nb::ndarray>& masters, + nb::ndarray>& coeffs, + nb::ndarray>& owners, + nb::ndarray>& offsets) + { + new (mpc) dolfinx_mpc::MultiPointConstraint( + V, std::span(slaves.data(), slaves.size()), + std::span(masters.data(), masters.size()), + std::span(coeffs.data(), coeffs.size()), + std::span(owners.data(), owners.size()), + std::span(offsets.data(), offsets.size())); + }) + .def_prop_ro("masters", &dolfinx_mpc::MultiPointConstraint::masters) .def("coefficients", [](dolfinx_mpc::MultiPointConstraint& self) { @@ -62,66 +77,61 @@ void declare_mpc(py::module& m, std::string type) const std::vector& offsets = adj->offsets(); const std::vector& data = adj->array(); - return std::pair, py::array_t>( - py::array_t(data.size(), data.data(), py::cast(self)), - py::array_t(offsets.size(), offsets.data(), - py::cast(self))); + return std::make_pair( + nb::ndarray>(data.data(), + {data.size()}), + nb::ndarray>( + offsets.data(), {offsets.size()})); }) - .def_property_readonly( - "constants", - [](dolfinx_mpc::MultiPointConstraint& self) - { - const std::vector& consts = self.constant_values(); - return py::array_t(consts.size(), consts.data(), py::cast(self)); - }) - .def_property_readonly("owners", - &dolfinx_mpc::MultiPointConstraint::owners) - .def_property_readonly( + .def_prop_ro("constants", + [](dolfinx_mpc::MultiPointConstraint& self) + { + const std::vector& consts = self.constant_values(); + return nb::ndarray>( + consts.data(), {consts.size()}); + }) + .def_prop_ro("owners", &dolfinx_mpc::MultiPointConstraint::owners) + .def_prop_ro( "slaves", [](dolfinx_mpc::MultiPointConstraint& self) { const std::vector& slaves = self.slaves(); - return py::array_t(slaves.size(), slaves.data(), - py::cast(self)); + return nb::ndarray>( + slaves.data(), {slaves.size()}); }) - .def_property_readonly( + .def_prop_ro( "is_slave", [](dolfinx_mpc::MultiPointConstraint& self) { std::span slaves = self.is_slave(); - return py::array_t(slaves.size(), slaves.data(), - py::cast(self)); + return nb::ndarray>( + slaves.data(), {slaves.size()}); }) - .def_property_readonly( - "cell_to_slaves", - &dolfinx_mpc::MultiPointConstraint::cell_to_slaves) - .def_property_readonly( - "num_local_slaves", - &dolfinx_mpc::MultiPointConstraint::num_local_slaves) - .def_property_readonly( - "function_space", - &dolfinx_mpc::MultiPointConstraint::function_space) - .def_property_readonly("owners", - &dolfinx_mpc::MultiPointConstraint::owners) + .def_prop_ro("cell_to_slaves", + &dolfinx_mpc::MultiPointConstraint::cell_to_slaves) + .def_prop_ro("num_local_slaves", + &dolfinx_mpc::MultiPointConstraint::num_local_slaves) + .def_prop_ro("function_space", + &dolfinx_mpc::MultiPointConstraint::function_space) + .def_prop_ro("owners", &dolfinx_mpc::MultiPointConstraint::owners) .def( "backsubstitution", [](dolfinx_mpc::MultiPointConstraint& self, - py::array_t u) - { self.backsubstitution(std::span(u.mutable_data(), u.size())); }, - py::arg("u"), "Backsubstitute slave values into vector") + nb::ndarray, nb::c_contig> u) + { self.backsubstitution(std::span(u.data(), u.size())); }, + "u"_a, "Backsubstitute slave values into vector") .def( "homogenize", [](dolfinx_mpc::MultiPointConstraint& self, - py::array_t u) - { self.homogenize(std::span(u.mutable_data(), u.size())); }, - py::arg("u"), "Homogenize (set to zero) values at slave DoF indices"); + nb::ndarray, nb::c_contig> u) + { self.homogenize(std::span(u.data(), u.size())); }, + "u"_a, "Homogenize (set to zero) values at slave DoF indices"); // .def("ghost_masters", &dolfinx_mpc::mpc_data::ghost_masters); } - template -void declare_functions(py::module& m) +void declare_functions(nb::module_& m) { m.def("compute_shared_indices", &dolfinx_mpc::compute_shared_indices); @@ -136,9 +146,10 @@ void declare_functions(py::module& m) m.def( "create_periodic_constraint_geometrical", [](const std::shared_ptr> V, - const std::function(const py::array_t&)>& - indicator, - const std::function(const py::array_t&)>& relation, + const std::function( + nb::ndarray, nb::numpy>&)>& indicator, + const std::function, nb::numpy>( + nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, T scale, bool collapse) @@ -152,9 +163,9 @@ void declare_functions(py::module& m) x) -> std::vector { assert(x.size() % 3 == 0); - const std::array shape = {3, x.size() / 3}; - const py::array_t _x(shape, x.data_handle(), py::none()); - py::array_t m = indicator(_x); + nb::ndarray, nb::numpy> x_view( + x.data_handle(), {3, x.size() / 3}); + auto m = indicator(x_view); std::vector s(m.data(), m.data() + m.size()); return s; }; @@ -162,15 +173,16 @@ void declare_functions(py::module& m) auto _relation = [&relation](std::span x) -> std::vector { assert(x.size() % 3 == 0); - const std::array shape = {3, x.size() / 3}; - const py::array_t _x(shape, x.data(), py::none()); - py::array_t v = relation(_x); + nb::ndarray, nb::numpy> x_view( + x.data(), {3, x.size() / 3}); + auto v = relation(x_view); std::vector output(v.data(), v.data() + v.size()); return output; }; return dolfinx_mpc::create_periodic_condition_geometrical( V, _indicator, _relation, bcs, scale, collapse); - }); + }, + "V"_a, "indicator"_a, "relation"_a, "bcs"_a, "scall"_a, "collapse"_a); m.def( "create_periodic_constraint_topological", @@ -178,79 +190,81 @@ void declare_functions(py::module& m) const std::shared_ptr>& meshtags, const int dim, - const std::function(const py::array_t&)>& relation, + const std::function, nb::numpy>( + nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, T scale, bool collapse) { auto _relation = [&relation](std::span x) -> std::vector { - const std::array shape = {3, x.size() / 3}; - py::array_t _x(shape, x.data(), py::none()); - py::array_t v = relation(_x); + nb::ndarray, nb::numpy> x_view( + x.data(), {3, x.size() / 3}); + auto v = relation(x_view); std::vector output(v.data(), v.data() + v.size()); return output; }; return dolfinx_mpc::create_periodic_condition_topological( V, meshtags, dim, _relation, bcs, scale, collapse); - }); + }, + "V"_a, "meshtags"_a, "dim"_a, "relation"_a, "bcs"_a, "scall"_a, + "collapse"_a); } template -void declare_mpc_data(py::module& m, std::string type) +void declare_mpc_data(nb::module_& m, std::string type) { - std::string pyclass_name = "mpc_data_" + type; - py::class_, - std::shared_ptr>> - mpc_data(m, pyclass_name.c_str(), "Object with data arrays for mpc"); - mpc_data - .def_property_readonly( + std::string nbclass_name = "mpc_data_" + type; + nb::class_>(m, nbclass_name.c_str(), + "Object with data arrays for mpc") + .def_prop_ro( "slaves", [](dolfinx_mpc::mpc_data& self) { const std::vector& slaves = self.slaves; - return py::array_t(slaves.size(), slaves.data(), - py::cast(self)); + return nb::ndarray>( + slaves.data(), {slaves.size()}); }) - .def_property_readonly( + .def_prop_ro( "masters", [](dolfinx_mpc::mpc_data& self) { const std::vector& masters = self.masters; - return py::array_t(masters.size(), masters.data(), - py::cast(self)); - }) - .def_property_readonly( - "coeffs", - [](dolfinx_mpc::mpc_data& self) - { - const std::vector& coeffs = self.coeffs; - return py::array_t(coeffs.size(), coeffs.data(), py::cast(self)); + return nb::ndarray>( + masters.data(), {masters.size()}); }) - .def_property_readonly( + .def_prop_ro("coeffs", + [](dolfinx_mpc::mpc_data& self) + { + const std::vector& coeffs = self.coeffs; + return nb::ndarray>( + coeffs.data(), {coeffs.size()}); + }) + .def_prop_ro( "owners", [](dolfinx_mpc::mpc_data& self) { const std::vector& owners = self.owners; - return py::array_t(owners.size(), owners.data(), - py::cast(self)); + return nb::ndarray>( + owners.data(), {owners.size()}); }) - .def_property_readonly( + .def_prop_ro( "offsets", [](dolfinx_mpc::mpc_data& self) { const std::vector& offsets = self.offsets; - return py::array_t(offsets.size(), offsets.data(), - py::cast(self)); + return nb::ndarray>( + offsets.data(), {offsets.size()}); }); } template -void declare_petsc_functions(py::module& m) +void declare_petsc_functions(nb::module_& m) { + import_petsc4py(); m.def("create_normal_approximation", [](std::shared_ptr> V, std::int32_t dim, - const py::array_t& entities) + const nb::ndarray, nb::c_contig>& entities) { return dolfinx_mpc::create_normal_approximation( V, dim, @@ -274,34 +288,32 @@ void declare_petsc_functions(py::module& m) }); m.def( "assemble_vector", - [](py::array_t b, const dolfinx::fem::Form& L, + [](nb::ndarray, nb::c_contig> b, + const dolfinx::fem::Form& L, const std::shared_ptr>& mpc) - { - dolfinx_mpc::assemble_vector(std::span(b.mutable_data(), b.size()), L, - mpc); - }, - py::arg("b"), py::arg("L"), py::arg("mpc"), - "Assemble linear form into an existing vector"); + { dolfinx_mpc::assemble_vector(std::span(b.data(), b.size()), L, mpc); }, + "b"_a, "L"_a, "mpc"_a, "Assemble linear form into an existing vector"); m.def( "apply_lifting", - [](py::array_t b, + [](nb::ndarray, nb::c_contig> b, std::vector>>& a, const std::vector>>>& bcs1, - const std::vector>& x0, U scale, + const std::vector, nb::c_contig>>& x0, + U scale, std::shared_ptr>& mpc) { std::vector> _x0; for (const auto& x : x0) _x0.emplace_back(x.data(), x.size()); - dolfinx_mpc::apply_lifting(std::span(b.mutable_data(), b.size()), a, - bcs1, _x0, scale, mpc); + dolfinx_mpc::apply_lifting(std::span(b.data(), b.size()), a, bcs1, _x0, + scale, mpc); }, - py::arg("b"), py::arg("a"), py::arg("bcs"), py::arg("x0"), - py::arg("scale"), py::arg("mpc"), + nb::arg("b"), nb::arg("a"), nb::arg("bcs"), nb::arg("x0"), + nb::arg("scale"), nb::arg("mpc"), "Assemble apply lifting from form a on vector b"); m.def( @@ -312,10 +324,10 @@ void declare_petsc_functions(py::module& m) auto A = dolfinx_mpc::create_matrix(a, mpc); Mat _A = A.mat(); PetscObjectReference((PetscObject)_A); + return _A; }, - py::return_value_policy::take_ownership, - "Create a PETSc Mat for bilinear form."); + nb::rv_policy::take_ownership, "Create a PETSc Mat for bilinear form."); m.def( "create_matrix", [](const dolfinx::fem::Form& a, @@ -327,8 +339,7 @@ void declare_petsc_functions(py::module& m) PetscObjectReference((PetscObject)_A); return _A; }, - py::return_value_policy::take_ownership, - "Create a PETSc Mat for bilinear form."); + nb::rv_policy::take_ownership, "Create a PETSc Mat for bilinear form."); } } // namespace @@ -336,7 +347,7 @@ void declare_petsc_functions(py::module& m) namespace dolfinx_mpc_wrappers { -void mpc(py::module& m) +void mpc(nb::module_& m) { declare_mpc(m, "float"); diff --git a/python/dolfinx_mpc/utils/test.py b/python/dolfinx_mpc/utils/test.py index 699d98b7..423722e5 100644 --- a/python/dolfinx_mpc/utils/test.py +++ b/python/dolfinx_mpc/utils/test.py @@ -49,7 +49,8 @@ def _gather_slaves_global(constraint): if num_local_slaves > 0: slave_blocks = _slaves[:num_local_slaves] // block_size slave_rems = _slaves[:num_local_slaves] % block_size - glob_slaves = imap.local_to_global(slave_blocks) * block_size + slave_rems + glob_slaves = np.asarray( + imap.local_to_global(slave_blocks), dtype=np.int64) * block_size + slave_rems else: glob_slaves = np.array([], dtype=np.int64) @@ -105,7 +106,8 @@ def gather_transformation_matrix(constraint, root=0): if num_local_slaves > 0: local_blocks = slaves // block_size local_rems = slaves % block_size - glob_slaves = imap.local_to_global(local_blocks) * block_size + local_rems + glob_slaves = np.asarray(imap.local_to_global(local_blocks), + dtype=np.int64) * block_size + local_rems else: glob_slaves = np.array([], dtype=np.int64) all_slaves = np.hstack(MPI.COMM_WORLD.allgather(glob_slaves)) @@ -119,8 +121,9 @@ def gather_transformation_matrix(constraint, root=0): # Add local contributions to K from local slaves for slave, global_slave in zip(slaves, glob_slaves): - masters_index = (imap.local_to_global(master_blocks[offsets[slave]: offsets[slave + 1]]) - * block_size + master_rems[offsets[slave]: offsets[slave + 1]]) + masters_index = (np.asarray(imap.local_to_global( + master_blocks[offsets[slave]: offsets[slave + 1]]), dtype=np.int64) + * block_size + master_rems[offsets[slave]: offsets[slave + 1]]) coeffs_index = coeffs[offsets[slave]: offsets[slave + 1]] # If we have a simply equality constraint (dirichletbc) if len(masters_index) > 0: diff --git a/python/pyproject.toml b/python/pyproject.toml new file mode 100644 index 00000000..6e72649f --- /dev/null +++ b/python/pyproject.toml @@ -0,0 +1,53 @@ +# The DOLFINx Python interface must be built without build isolation (PEP517) +# due to its runtime and build time dependency on system built petsc4py and +# mpi4py. +# pip install -r build-requirements.txt +[build-system] +requires = [ + "scikit-build-core[pyproject]", + "nanobind>=1.8.0", + "petsc4py", + "mpi4py", +] +build-backend = "scikit_build_core.build" + +[project] +name = "dolfinx_mpc" +version = "0.7.0.dev0" +description = "DOLFINx_MPC Python interface" +readme = "README.md" +requires-python = ">=3.8.0" +license = { file = "../LICENSE" } +authors = [ + { email = "dokken@simula.no", name = "Jørgen S. Dokken" }, +] +packages = ["dolfinx_mpc"] +dependencies = [ + "numpy>=1.21", + "cffi", + "petsc4py", + "mpi4py", + "fenics-dolfinx>=0.8.0.dev0,<0.9.0", +] + +[project.optional-dependencies] +docs = ['jupyter-book', 'jupytext'] +lint = ["flake8", "mypy"] +optional = ["numba"] +test = ["pytest", "coverage"] +all = ["dolfinx_mpc[docs]", +"dolfinx_mpc[optional]", +"dolfinx_mpc[lint]", +"dolfinx_mpc[test]" +] + +[tool.scikit-build] +wheel.packages = ["dolfinx_mpc"] +sdist.exclude = ["*.cpp"] +cmake.build-type = "Release" + +[tool.pytest] +junit_family = "xunit2" + +[tool.mypy] +ignore_missing_imports = true diff --git a/python/setup.cfg b/python/setup.cfg deleted file mode 100644 index eec6b8ba..00000000 --- a/python/setup.cfg +++ /dev/null @@ -1,10 +0,0 @@ -[flake8] -max-line-length = 120 -ignore = W503 - -[tool:pytest] -junit_family = xunit2 - - -[mypy] -ignore_missing_imports = True \ No newline at end of file diff --git a/python/setup.py b/python/setup.py deleted file mode 100644 index 718a79b8..00000000 --- a/python/setup.py +++ /dev/null @@ -1,94 +0,0 @@ -import os -import shlex -import platform -import re -import subprocess -import sys -from itertools import chain -from distutils.version import LooseVersion - -from setuptools import Extension, setup -from setuptools.command.build_ext import build_ext - -if sys.version_info < (3, 8): - print("Python 3.8 or higher required, please upgrade.") - sys.exit(1) - -VERSION = "0.7.0.dev0" - -REQUIREMENTS = ["numpy>=1.21", "fenics-dolfinx>0.7.0"] - -extras = { - 'docs': ['jupyter-book', 'jupytext'], - "test": ["pytest", "coverage"] -} - -# 'all' includes all of the above -extras['all'] = list(chain(*extras.values())) - - -class CMakeExtension(Extension): - def __init__(self, name, sourcedir=''): - Extension.__init__(self, name, sources=[]) - self.sourcedir = os.path.abspath(sourcedir) - - -class CMakeBuild(build_ext): - def run(self): - try: - out = subprocess.check_output(['cmake', '--version']) - except OSError: - raise RuntimeError("CMake must be installed to build the" - + "following extensions:" - + ", ".join(e.name for e in self.extensions)) - - if platform.system() == "Windows": - cmake_version = LooseVersion(re.search(r'version\s*([\d.]+)', - out.decode()).group(1)) - if cmake_version < '3.1.0': - raise RuntimeError("CMake >= 3.1.0 is required on Windows") - - for ext in self.extensions: - self.build_extension(ext) - - def build_extension(self, ext): - extdir = os.path.abspath(os.path.dirname( - self.get_ext_fullpath(ext.name))) - cmake_args = shlex.split(os.environ.get("CMAKE_ARGS", "")) - cmake_args += ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, - '-DPYTHON_EXECUTABLE=' + sys.executable] - - cfg = 'Debug' if self.debug else 'Release' - build_args = ['--config', cfg] - - env = os.environ.copy() - cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] - # default to 3 build threads - if "CMAKE_BUILD_PARALLEL_LEVEL" not in env: - env["CMAKE_BUILD_PARALLEL_LEVEL"] = "3" - - import pybind11 - env['pybind11_DIR'] = pybind11.get_cmake_dir() - env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format( - env.get('CXXFLAGS', ''), self.distribution.get_version()) - if not os.path.exists(self.build_temp): - os.makedirs(self.build_temp) - subprocess.check_call(['cmake', ext.sourcedir] + cmake_args, - cwd=self.build_temp, env=env) - subprocess.check_call(['cmake', '--build', '.'] + build_args, - cwd=self.build_temp, env=env) - - -setup(name='dolfinx-mpc', - version=VERSION, - author='Jørgen S. Dokken', - description='Python interface for multipointconstraints in dolfinx', - long_description='', - packages=["dolfinx_mpc", "dolfinx_mpc.utils", "dolfinx_mpc.numba"], - ext_modules=[CMakeExtension('dolfinx_mpc.cpp')], - package_data={'dolfinx_mpc.wrappers': ['*.h'], 'dolfinx_mpc': ['py.typed'], - 'dolfinx_mpc.numba': ['py.typed'], 'dolfinx_mpc.utils': ['py.typed']}, - cmdclass=dict(build_ext=CMakeBuild), - install_requires=REQUIREMENTS, - zip_safe=False, - extras_require=extras) diff --git a/python/tests/test_integration_domains.py b/python/tests/test_integration_domains.py index cafab4d8..cf6b2216 100644 --- a/python/tests/test_integration_domains.py +++ b/python/tests/test_integration_domains.py @@ -35,11 +35,12 @@ def left_side(x): tdim = mesh.topology.dim num_cells = mesh.topology.index_map(tdim).size_local - cell_midpoints = compute_midpoints(mesh, tdim, range(num_cells)) - values = np.ones(num_cells, dtype=np.intc) + cells = np.arange(num_cells, dtype=np.int32) + cell_midpoints = compute_midpoints(mesh, tdim, cells) + values = np.ones_like(cells) # All cells on right side marked one, all other with 1 values += left_side(cell_midpoints.T) - ct = meshtags(mesh, mesh.topology.dim, np.arange(num_cells, dtype=np.int32), values) + ct = meshtags(mesh, mesh.topology.dim, cells, values) # Solve Problem without MPC for reference u = ufl.TrialFunction(V) From 8438260fc5c4d7b68ef3217033619538286b4667 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 21:54:01 +0000 Subject: [PATCH 02/30] Fixes of implicit conversion of integers to garbage --- python/dolfinx_mpc/mpc.cpp | 11 +++++------ python/tests/test_linear_problem.py | 2 +- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index 58af9fd7..56881951 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -142,11 +142,10 @@ void declare_functions(nb::module_& m) m.def("create_slip_condition", &dolfinx_mpc::create_slip_condition); m.def("create_contact_inelastic_condition", &dolfinx_mpc::create_contact_inelastic_condition); - m.def( "create_periodic_constraint_geometrical", [](const std::shared_ptr> V, - const std::function( + const std::function, nb::c_contig>( nb::ndarray, nb::numpy>&)>& indicator, const std::function, nb::numpy>( nb::ndarray, nb::numpy>&)>& relation, @@ -182,8 +181,8 @@ void declare_functions(nb::module_& m) return dolfinx_mpc::create_periodic_condition_geometrical( V, _indicator, _relation, bcs, scale, collapse); }, - "V"_a, "indicator"_a, "relation"_a, "bcs"_a, "scall"_a, "collapse"_a); - + "V"_a, "indicator"_a, "relation"_a, "bcs"_a, nb::arg("scale").noconvert(), + nb::arg("collapse").noconvert()); m.def( "create_periodic_constraint_topological", [](const std::shared_ptr>& V, @@ -207,8 +206,8 @@ void declare_functions(nb::module_& m) return dolfinx_mpc::create_periodic_condition_topological( V, meshtags, dim, _relation, bcs, scale, collapse); }, - "V"_a, "meshtags"_a, "dim"_a, "relation"_a, "bcs"_a, "scall"_a, - "collapse"_a); + "V"_a, "meshtags"_a, "dim"_a, "relation"_a, "bcs"_a, + nb::arg("scale").noconvert(), nb::arg("collapse").noconvert()); } template diff --git a/python/tests/test_linear_problem.py b/python/tests/test_linear_problem.py index 2942fb16..f1757561 100644 --- a/python/tests/test_linear_problem.py +++ b/python/tests/test_linear_problem.py @@ -56,7 +56,7 @@ def PeriodicBoundary(x): mt = meshtags(mesh, mesh.topology.dim - 1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int32)) mpc = dolfinx_mpc.MultiPointConstraint(V) - mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, [], 1) + mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, [], 1.) mpc.finalize() if u_from_mpc: From 58d7a03cb3a15eac0e2ff7eea4b7550ff930fa28 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 22:18:10 +0000 Subject: [PATCH 03/30] Various conversion for nanobind --- python/dolfinx_mpc/assemble_vector.py | 4 +++- python/dolfinx_mpc/mpc.cpp | 2 +- python/dolfinx_mpc/multipointconstraint.py | 4 ++++ 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index d30ac265..1ed89601 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -14,7 +14,7 @@ import ufl from dolfinx.common import Timer from petsc4py import PETSc as _PETSc - +import numpy import dolfinx_mpc.cpp from .multipointconstraint import MultiPointConstraint @@ -34,6 +34,8 @@ def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.Diri x0: List of vectors scale: Scaling for lifting """ + if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types + scale = scale.item() t = Timer("~MPC: Apply lifting (C++)") with contextlib.ExitStack() as stack: x0 = [stack.enter_context(x.localForm()) for x in x0] diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index 56881951..a8a75fd9 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -300,7 +300,7 @@ void declare_petsc_functions(nb::module_& m) std::vector>>& a, const std::vector>>>& bcs1, - const std::vector, nb::c_contig>>& x0, + const std::vector, nb::c_contig>>& x0, U scale, std::shared_ptr>& mpc) { diff --git a/python/dolfinx_mpc/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index 6fb1ec7d..08161fae 100644 --- a/python/dolfinx_mpc/multipointconstraint.py +++ b/python/dolfinx_mpc/multipointconstraint.py @@ -175,6 +175,8 @@ def create_periodic_constraint_topological(self, V: _fem.FunctionSpace, meshtag: scale: Float for scaling bc """ bcs_ = [bc._cpp_object for bc in bcs] + if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types + scale = scale.item() if (V is self.V): mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological( self.V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, False) @@ -203,6 +205,8 @@ def create_periodic_constraint_geometrical(self, V: _fem.FunctionSpace, (Periodic constraints will be ignored for these dofs) scale: Float for scaling bc """ + if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types + scale = scale.item() bcs = [] if bcs is None else [bc._cpp_object for bc in bcs] if (V is self.V): mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical( From 6fe20bf490f42e575b6fdd20b80e5d8721d989f7 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 22:35:56 +0000 Subject: [PATCH 04/30] Remove consts --- python/dolfinx_mpc/mpc.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index a8a75fd9..d304c3c8 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -145,9 +145,9 @@ void declare_functions(nb::module_& m) m.def( "create_periodic_constraint_geometrical", [](const std::shared_ptr> V, - const std::function, nb::c_contig>( + std::function, nb::c_contig>( nb::ndarray, nb::numpy>&)>& indicator, - const std::function, nb::numpy>( + std::function, nb::numpy>( nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, @@ -189,7 +189,7 @@ void declare_functions(nb::module_& m) const std::shared_ptr>& meshtags, const int dim, - const std::function, nb::numpy>( + std::function, nb::numpy>( nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, From c2debf5706fd89a79ca3225659ee3e310547a4a2 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 23:13:24 +0000 Subject: [PATCH 05/30] Typing --- python/dolfinx_mpc/multipointconstraint.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/dolfinx_mpc/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index 08161fae..fc7f5ba9 100644 --- a/python/dolfinx_mpc/multipointconstraint.py +++ b/python/dolfinx_mpc/multipointconstraint.py @@ -176,7 +176,7 @@ def create_periodic_constraint_topological(self, V: _fem.FunctionSpace, meshtag: """ bcs_ = [bc._cpp_object for bc in bcs] if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types - scale = scale.item() + scale = scale.item() # type: ignore if (V is self.V): mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological( self.V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, False) @@ -206,7 +206,7 @@ def create_periodic_constraint_geometrical(self, V: _fem.FunctionSpace, scale: Float for scaling bc """ if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types - scale = scale.item() + scale = scale.item() # type: ignore bcs = [] if bcs is None else [bc._cpp_object for bc in bcs] if (V is self.V): mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical( From 9726797d8f1d9b65d750e04697eda6cbe0f257e4 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 23:39:59 +0000 Subject: [PATCH 06/30] Deactivate nonlinear test for now --- cpp/ContactConstraint.h | 2 +- python/dolfinx_mpc/mpc.cpp | 27 ++-- python/dolfinx_mpc/multipointconstraint.py | 4 + python/tests/test_nonlinear_assembly.py | 171 ++++++++++----------- 4 files changed, 103 insertions(+), 101 deletions(-) diff --git a/cpp/ContactConstraint.h b/cpp/ContactConstraint.h index 9d3ffff5..eb8c8e83 100644 --- a/cpp/ContactConstraint.h +++ b/cpp/ContactConstraint.h @@ -884,7 +884,7 @@ template mpc_data create_contact_inelastic_condition( const dolfinx::fem::FunctionSpace& V, dolfinx::mesh::MeshTags meshtags, std::int32_t slave_marker, - std::int32_t master_marker, const double eps2 = 1e-20) + std::int32_t master_marker, const U eps2 = 1e-20) { dolfinx::common::Timer timer("~MPC: Inelastic condition"); diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index d304c3c8..8657a217 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -1,8 +1,8 @@ -// // Copyright (C) 2020 Jørgen S. Dokken -// // -// // This file is part of DOLFINX-MPC -// // -// // SPDX-License-Identifier: MIT +// Copyright (C) 2020 Jørgen S. Dokken +// +// This file is part of DOLFINX-MPC +// +// SPDX-License-Identifier: MIT #include #include @@ -144,10 +144,10 @@ void declare_functions(nb::module_& m) &dolfinx_mpc::create_contact_inelastic_condition); m.def( "create_periodic_constraint_geometrical", - [](const std::shared_ptr> V, - std::function, nb::c_contig>( + [](std::shared_ptr> V, + const std::function, nb::c_contig>( nb::ndarray, nb::numpy>&)>& indicator, - std::function, nb::numpy>( + const std::function, nb::numpy>( nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, @@ -185,11 +185,10 @@ void declare_functions(nb::module_& m) nb::arg("collapse").noconvert()); m.def( "create_periodic_constraint_topological", - [](const std::shared_ptr>& V, - const std::shared_ptr>& - meshtags, + [](std::shared_ptr>& V, + std::shared_ptr>& meshtags, const int dim, - std::function, nb::numpy>( + const std::function, nb::numpy>( nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, @@ -299,9 +298,9 @@ void declare_petsc_functions(nb::module_& m) [](nb::ndarray, nb::c_contig> b, std::vector>>& a, const std::vector>>>& bcs1, + std::shared_ptr>>>& bcs1, const std::vector, nb::c_contig>>& x0, - U scale, + T scale, std::shared_ptr>& mpc) { std::vector> _x0; diff --git a/python/dolfinx_mpc/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index fc7f5ba9..ab6b6aab 100644 --- a/python/dolfinx_mpc/multipointconstraint.py +++ b/python/dolfinx_mpc/multipointconstraint.py @@ -327,6 +327,8 @@ def create_contact_slip_condition(self, meshtags: _mesh.MeshTags, slave_marker: normal: The function used in the dot-product of the constraint eps2: The tolerance for the squared distance between cells to be considered as a collision """ + if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types + eps2 = eps2.item() # type: ignore mpc_data = dolfinx_mpc.cpp.mpc.create_contact_slip_condition( self.V._cpp_object, meshtags._cpp_object, slave_marker, master_marker, normal._cpp_object, eps2) self.add_constraint_from_mpc_data(self.V, mpc_data) @@ -345,6 +347,8 @@ def create_contact_inelastic_condition(self, meshtags: _cpp.mesh.MeshTags_int32, master_marker: The marker of the master facets eps2: The tolerance for the squared distance between cells to be considered as a collision """ + if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types + eps2 = eps2.item() # type: ignore mpc_data = dolfinx_mpc.cpp.mpc.create_contact_inelastic_condition( self.V._cpp_object, meshtags._cpp_object, slave_marker, master_marker, eps2) self.add_constraint_from_mpc_data(self.V, mpc_data) diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index b14d7f9e..72f6ff55 100644 --- a/python/tests/test_nonlinear_assembly.py +++ b/python/tests/test_nonlinear_assembly.py @@ -97,92 +97,91 @@ def b(self) -> PETSc.Vec: # type: ignore return self._b -@pytest.mark.skipif(np.issubdtype(dolfinx.default_scalar_type, np.complexfloating), - reason="This test does not work in complex mode.") -@pytest.mark.parametrize("poly_order", [1, 2, 3]) -def test_nonlinear_poisson(poly_order): - # Solve a standard Poisson problem with known solution which has - # rotational symmetry of pi/2 at (x, y) = (0.5, 0.5). Therefore we may - # impose MPCs on those DoFs which lie on the symmetry plane(s) and test - # our numerical approximation. We do not impose any constraints at the - # rotationally degenerate point (x, y) = (0.5, 0.5). - - N_vals = np.array([4, 8, 16], dtype=np.int32) - l2_error = np.zeros_like(N_vals, dtype=np.double) - for run_no, N in enumerate(N_vals): - mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) - - V = dolfinx.fem.functionspace(mesh, ("Lagrange", poly_order)) - - def boundary(x): - return np.ones_like(x[0], dtype=np.int8) - - u_bc = dolfinx.fem.Function(V) - u_bc.x.array[:] = 0.0 - - facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary) - topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets) - zero = np.array(0, dtype=dolfinx.default_scalar_type) - bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V) - bcs = [bc] - - # Define variational problem - u = dolfinx.fem.Function(V) - v = ufl.TestFunction(V) - x = ufl.SpatialCoordinate(mesh) - - u_soln = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) - f = -ufl.div((1 + u_soln**2) * ufl.grad(u_soln)) - - F = ufl.inner((1 + u**2) * ufl.grad(u), ufl.grad(v)) * ufl.dx \ - - ufl.inner(f, v) * ufl.dx - J = ufl.derivative(F, u) - - # -- Impose the pi/2 rotational symmetry of the solution as a constraint, - # -- except at the centre DoF - def periodic_boundary(x): - eps = 1000 * np.finfo(x.dtype).resolution - return np.isclose(x[0], 0.5, atol=eps) & ( - (x[1] < 0.5 - eps) | (x[1] > 0.5 + eps)) - - def periodic_relation(x): - out_x = np.zeros_like(x) - out_x[0] = x[1] - out_x[1] = x[0] - out_x[2] = x[2] - return out_x - - mpc = dolfinx_mpc.MultiPointConstraint(V) - mpc.create_periodic_constraint_geometrical( - V, periodic_boundary, periodic_relation, bcs) - mpc.finalize() - - # Sanity check that the MPC class has some constraints to impose - num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM) - num_masters_global = mesh.comm.allreduce( - len(mpc.masters.array), op=MPI.SUM) - - assert num_slaves_global > 0 - assert num_masters_global == num_slaves_global - - problem = NonlinearMPCProblem(F, u, mpc, bcs=bcs, J=J) - solver = NewtonSolverMPC(mesh.comm, problem, mpc) - solver.atol = 1e1 * np.finfo(u.x.array.dtype).resolution - solver.rtol = 1e1 * np.finfo(u.x.array.dtype).resolution - - # Ensure the solver works with nonzero initial guess - u.interpolate(lambda x: x[0]**2 * x[1]**2) - solver.solve(u) - - l2_error_local = dolfinx.fem.assemble_scalar( - dolfinx.fem.form((u - u_soln) ** 2 * ufl.dx)) - l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM) - - l2_error[run_no] = l2_error_global**0.5 - - rates = np.log(l2_error[:-1] / l2_error[1:]) / np.log(2.0) - - assert np.all(rates > poly_order + 0.9) +# @pytest.mark.skipif(np.issubdtype(dolfinx.default_scalar_type, np.complexfloating), +# reason="This test does not work in complex mode.") +# @pytest.mark.parametrize("poly_order", [1, 2, 3]) +# def test_nonlinear_poisson(poly_order): +# # Solve a standard Poisson problem with known solution which has +# # rotational symmetry of pi/2 at (x, y) = (0.5, 0.5). Therefore we may +# # impose MPCs on those DoFs which lie on the symmetry plane(s) and test +# # our numerical approximation. We do not impose any constraints at the +# # rotationally degenerate point (x, y) = (0.5, 0.5). + +# N_vals = np.array([4, 8, 16], dtype=np.int32) +# l2_error = np.zeros_like(N_vals, dtype=np.double) +# for run_no, N in enumerate(N_vals): +# mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) + +# V = dolfinx.fem.functionspace(mesh, ("Lagrange", poly_order)) + +# def boundary(x): +# return np.ones_like(x[0], dtype=np.int8) + +# u_bc = dolfinx.fem.Function(V) +# u_bc.x.array[:] = 0.0 + +# facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary) +# topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets) +# zero = np.array(0, dtype=dolfinx.default_scalar_type) +# bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V) +# bcs = [bc] + +# # Define variational problem +# u = dolfinx.fem.Function(V) +# v = ufl.TestFunction(V) +# x = ufl.SpatialCoordinate(mesh) + +# u_soln = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) +# f = -ufl.div((1 + u_soln**2) * ufl.grad(u_soln)) + +# F = ufl.inner((1 + u**2) * ufl.grad(u), ufl.grad(v)) * ufl.dx \ +# - ufl.inner(f, v) * ufl.dx +# J = ufl.derivative(F, u) + +# # -- Impose the pi/2 rotational symmetry of the solution as a constraint, +# # -- except at the centre DoF +# def periodic_boundary(x): +# eps = 1000 * np.finfo(x.dtype).resolution +# return np.isclose(x[0], 0.5, atol=eps) & ( +# (x[1] < 0.5 - eps) | (x[1] > 0.5 + eps)) + +# def periodic_relation(x): +# out_x = np.zeros_like(x) +# out_x[0] = x[1] +# out_x[1] = x[0] +# out_x[2] = x[2] +# return out_x + +# mpc = dolfinx_mpc.MultiPointConstraint(V) +# mpc.create_periodic_constraint_geometrical( +# V, periodic_boundary, periodic_relation, bcs) +# mpc.finalize() + +# # Sanity check that the MPC class has some constraints to impose +# num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM) +# num_masters_global = mesh.comm.allreduce( +# len(mpc.masters.array), op=MPI.SUM) + +# assert num_slaves_global > 0 +# assert num_masters_global == num_slaves_global + +# problem = NonlinearMPCProblem(F, u, mpc, bcs=bcs, J=J) +# solver = NewtonSolverMPC(mesh.comm, problem, mpc) +# solver.atol = 1e1 * np.finfo(u.x.array.dtype).resolution +# solver.rtol = 1e1 * np.finfo(u.x.array.dtype).resolution + +# # Ensure the solver works with nonzero initial guess +# u.interpolate(lambda x: x[0]**2 * x[1]**2) +# solver.solve(u) +# l2_error_local = dolfinx.fem.assemble_scalar( +# dolfinx.fem.form((u - u_soln) ** 2 * ufl.dx)) +# l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM) + +# l2_error[run_no] = l2_error_global**0.5 + +# rates = np.log(l2_error[:-1] / l2_error[1:]) / np.log(2.0) + +# assert np.all(rates > poly_order + 0.9) @pytest.mark.parametrize("tensor_order", [0, 1, 2]) From 348aee0006331ee88396f46825f84a9cc05f812d Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 23:41:59 +0000 Subject: [PATCH 07/30] Update CI --- .github/workflows/test_mpc.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test_mpc.yml b/.github/workflows/test_mpc.yml index 6107ebe9..baf680d5 100644 --- a/.github/workflows/test_mpc.yml +++ b/.github/workflows/test_mpc.yml @@ -99,7 +99,7 @@ jobs: cmake --install build-dir - name: Install DOLFINx-MPC (Python) - run: CXX_FLAGS="${MPC_CMAKE_CXX_FLAGS}" python3 -m pip -v install -e python/[test] + run: python3 -m pip -v install --config-settings=cmake.build-type=${MPC_BUILD_MODE} --no-build-isolation -e python/[test] - name: Run tests (serial) From 8a120d437b7bce363c23a68096e05af614bf758d Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sat, 18 Nov 2023 23:49:15 +0000 Subject: [PATCH 08/30] Turn of ufl import --- python/tests/test_nonlinear_assembly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index 72f6ff55..7ec8cdb7 100644 --- a/python/tests/test_nonlinear_assembly.py +++ b/python/tests/test_nonlinear_assembly.py @@ -12,7 +12,7 @@ import dolfinx.nls.petsc import numpy as np import pytest -import ufl +# import ufl from mpi4py import MPI from petsc4py import PETSc From 4f210deaf0b1d39b298118747cde9cfaf5db0323 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 09:53:27 +0100 Subject: [PATCH 09/30] Fix lifting default arg --- python/dolfinx_mpc/assemble_vector.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index 1ed89601..99df3955 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -16,12 +16,12 @@ from petsc4py import PETSc as _PETSc import numpy import dolfinx_mpc.cpp - -from .multipointconstraint import MultiPointConstraint +from dolfinx import default_scalar_type +from .multipointconstraint import MultiPointConstraint, _float_classes def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.DirichletBC]], # type: ignore - constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], scale: float = 1.0): # type: ignore + constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], scale: _float_classes = default_scalar_type(1.0)): # type: ignore """ Apply lifting to vector b, i.e. :math:`b = b - scale \\cdot K^T (A_j (g_j - x0_j))` From ce45a24a05744964f5344ffdd80d1c88002e7d50 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 10:06:10 +0100 Subject: [PATCH 10/30] flake8 --- python/dolfinx_mpc/assemble_vector.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index 99df3955..0d5ab12d 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -21,7 +21,8 @@ def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.DirichletBC]], # type: ignore - constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], scale: _float_classes = default_scalar_type(1.0)): # type: ignore + constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], + scale: _float_classes = default_scalar_type(1.0)): # type: ignore """ Apply lifting to vector b, i.e. :math:`b = b - scale \\cdot K^T (A_j (g_j - x0_j))` From 5e893bfcaf6a434e761aecbeb9a6c0c172a7ca57 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 09:40:09 +0000 Subject: [PATCH 11/30] Remove muddling build types --- python/CMakeLists.txt | 5 ----- 1 file changed, 5 deletions(-) diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 34efc54d..e61b726c 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -7,11 +7,6 @@ set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) -if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) - set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE) - set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo") -endif() - find_package(Python COMPONENTS Interpreter Development REQUIRED) # Detect the installed nanobind package and import it into CMake From 5132a7d822fa0ab0cdde768b70ae34d5971ab814 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 09:51:51 +0000 Subject: [PATCH 12/30] Flake8 + mypy --- python/dolfinx_mpc/assemble_vector.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index 0d5ab12d..592e5827 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -21,7 +21,7 @@ def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.DirichletBC]], # type: ignore - constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], + constraint: MultiPointConstraint, x0: List[_PETSc.Vec] = [], # type: ignore scale: _float_classes = default_scalar_type(1.0)): # type: ignore """ Apply lifting to vector b, i.e. @@ -36,7 +36,7 @@ def apply_lifting(b: _PETSc.Vec, form: List[_fem.Form], bcs: List[List[_fem.Diri scale: Scaling for lifting """ if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types - scale = scale.item() + scale = scale.item() # type: ignore t = Timer("~MPC: Apply lifting (C++)") with contextlib.ExitStack() as stack: x0 = [stack.enter_context(x.localForm()) for x in x0] From f20a9142c832a6a095908e514b80116b0953c041 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 11:36:40 +0000 Subject: [PATCH 13/30] Mypy --- python/tests/test_nonlinear_assembly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index 7ec8cdb7..994fc427 100644 --- a/python/tests/test_nonlinear_assembly.py +++ b/python/tests/test_nonlinear_assembly.py @@ -35,7 +35,7 @@ def F(self, x: PETSc.Vec, F: PETSc.Vec): # type: ignore # Apply boundary condition dolfinx_mpc.apply_lifting(F, [self._a], bcs=[self.bcs], - constraint=self.mpc, x0=[x], scale=-1.0) + constraint=self.mpc, x0=[x], scale=dolfinx.default_scalar_type(-1.0)) F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # type: ignore dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0) From 67c05e3c7c0257978cc5cd08ff5ffa3ee774cc45 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 12:29:36 +0000 Subject: [PATCH 14/30] Lifting type fixes --- cpp/lifting.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/lifting.h b/cpp/lifting.h index 032f9911..0f97952c 100644 --- a/cpp/lifting.h +++ b/cpp/lifting.h @@ -470,7 +470,8 @@ void apply_lifting( a, const std::vector>>>>& bcs1, - const std::vector>>& x0, double scale, + const std::vector>>& x0, + std::complex scale, const std::shared_ptr< const dolfinx_mpc::MultiPointConstraint, double>>& @@ -583,7 +584,8 @@ void apply_lifting( const std::vector>>>>& bcs1, - const std::vector>>& x0, float scale, + const std::vector>>& x0, + std::complex scale, const std::shared_ptr< const dolfinx_mpc::MultiPointConstraint, float>>& From ecbf72a177b5b31a696125d821bc3d2a8e6c68f9 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 18:51:01 +0000 Subject: [PATCH 15/30] Switch if order --- python/dolfinx_mpc/assemble_matrix.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/python/dolfinx_mpc/assemble_matrix.py b/python/dolfinx_mpc/assemble_matrix.py index e7ba5e6c..bcff4e9d 100644 --- a/python/dolfinx_mpc/assemble_matrix.py +++ b/python/dolfinx_mpc/assemble_matrix.py @@ -72,18 +72,16 @@ def create_sparsity_pattern(form: _fem.Form, mpc: For square forms, the MPC. For rectangular forms a list of 2 MPCs on axis 0 & 1, respectively """ - if isinstance(mpc, MultiPointConstraint): - mpc._not_finalized() - return cpp.mpc.create_sparsity_pattern(form._cpp_object, mpc._cpp_object, - mpc._cpp_object) - elif isinstance(mpc, list): + if isinstance(mpc, list): assert len(mpc) == 2 for mpc_ in mpc: mpc_._not_finalized() return cpp.mpc.create_sparsity_pattern(form._cpp_object, mpc[0]._cpp_object, mpc[1]._cpp_object) else: - raise ValueError(f"Unknown input type {type(mpc)}") + mpc._not_finalized() + return cpp.mpc.create_sparsity_pattern(form._cpp_object, mpc._cpp_object, + mpc._cpp_object) def create_matrix_nest( From 271d598070eecb4da829233042b65f4c9a3ae17e Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 20:44:41 +0000 Subject: [PATCH 16/30] mypy --- python/dolfinx_mpc/assemble_matrix.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/dolfinx_mpc/assemble_matrix.py b/python/dolfinx_mpc/assemble_matrix.py index bcff4e9d..0d176691 100644 --- a/python/dolfinx_mpc/assemble_matrix.py +++ b/python/dolfinx_mpc/assemble_matrix.py @@ -75,13 +75,13 @@ def create_sparsity_pattern(form: _fem.Form, if isinstance(mpc, list): assert len(mpc) == 2 for mpc_ in mpc: - mpc_._not_finalized() + mpc_._not_finalized() # type: ignore return cpp.mpc.create_sparsity_pattern(form._cpp_object, mpc[0]._cpp_object, mpc[1]._cpp_object) else: - mpc._not_finalized() - return cpp.mpc.create_sparsity_pattern(form._cpp_object, mpc._cpp_object, - mpc._cpp_object) + mpc._not_finalized() # type: ignore + return cpp.mpc.create_sparsity_pattern(form._cpp_object, mpc._cpp_object, # type: ignore + mpc._cpp_object) # type: ignore def create_matrix_nest( From f23c729df2402fa6fa23f077321e9bdad2ae0c18 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 20:50:37 +0000 Subject: [PATCH 17/30] Renable tests due to fixes in https://github.com/FEniCS/dolfinx/pull/2899/ --- python/tests/test_nonlinear_assembly.py | 172 ++++++++++++------------ 1 file changed, 86 insertions(+), 86 deletions(-) diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index 994fc427..6a137642 100644 --- a/python/tests/test_nonlinear_assembly.py +++ b/python/tests/test_nonlinear_assembly.py @@ -12,7 +12,7 @@ import dolfinx.nls.petsc import numpy as np import pytest -# import ufl +import ufl from mpi4py import MPI from petsc4py import PETSc @@ -97,91 +97,91 @@ def b(self) -> PETSc.Vec: # type: ignore return self._b -# @pytest.mark.skipif(np.issubdtype(dolfinx.default_scalar_type, np.complexfloating), -# reason="This test does not work in complex mode.") -# @pytest.mark.parametrize("poly_order", [1, 2, 3]) -# def test_nonlinear_poisson(poly_order): -# # Solve a standard Poisson problem with known solution which has -# # rotational symmetry of pi/2 at (x, y) = (0.5, 0.5). Therefore we may -# # impose MPCs on those DoFs which lie on the symmetry plane(s) and test -# # our numerical approximation. We do not impose any constraints at the -# # rotationally degenerate point (x, y) = (0.5, 0.5). - -# N_vals = np.array([4, 8, 16], dtype=np.int32) -# l2_error = np.zeros_like(N_vals, dtype=np.double) -# for run_no, N in enumerate(N_vals): -# mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) - -# V = dolfinx.fem.functionspace(mesh, ("Lagrange", poly_order)) - -# def boundary(x): -# return np.ones_like(x[0], dtype=np.int8) - -# u_bc = dolfinx.fem.Function(V) -# u_bc.x.array[:] = 0.0 - -# facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary) -# topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets) -# zero = np.array(0, dtype=dolfinx.default_scalar_type) -# bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V) -# bcs = [bc] - -# # Define variational problem -# u = dolfinx.fem.Function(V) -# v = ufl.TestFunction(V) -# x = ufl.SpatialCoordinate(mesh) - -# u_soln = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) -# f = -ufl.div((1 + u_soln**2) * ufl.grad(u_soln)) - -# F = ufl.inner((1 + u**2) * ufl.grad(u), ufl.grad(v)) * ufl.dx \ -# - ufl.inner(f, v) * ufl.dx -# J = ufl.derivative(F, u) - -# # -- Impose the pi/2 rotational symmetry of the solution as a constraint, -# # -- except at the centre DoF -# def periodic_boundary(x): -# eps = 1000 * np.finfo(x.dtype).resolution -# return np.isclose(x[0], 0.5, atol=eps) & ( -# (x[1] < 0.5 - eps) | (x[1] > 0.5 + eps)) - -# def periodic_relation(x): -# out_x = np.zeros_like(x) -# out_x[0] = x[1] -# out_x[1] = x[0] -# out_x[2] = x[2] -# return out_x - -# mpc = dolfinx_mpc.MultiPointConstraint(V) -# mpc.create_periodic_constraint_geometrical( -# V, periodic_boundary, periodic_relation, bcs) -# mpc.finalize() - -# # Sanity check that the MPC class has some constraints to impose -# num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM) -# num_masters_global = mesh.comm.allreduce( -# len(mpc.masters.array), op=MPI.SUM) - -# assert num_slaves_global > 0 -# assert num_masters_global == num_slaves_global - -# problem = NonlinearMPCProblem(F, u, mpc, bcs=bcs, J=J) -# solver = NewtonSolverMPC(mesh.comm, problem, mpc) -# solver.atol = 1e1 * np.finfo(u.x.array.dtype).resolution -# solver.rtol = 1e1 * np.finfo(u.x.array.dtype).resolution - -# # Ensure the solver works with nonzero initial guess -# u.interpolate(lambda x: x[0]**2 * x[1]**2) -# solver.solve(u) -# l2_error_local = dolfinx.fem.assemble_scalar( -# dolfinx.fem.form((u - u_soln) ** 2 * ufl.dx)) -# l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM) - -# l2_error[run_no] = l2_error_global**0.5 - -# rates = np.log(l2_error[:-1] / l2_error[1:]) / np.log(2.0) - -# assert np.all(rates > poly_order + 0.9) +@pytest.mark.skipif(np.issubdtype(dolfinx.default_scalar_type, np.complexfloating), + reason="This test does not work in complex mode.") +@pytest.mark.parametrize("poly_order", [1, 2, 3]) +def test_nonlinear_poisson(poly_order): + # Solve a standard Poisson problem with known solution which has + # rotational symmetry of pi/2 at (x, y) = (0.5, 0.5). Therefore we may + # impose MPCs on those DoFs which lie on the symmetry plane(s) and test + # our numerical approximation. We do not impose any constraints at the + # rotationally degenerate point (x, y) = (0.5, 0.5). + + N_vals = np.array([4, 8, 16], dtype=np.int32) + l2_error = np.zeros_like(N_vals, dtype=np.double) + for run_no, N in enumerate(N_vals): + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) + + V = dolfinx.fem.functionspace(mesh, ("Lagrange", poly_order)) + + def boundary(x): + return np.ones_like(x[0], dtype=np.int8) + + u_bc = dolfinx.fem.Function(V) + u_bc.x.array[:] = 0.0 + + facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary) + topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets) + zero = np.array(0, dtype=dolfinx.default_scalar_type) + bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V) + bcs = [bc] + + # Define variational problem + u = dolfinx.fem.Function(V) + v = ufl.TestFunction(V) + x = ufl.SpatialCoordinate(mesh) + + u_soln = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) + f = -ufl.div((1 + u_soln**2) * ufl.grad(u_soln)) + + F = ufl.inner((1 + u**2) * ufl.grad(u), ufl.grad(v)) * ufl.dx \ + - ufl.inner(f, v) * ufl.dx + J = ufl.derivative(F, u) + + # -- Impose the pi/2 rotational symmetry of the solution as a constraint, + # -- except at the centre DoF + def periodic_boundary(x): + eps = 1000 * np.finfo(x.dtype).resolution + return np.isclose(x[0], 0.5, atol=eps) & ( + (x[1] < 0.5 - eps) | (x[1] > 0.5 + eps)) + + def periodic_relation(x): + out_x = np.zeros_like(x) + out_x[0] = x[1] + out_x[1] = x[0] + out_x[2] = x[2] + return out_x + + mpc = dolfinx_mpc.MultiPointConstraint(V) + mpc.create_periodic_constraint_geometrical( + V, periodic_boundary, periodic_relation, bcs) + mpc.finalize() + + # Sanity check that the MPC class has some constraints to impose + num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM) + num_masters_global = mesh.comm.allreduce( + len(mpc.masters.array), op=MPI.SUM) + + assert num_slaves_global > 0 + assert num_masters_global == num_slaves_global + + problem = NonlinearMPCProblem(F, u, mpc, bcs=bcs, J=J) + solver = NewtonSolverMPC(mesh.comm, problem, mpc) + solver.atol = 1e1 * np.finfo(u.x.array.dtype).resolution + solver.rtol = 1e1 * np.finfo(u.x.array.dtype).resolution + + # Ensure the solver works with nonzero initial guess + u.interpolate(lambda x: x[0]**2 * x[1]**2) + solver.solve(u) + l2_error_local = dolfinx.fem.assemble_scalar( + dolfinx.fem.form((u - u_soln) ** 2 * ufl.dx)) + l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM) + + l2_error[run_no] = l2_error_global**0.5 + + rates = np.log(l2_error[:-1] / l2_error[1:]) / np.log(2.0) + + assert np.all(rates > poly_order + 0.9) @pytest.mark.parametrize("tensor_order", [0, 1, 2]) From baf3f79b29e650574bfde3083a7a08d26f999e95 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 21:20:21 +0000 Subject: [PATCH 18/30] Fix assemble_vector --- python/dolfinx_mpc/assemble_vector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index 592e5827..93f0b151 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -69,7 +69,7 @@ def assemble_vector(form: ufl.form.Form, constraint: MultiPointConstraint, t = Timer("~MPC: Assemble vector (C++)") with b.localForm() as b_local: b_local.set(0.0) - dolfinx_mpc.cpp.mpc.assemble_vector(b_local, form._cpp_object, + dolfinx_mpc.cpp.mpc.assemble_vector(b_local.array_w, form._cpp_object, constraint._cpp_object) t.stop() return b From 8ce3d7905c2c3d813ff638a577ea618e39850946 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 21:32:03 +0000 Subject: [PATCH 19/30] update action --- .github/actions/install-dolfinx/action.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/actions/install-dolfinx/action.yaml b/.github/actions/install-dolfinx/action.yaml index 694ffd24..bd2b8f70 100644 --- a/.github/actions/install-dolfinx/action.yaml +++ b/.github/actions/install-dolfinx/action.yaml @@ -70,4 +70,4 @@ runs: - name: Build Python interface (dolfinx) shell: bash -el {0} - run: PETSC_DIR=${{ inputs.petsc_dir }} PETSC_ARCH=${{ inputs.petsc_arch }} BUILD_TYPE="Release" python3 -m pip -v install ./dolfinx/python/ + run: PETSC_DIR=${{ inputs.petsc_dir }} PETSC_ARCH=${{ inputs.petsc_arch }} python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./dolfinx/python From 984b34e9b89f69c03e5f385020cb5f3a8040cb04 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Sun, 19 Nov 2023 21:37:11 +0000 Subject: [PATCH 20/30] More workflows --- .github/workflows/build_docs.yml | 2 +- .github/workflows/sonarcloud.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml index c9f14e09..a8aedf74 100644 --- a/.github/workflows/build_docs.yml +++ b/.github/workflows/build_docs.yml @@ -36,7 +36,7 @@ jobs: cmake --install build-dir - name: Install DOLFINx-MPC (Python) - run: python3 -m pip -v install python/[docs] + run: python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./python/[docs] - name: Build docs run: jupyter book build . diff --git a/.github/workflows/sonarcloud.yml b/.github/workflows/sonarcloud.yml index 83529a19..41a4bd43 100644 --- a/.github/workflows/sonarcloud.yml +++ b/.github/workflows/sonarcloud.yml @@ -81,7 +81,7 @@ jobs: cmake --install build-dir - name: Install DOLFINx-MPC (Python) - run: python3 -m pip -v install python/ + run: python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./python - name: Run sonar-scanner env: From 138d0f6cfa5207575084194425492b5d15349086 Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 08:55:12 +0000 Subject: [PATCH 21/30] FIx connectivity --- cpp/PeriodicConstraint.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/cpp/PeriodicConstraint.h b/cpp/PeriodicConstraint.h index 62c38ed7..443d8efb 100644 --- a/cpp/PeriodicConstraint.h +++ b/cpp/PeriodicConstraint.h @@ -532,6 +532,11 @@ dolfinx_mpc::mpc_data geometrical_condition( { std::vector slave_blocks = dolfinx::fem::locate_dofs_geometrical(*V, indicator); + std::stringstream cc; + for (auto b : slave_blocks) + cc << b << " "; + cc << "\n"; + std::cout << cc.str(); reduced_blocks.reserve(slave_blocks.size()); // Remove blocks in Dirichlet bcs std::vector bc_marker @@ -568,7 +573,8 @@ dolfinx_mpc::mpc_data topological_condition( T scale, bool collapse) { std::vector entities = meshtag->find(tag); - + V->mesh()->topology_mutable()->create_connectivity( + meshtag->dim(), V->mesh()->topology()->dim()); if (collapse) { // Locate dofs in sub and parent space From 005b8eb164e133a24a5b7c63b61523333a6da935 Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 12:29:41 +0000 Subject: [PATCH 22/30] Fix axtion and tolerance in periodic constraint --- .github/actions/install-dolfinx/action.yaml | 4 +++- cpp/PeriodicConstraint.h | 10 +++------- python/tests/test_nonlinear_assembly.py | 7 ++----- 3 files changed, 8 insertions(+), 13 deletions(-) diff --git a/.github/actions/install-dolfinx/action.yaml b/.github/actions/install-dolfinx/action.yaml index bd2b8f70..ec595b6d 100644 --- a/.github/actions/install-dolfinx/action.yaml +++ b/.github/actions/install-dolfinx/action.yaml @@ -70,4 +70,6 @@ runs: - name: Build Python interface (dolfinx) shell: bash -el {0} - run: PETSC_DIR=${{ inputs.petsc_dir }} PETSC_ARCH=${{ inputs.petsc_arch }} python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./dolfinx/python + run: | + python3 -m pip -r install python/build-requirements.txt + PETSC_DIR=${{ inputs.petsc_dir }} PETSC_ARCH=${{ inputs.petsc_arch }} python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./dolfinx/python diff --git a/cpp/PeriodicConstraint.h b/cpp/PeriodicConstraint.h index 443d8efb..a4f62dd0 100644 --- a/cpp/PeriodicConstraint.h +++ b/cpp/PeriodicConstraint.h @@ -129,13 +129,13 @@ dolfinx_mpc::mpc_data _create_periodic_condition( // Create bounding-box tree over owned cells std::vector r(num_cells_local); std::iota(r.begin(), r.end(), 0); - dolfinx::geometry::BoundingBoxTree tree(*mesh.get(), tdim, r, 1e-15); + dolfinx::geometry::BoundingBoxTree tree(*mesh.get(), tdim, r, tol); auto process_tree = tree.create_global_tree(mesh->comm()); auto colliding_bbox_processes = dolfinx::geometry::compute_collisions(process_tree, mapped_T_b); std::vector local_cell_collisions - = dolfinx_mpc::find_local_collisions(*mesh, tree, mapped_T_b, 1e-20); + = dolfinx_mpc::find_local_collisions(*mesh, tree, mapped_T_b, tol); dolfinx::common::Timer t0("~~Periodic: Local cell and eval basis"); auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( V, mapped_T_b, local_cell_collisions); @@ -532,11 +532,7 @@ dolfinx_mpc::mpc_data geometrical_condition( { std::vector slave_blocks = dolfinx::fem::locate_dofs_geometrical(*V, indicator); - std::stringstream cc; - for (auto b : slave_blocks) - cc << b << " "; - cc << "\n"; - std::cout << cc.str(); + reduced_blocks.reserve(slave_blocks.size()); // Remove blocks in Dirichlet bcs std::vector bc_marker diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index 6a137642..474aa010 100644 --- a/python/tests/test_nonlinear_assembly.py +++ b/python/tests/test_nonlinear_assembly.py @@ -113,14 +113,11 @@ def test_nonlinear_poisson(poly_order): mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N) V = dolfinx.fem.functionspace(mesh, ("Lagrange", poly_order)) - - def boundary(x): - return np.ones_like(x[0], dtype=np.int8) - u_bc = dolfinx.fem.Function(V) u_bc.x.array[:] = 0.0 - facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary) + mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) + facets = dolfinx.mesh.exterior_facet_indices(mesh.topology) topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets) zero = np.array(0, dtype=dolfinx.default_scalar_type) bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V) From bed25e8c76e493022e594a33106c34d58a46806a Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 12:37:11 +0000 Subject: [PATCH 23/30] Tweak action --- .github/actions/install-dolfinx/action.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/actions/install-dolfinx/action.yaml b/.github/actions/install-dolfinx/action.yaml index ec595b6d..cbbcb0fd 100644 --- a/.github/actions/install-dolfinx/action.yaml +++ b/.github/actions/install-dolfinx/action.yaml @@ -71,5 +71,5 @@ runs: - name: Build Python interface (dolfinx) shell: bash -el {0} run: | - python3 -m pip -r install python/build-requirements.txt + python3 -m pip install -r ./python/build-requirements.txt PETSC_DIR=${{ inputs.petsc_dir }} PETSC_ARCH=${{ inputs.petsc_arch }} python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./dolfinx/python From 76517a831e40419133ec232da08c4ba57a4880ad Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 12:41:02 +0000 Subject: [PATCH 24/30] More tweaks --- .github/actions/install-dolfinx/action.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/actions/install-dolfinx/action.yaml b/.github/actions/install-dolfinx/action.yaml index cbbcb0fd..aef39fd8 100644 --- a/.github/actions/install-dolfinx/action.yaml +++ b/.github/actions/install-dolfinx/action.yaml @@ -71,5 +71,5 @@ runs: - name: Build Python interface (dolfinx) shell: bash -el {0} run: | - python3 -m pip install -r ./python/build-requirements.txt + python3 -m pip install -r ./dolfinx/python/build-requirements.txt PETSC_DIR=${{ inputs.petsc_dir }} PETSC_ARCH=${{ inputs.petsc_arch }} python3 -m pip -v install --config-settings=cmake.build-type="Release" --no-build-isolation ./dolfinx/python From e0f91b7e6e729a95d57600af176ef64a1df45749 Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 13:34:41 +0000 Subject: [PATCH 25/30] Various fixes --- python/benchmarks/bench_periodic.py | 3 ++- python/demos/demo_elasticity.py | 3 ++- python/demos/demo_elasticity_disconnect.py | 1 + python/demos/demo_elasticity_disconnect_2D.py | 1 + python/demos/demo_periodic3d_topological.py | 23 +++++++++++-------- python/dolfinx_mpc/utils/mpc_utils.py | 4 ++-- 6 files changed, 21 insertions(+), 14 deletions(-) diff --git a/python/benchmarks/bench_periodic.py b/python/benchmarks/bench_periodic.py index a7b446b0..aebe3a35 100644 --- a/python/benchmarks/bench_periodic.py +++ b/python/benchmarks/bench_periodic.py @@ -20,7 +20,8 @@ from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.fem import (Function, dirichletbc, form, functionspace, - locate_dofs_geometrical, set_bc) + locate_dofs_geometrical) +from dolfinx.fem.petsc import set_bc from dolfinx.io import XDMFFile from dolfinx.mesh import (CellType, create_unit_cube, locate_entities_boundary, meshtags) diff --git a/python/demos/demo_elasticity.py b/python/demos/demo_elasticity.py index 9e38cd2e..06718fcc 100644 --- a/python/demos/demo_elasticity.py +++ b/python/demos/demo_elasticity.py @@ -135,7 +135,8 @@ def l2b(li): master_owner = mpc._cpp_object.owners.links(slave)[0] _masters = mpc.masters master = _masters.links(slave)[0] - glob_master = mpc.function_space.dofmap.index_map.local_to_global([master // bs])[0] + glob_master = mpc.function_space.dofmap.index_map.local_to_global( + np.array([master // bs], dtype=np.int32))[0] coeffs, offs = mpc.coefficients() master_data = [glob_master * bs + master % bs, coeffs[offs[slave]:offs[slave + 1]][0]] diff --git a/python/demos/demo_elasticity_disconnect.py b/python/demos/demo_elasticity_disconnect.py index 8510ff1a..e9b5bc1b 100644 --- a/python/demos/demo_elasticity_disconnect.py +++ b/python/demos/demo_elasticity_disconnect.py @@ -112,6 +112,7 @@ fdim = tdim - 1 DG0 = functionspace(mesh, ("DG", 0)) +mesh.topology.create_connectivity(tdim, tdim) outer_dofs = locate_dofs_topological(DG0, tdim, ct.find(outer_tag)) inner_dofs = locate_dofs_topological(DG0, tdim, ct.find(inner_tag)) diff --git a/python/demos/demo_elasticity_disconnect_2D.py b/python/demos/demo_elasticity_disconnect_2D.py index d205828f..a9c96410 100644 --- a/python/demos/demo_elasticity_disconnect_2D.py +++ b/python/demos/demo_elasticity_disconnect_2D.py @@ -65,6 +65,7 @@ # Locate cells with different elasticity parameters DG0 = functionspace(mesh, ("DG", 0)) +mesh.topology.create_connectivity(tdim, tdim) left_dofs = locate_dofs_topological(DG0, tdim, ct.find(left_tag)) right_dofs = locate_dofs_topological(DG0, tdim, ct.find(right_tag)) diff --git a/python/demos/demo_periodic3d_topological.py b/python/demos/demo_periodic3d_topological.py index 1c9f21f2..5aa11454 100644 --- a/python/demos/demo_periodic3d_topological.py +++ b/python/demos/demo_periodic3d_topological.py @@ -51,10 +51,11 @@ def demo_periodic3D(celltype: CellType): N = 10 mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, CellType.hexahedron) V = fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, ))) + tol = float(5e2 * np.finfo(default_scalar_type).resolution) def dirichletboundary(x: NDArray[Union[np.float32, np.float64]]) -> NDArray[np.bool_]: - return np.logical_or(np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)), - np.logical_or(np.isclose(x[2], 0), np.isclose(x[2], 1))) + return np.logical_or(np.logical_or(np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol)), + np.logical_or(np.isclose(x[2], 0, atol=tol), np.isclose(x[2], 1, atol=tol))) # Create Dirichlet boundary condition zero = default_scalar_type([0, 0, 0]) @@ -63,7 +64,10 @@ def dirichletboundary(x: NDArray[Union[np.float32, np.float64]]) -> NDArray[np.b bcs = [bc] def PeriodicBoundary(x): - return np.isclose(x[0], 1) + """ + Full surface minus dofs constrained by BCs + """ + return np.isclose(x[0], 1, atol=tol) facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, PeriodicBoundary) arg_sort = np.argsort(facets) @@ -94,7 +98,6 @@ def periodic_relation(x): rhs = inner(f, v) * dx petsc_options: Dict[str, Union[str, float, int]] - tol = float(5e2 * np.finfo(default_scalar_type).resolution) if complex_mode or default_scalar_type == np.float32: petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} else: @@ -138,21 +141,21 @@ def periodic_relation(x): with Timer("~Demo: Verification"): dolfinx_mpc.utils.compare_mpc_lhs(org_problem.A, problem.A, mpc, root=root) dolfinx_mpc.utils.compare_mpc_rhs(org_problem.b, problem.b, mpc, root=root) - + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 # Gather LHS, RHS and solution on one process A_csr = dolfinx_mpc.utils.gather_PETScMatrix(org_problem.A, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) L_np = dolfinx_mpc.utils.gather_PETScVector(org_problem.b, root=root) u_mpc = dolfinx_mpc.utils.gather_PETScVector(u_h.vector, root=root) - if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc, rtol=tol, atol=tol) + uh_numpy = K.astype(scipy_dtype) @ d.astype(scipy_dtype) + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, rtol=tol, atol=tol) if __name__ == "__main__": diff --git a/python/dolfinx_mpc/utils/mpc_utils.py b/python/dolfinx_mpc/utils/mpc_utils.py index bce7ef14..62c3739e 100644 --- a/python/dolfinx_mpc/utils/mpc_utils.py +++ b/python/dolfinx_mpc/utils/mpc_utils.py @@ -299,7 +299,7 @@ def create_point_to_point_constraint(V, slave_point, master_point, vector=None): slave_index = np.argmax(np.abs(vector)) if is_slave_proc: assert len(slave_block) == 1 - slave_block_g = imap.local_to_global(slave_block)[0] + slave_block_g = imap.local_to_global(np.asarray(slave_block, dtype=np.int32))[0] if vector is None: slaves = np.arange(slave_block[0] * block_size, slave_block[0] * block_size + block_size, dtype=np.int32) @@ -319,7 +319,7 @@ def create_point_to_point_constraint(V, slave_point, master_point, vector=None): global_masters = None if is_master_proc: assert len(master_block) == 1 - master_block_g = imap.local_to_global(master_block)[0] + master_block_g = imap.local_to_global(np.asarray(master_block, dtype=np.int32))[0] masters_as_glob = np.arange(master_block_g * block_size, master_block_g * block_size + block_size, dtype=np.int64) else: From 61b690ff8aefc07309ba03bae9c8f9196ccd3e36 Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 13:49:55 +0000 Subject: [PATCH 26/30] More demo fixes --- python/demos/demo_periodic_geometrical.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/python/demos/demo_periodic_geometrical.py b/python/demos/demo_periodic_geometrical.py index d283c109..4560f635 100644 --- a/python/demos/demo_periodic_geometrical.py +++ b/python/demos/demo_periodic_geometrical.py @@ -38,10 +38,11 @@ NY = 100 mesh = create_unit_square(MPI.COMM_WORLD, NX, NY) V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) +tol = 250 * np.finfo(default_scalar_type).resolution def dirichletboundary(x): - return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)) + return np.logical_or(np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol)) # Create Dirichlet boundary condition @@ -53,11 +54,11 @@ def dirichletboundary(x): def periodic_boundary(x): - return np.isclose(x[0], 1) + return np.isclose(x[0], 1, atol=tol) def periodic_relation(x): - out_x = np.zeros(x.shape) + out_x = np.zeros_like(x) out_x[0] = 1 - x[0] out_x[1] = x[1] out_x[2] = x[2] @@ -152,7 +153,8 @@ def periodic_relation(x): with Timer("~Demo: Verification"): dolfinx_mpc.utils.compare_mpc_lhs(A_org, problem._A, mpc, root=root) dolfinx_mpc.utils.compare_mpc_rhs(L_org, problem._b, mpc, root=root) - + is_complex = np.issubdtype(default_scalar_type, np.complexfloating) # type: ignore + scipy_dtype = np.complex128 if is_complex else np.float64 # Gather LHS, RHS and solution on one process A_csr = dolfinx_mpc.utils.gather_PETScMatrix(A_org, root=root) K = dolfinx_mpc.utils.gather_transformation_matrix(mpc, root=root) @@ -160,12 +162,12 @@ def periodic_relation(x): u_mpc = dolfinx_mpc.utils.gather_PETScVector(uh.vector, root=root) if MPI.COMM_WORLD.rank == root: - KTAK = K.T * A_csr * K - reduced_L = K.T @ L_np + KTAK = K.T.astype(scipy_dtype) * A_csr.astype(scipy_dtype) * K.astype(scipy_dtype) + reduced_L = K.T.astype(scipy_dtype) @ L_np.astype(scipy_dtype) # Solve linear system d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector - uh_numpy = K @ d - assert np.allclose(uh_numpy, u_mpc, atol=float(np.finfo(u_mpc.dtype).resolution)) + uh_numpy = K.astype(scipy_dtype) @ d.astype(scipy_dtype) + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, atol=tol) list_timings(MPI.COMM_WORLD, [TimingType.wall]) L_org.destroy() From 3c4867076ada7de35fdd232e682420a6c75dcf9f Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 13:54:30 +0000 Subject: [PATCH 27/30] Mypy --- python/demos/demo_periodic_geometrical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/demos/demo_periodic_geometrical.py b/python/demos/demo_periodic_geometrical.py index 4560f635..93fba970 100644 --- a/python/demos/demo_periodic_geometrical.py +++ b/python/demos/demo_periodic_geometrical.py @@ -168,6 +168,6 @@ def periodic_relation(x): d = scipy.sparse.linalg.spsolve(KTAK, reduced_L) # Back substitution to full solution vector uh_numpy = K.astype(scipy_dtype) @ d.astype(scipy_dtype) - assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, atol=tol) + assert np.allclose(uh_numpy.astype(u_mpc.dtype), u_mpc, atol=float(tol)) list_timings(MPI.COMM_WORLD, [TimingType.wall]) L_org.destroy() From c405ff0fe66dda505fe85f8234592a488025a16f Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 14:10:56 +0000 Subject: [PATCH 28/30] Fix edge example --- python/benchmarks/bench_elasticity_edge.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/benchmarks/bench_elasticity_edge.py b/python/benchmarks/bench_elasticity_edge.py index dbe663c4..7d2706fc 100644 --- a/python/benchmarks/bench_elasticity_edge.py +++ b/python/benchmarks/bench_elasticity_edge.py @@ -15,7 +15,8 @@ from dolfinx import default_scalar_type from dolfinx.common import Timer, TimingType, list_timings from dolfinx.fem import (Constant, Function, functionspace, dirichletbc, form, - locate_dofs_topological, set_bc) + locate_dofs_topological) +from dolfinx.fem.petsc import set_bc from dolfinx.io import XDMFFile from dolfinx.mesh import (CellType, create_unit_cube, locate_entities_boundary, meshtags) From 4b9c660daf5c6cfcec1042b2f0f414753d4a7cb6 Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 14:56:56 +0000 Subject: [PATCH 29/30] Remove reshape --- python/benchmarks/bench_contact_3D.py | 8 ++++---- python/demos/create_and_export_mesh.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/python/benchmarks/bench_contact_3D.py b/python/benchmarks/bench_contact_3D.py index 0e71b759..1bf6454c 100644 --- a/python/benchmarks/bench_contact_3D.py +++ b/python/benchmarks/bench_contact_3D.py @@ -80,12 +80,12 @@ def over_plane(p0, p1, p2): tdim0 = mesh0.topology.dim num_cells0 = mesh0.topology.index_map(tdim0).size_local - cells0 = entities_to_geometry(mesh0._cpp_object, tdim0, np.arange( - num_cells0, dtype=np.int32).reshape((-1, 1)), False) + cells0 = entities_to_geometry(mesh0._cpp_object, tdim0, + np.arange(num_cells0, dtype=np.int32), False) tdim1 = mesh1.topology.dim num_cells1 = mesh1.topology.index_map(tdim1).size_local - cells1 = entities_to_geometry(mesh1._cpp_object, tdim1, np.arange( - num_cells1, dtype=np.int32).reshape((-1, 1)), False) + cells1 = entities_to_geometry(mesh1._cpp_object, tdim1, + np.arange(num_cells1, dtype=np.int32), False) cells1 += mesh0.geometry.x.shape[0] # Concatenate points and cells diff --git a/python/demos/create_and_export_mesh.py b/python/demos/create_and_export_mesh.py index ce97c817..a6bfd083 100644 --- a/python/demos/create_and_export_mesh.py +++ b/python/demos/create_and_export_mesh.py @@ -384,11 +384,11 @@ def over_line(p0, p1): tdim0 = mesh0.topology.dim num_cells0 = mesh0.topology.index_map(tdim0).size_local cells0 = _cpp.mesh.entities_to_geometry( - mesh0._cpp_object, tdim0, np.arange(num_cells0, dtype=np.int32).reshape((-1, 1)), False) + mesh0._cpp_object, tdim0, np.arange(num_cells0, dtype=np.int32), False) tdim1 = mesh1.topology.dim num_cells1 = mesh1.topology.index_map(tdim1).size_local cells1 = _cpp.mesh.entities_to_geometry( - mesh1._cpp_object, tdim1, np.arange(num_cells1, dtype=np.int32).reshape((-1, 1)), False) + mesh1._cpp_object, tdim1, np.arange(num_cells1, dtype=np.int32), False) cells1 += mesh0.geometry.x.shape[0] cells = np.vstack([cells0, cells1]) @@ -504,11 +504,11 @@ def over_plane(p0, p1, p2): tdim0 = mesh0.topology.dim num_cells0 = mesh0.topology.index_map(tdim0).size_local cells0 = _cpp.mesh.entities_to_geometry( - mesh0._cpp_object, tdim0, np.arange(num_cells0, dtype=np.int32).reshape((-1, 1)), False) + mesh0._cpp_object, tdim0, np.arange(num_cells0, dtype=np.int32), False) tdim1 = mesh1.topology.dim num_cells1 = mesh1.topology.index_map(tdim1).size_local cells1 = _cpp.mesh.entities_to_geometry( - mesh1._cpp_object, tdim1, np.arange(num_cells1, dtype=np.int32).reshape((-1, 1)), False) + mesh1._cpp_object, tdim1, np.arange(num_cells1, dtype=np.int32), False) cells1 += mesh0.geometry.x.shape[0] cells = np.vstack([cells0, cells1]) From facbe47397f709bf7f5820530f0740543744c7a2 Mon Sep 17 00:00:00 2001 From: "Jorgen S. Dokken" Date: Mon, 20 Nov 2023 15:24:39 +0000 Subject: [PATCH 30/30] Arange for cells in midpoints --- python/benchmarks/bench_contact_3D.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/benchmarks/bench_contact_3D.py b/python/benchmarks/bench_contact_3D.py index 1bf6454c..05db7b04 100644 --- a/python/benchmarks/bench_contact_3D.py +++ b/python/benchmarks/bench_contact_3D.py @@ -128,7 +128,8 @@ def over_plane(p0, p1, p2): num_cells = mesh.topology.index_map(tdim).size_local # Find top and bottom interface facets - cell_midpoints = compute_midpoints(mesh, tdim, range(num_cells)) + cells = np.arange(num_cells, dtype=np.int32) + cell_midpoints = compute_midpoints(mesh, tdim, cells) top_cube = over_plane(if_points[:, 0], if_points[:, 1], if_points[:, 2]) for facet in i_facets: i_cells = facet_to_cell.links(facet) @@ -140,8 +141,6 @@ def over_plane(p0, p1, p2): bottom_interface.append(facet) # Create cell tags - num_cells = mesh.topology.index_map(tdim).size_local - cell_midpoints = compute_midpoints(mesh, tdim, range(num_cells)) top_cube_marker = 2 indices = [] values = []