diff --git a/.github/actions/install-dolfinx/action.yaml b/.github/actions/install-dolfinx/action.yaml index 694ffd24..aef39fd8 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 }} BUILD_TYPE="Release" python3 -m pip -v install ./dolfinx/python/ + run: | + 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 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: 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) 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..eb8c8e83 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) + std::int32_t master_marker, const U 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/PeriodicConstraint.h b/cpp/PeriodicConstraint.h index 62c38ed7..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,6 +532,7 @@ dolfinx_mpc::mpc_data geometrical_condition( { std::vector slave_blocks = dolfinx::fem::locate_dofs_geometrical(*V, indicator); + reduced_blocks.reserve(slave_blocks.size()); // Remove blocks in Dirichlet bcs std::vector bc_marker @@ -568,7 +569,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 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..0f97952c 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) @@ -469,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>>& @@ -582,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>>& 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..e61b726c 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -2,25 +2,77 @@ 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) +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 +85,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/benchmarks/bench_contact_3D.py b/python/benchmarks/bench_contact_3D.py index 0e71b759..05db7b04 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 @@ -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 = [] 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) 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/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]) 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/demos/demo_periodic_geometrical.py b/python/demos/demo_periodic_geometrical.py index d283c109..93fba970 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=float(tol)) list_timings(MPI.COMM_WORLD, [TimingType.wall]) L_org.destroy() diff --git a/python/dolfinx_mpc/assemble_matrix.py b/python/dolfinx_mpc/assemble_matrix.py index e7ba5e6c..0d176691 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() + mpc_._not_finalized() # type: ignore 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() # 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( diff --git a/python/dolfinx_mpc/assemble_vector.py b/python/dolfinx_mpc/assemble_vector.py index d30ac265..93f0b151 100644 --- a/python/dolfinx_mpc/assemble_vector.py +++ b/python/dolfinx_mpc/assemble_vector.py @@ -14,14 +14,15 @@ import ufl from dolfinx.common import Timer 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] = [], # type: ignore + 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))` @@ -34,6 +35,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() # type: ignore t = Timer("~MPC: Apply lifting (C++)") with contextlib.ExitStack() as stack: x0 = [stack.enter_context(x.localForm()) for x in x0] @@ -66,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 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..8657a217 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -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); @@ -132,13 +142,13 @@ void declare_functions(py::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 py::array_t&)>& - indicator, - const std::function(const py::array_t&)>& relation, + [](std::shared_ptr> V, + const std::function, nb::c_contig>( + 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 +162,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,95 +172,97 @@ 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, nb::arg("scale").noconvert(), + 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, - 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, + nb::arg("scale").noconvert(), nb::arg("collapse").noconvert()); } 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)); + return nb::ndarray>( + masters.data(), {masters.size()}); }) - .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)); - }) - .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 +286,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, + std::shared_ptr>>>& bcs1, + const std::vector, nb::c_contig>>& x0, + T 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 +322,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 +337,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 +345,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/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index 6fb1ec7d..ab6b6aab 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() # 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) @@ -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() # 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( @@ -323,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) @@ -341,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/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: 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) 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: diff --git a/python/tests/test_nonlinear_assembly.py b/python/tests/test_nonlinear_assembly.py index b14d7f9e..474aa010 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) @@ -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) @@ -173,7 +170,6 @@ def periodic_relation(x): # 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)