From 17c532e9f97b26ed1e66875e719f4c3291dc25bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 19 Sep 2023 17:00:39 +0200 Subject: [PATCH] Dokken/mpc generalize floating (#75) * mdspan updates + templating over mesh dtype to a larger extent * Function-space typing updates * Correct templating for things that depend on PETSc * Update changelog and fix dtypes in python * Fix bug (thanks sonarcloud) and flake8 * Typing + flake8 * Add back code --- Changelog.md | 8 +- cpp/ContactConstraint.h | 216 +++++----- cpp/PeriodicConstraint.h | 162 +++++--- cpp/SlipConstraint.h | 33 +- cpp/assemble_matrix.cpp | 101 +++-- cpp/assemble_matrix.h | 54 +++ cpp/assemble_vector.cpp | 56 ++- cpp/assemble_vector.h | 28 +- cpp/lifting.h | 130 +++++- cpp/mpc_helpers.h | 7 +- cpp/utils.cpp | 36 -- cpp/utils.h | 182 +++++---- python/benchmarks/bench_elasticity_edge.py | 3 +- python/demos/demo_elasticity_disconnect_2D.py | 17 +- python/demos/demo_periodic3d_topological.py | 2 +- python/dolfinx_mpc/mpc.cpp | 371 ++++++++++-------- python/dolfinx_mpc/multipointconstraint.py | 106 +++-- python/dolfinx_mpc/utils/mpc_utils.py | 4 +- python/tests/test_cube_contact.py | 2 +- python/tests/test_rectangular_assembly.py | 6 +- python/tests/test_surface_integral.py | 4 +- python/tests/test_vector_poisson.py | 2 +- 22 files changed, 970 insertions(+), 560 deletions(-) diff --git a/Changelog.md b/Changelog.md index 603093c9..8c70a7fa 100644 --- a/Changelog.md +++ b/Changelog.md @@ -1,7 +1,13 @@ # Changelog ## main -- Change input of `dolfinx_mpc.MultiPointConstraint.homogenize` and `dolfinx_mpc.backsubstitution` to `dolfinx.fem.Function` instead of `PETSc.Vec`. +- **API**: + - Change input of `dolfinx_mpc.MultiPointConstraint.homogenize` and `dolfinx_mpc.backsubstitution` to `dolfinx.fem.Function` instead of `PETSc.Vec`. + - Add support for more floating types (`dtype` can now be specified as input for coefficients). In theory this means one could support single precision. Not thoroughly tested. + - 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}` +- **DOLFINX API-changes**: + - Use `dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))` instead of `dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 1))` as the latter is being deprecated. + - Use `basix.ufl.element` in favor of `ufl.FiniteElement` as the latter is deprecated in DOLFINx. ## v0.6.1 (30.01.2023) - Fixes for CI diff --git a/cpp/ContactConstraint.h b/cpp/ContactConstraint.h index 13f425e6..84574f59 100644 --- a/cpp/ContactConstraint.h +++ b/cpp/ContactConstraint.h @@ -59,13 +59,13 @@ template dolfinx_mpc::mpc_data compute_master_contributions( std::span local_rems, std::span local_colliding_cell, - std::experimental::mdspan< - double, std::experimental::extents< - std::size_t, std::experimental::dynamic_extent, 3>> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> normals, std::shared_ptr> V, - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> tabulated_basis_values) { const double tol = 1e-6; @@ -204,11 +204,11 @@ locate_slave_dofs(std::shared_ptr> V, /// @param[in] block_size The block size of the index map /// @param[in] rank The rank of current process /// @returns A mpc_data struct with slaves, masters, coeffs and owners -template +template dolfinx_mpc::mpc_data compute_block_contributions( const std::vector& local_slaves, const std::vector& local_slave_blocks, - std::span normals, + std::span normals, const std::shared_ptr imap, std::int32_t block_size, int rank) { @@ -345,13 +345,12 @@ namespace dolfinx_mpc /// @param[in] nh Function containing the normal at the slave marker interface /// @param[in] eps2 The tolerance for the squared distance to be considered a /// collision -template -mpc_data create_contact_slip_condition( +template +mpc_data create_contact_slip_condition( std::shared_ptr> 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) + std::shared_ptr> nh, const U eps2 = 1e-20) { dolfinx::common::Timer timer("~MPC: Create slip constraint"); @@ -397,7 +396,7 @@ mpc_data create_contact_slip_condition( // Data structures to hold information about slave data local to process std::vector local_slaves(local_slave_blocks.size()); std::vector local_rems(local_slave_blocks.size()); - dolfinx_mpc::mpc_data mpc_local; + dolfinx_mpc::mpc_data mpc_local; // Find all local contributions to MPC, meaning: // 1. Degrees of freedom from the same block as the slave @@ -406,17 +405,19 @@ mpc_data create_contact_slip_condition( // Helper function // Determine component of each block has the largest normal value, and use // it as slave dofs to avoid zero division in constraint + // 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) + std::span normal) { std::iota(dofs.begin(), dofs.end(), block * block_size); for (int j = 0; j < gdim; ++j) normal[j] = std::real(normal_array[dofs[j]]); - double norm = std::sqrt(normal[0] * normal[0] + normal[1] * normal[1] - + normal[2] * normal[2]); + U norm = std::sqrt(normal[0] * normal[0] + normal[1] * normal[1] + + normal[2] * normal[2]); std::for_each(normal.begin(), normal.end(), [norm](auto& n) { return std::abs(n / norm); }); return std::distance(normal.begin(), @@ -424,56 +425,55 @@ mpc_data create_contact_slip_condition( }; // Determine which dof in local slave block is the actual slave - std::vector normals(3 * local_slave_blocks.size(), 0); + std::vector normals(3 * local_slave_blocks.size(), 0); assert(block_size == gdim); for (std::size_t i = 0; i < local_slave_blocks.size(); ++i) { const std::int32_t slave = local_slave_blocks[i]; const auto block = largest_normal_component( - slave, std::span(std::next(normals.begin(), 3 * i), 3)); + slave, std::span(std::next(normals.begin(), 3 * i), 3)); local_slaves[i] = block_size * slave + block; local_rems[i] = block; } // Compute local contributions to constraint using helper function // i.e. compute dot(u, n) on slave side - mpc_data mpc_in_cell - = impl::compute_block_contributions( - local_slaves, local_slave_blocks, normals, imap, block_size, rank); + mpc_data mpc_in_cell = impl::compute_block_contributions( + local_slaves, local_slave_blocks, normals, imap, block_size, rank); dolfinx::geometry::BoundingBoxTree bb_tree = impl::create_boundingbox_tree( *mesh, meshtags, master_marker, std::sqrt(eps2)); // Compute contributions on other side local to process - mpc_data mpc_master_local; + mpc_data mpc_master_local; // 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); - std::vector slave_coordinates; + std::vector slave_coordinates; std::array coord_shape; { std::tie(slave_coordinates, coord_shape) - = dolfinx_mpc::tabulate_dof_coordinates(*V, local_slave_blocks, - slave_cells); + = 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( + auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( *V, slave_coordinates, local_cell_collisions); assert(basis_shape.back() == 1); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> basis_span(basis.data(), basis_shape[0], basis_shape[1]); - std::experimental::mdspan< - double, std::experimental::extents< - std::size_t, std::experimental::dynamic_extent, 3>> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> normal_span(normals.data(), local_slave_blocks.size(), 3); - mpc_master_local = impl::compute_master_contributions( + mpc_master_local = impl::compute_master_contributions( local_rems, local_cell_collisions, normal_span, V, basis_span); } - // Find slave indices were contributions are not found on the process + // // Find slave indices were contributions are not found on the process std::vector& l_offsets = mpc_master_local.offsets; std::vector slave_indices_remote; slave_indices_remote.reserve(local_rems.size()); @@ -517,8 +517,8 @@ mpc_data create_contact_slip_condition( // and prepare data (coordinates and normals) to send to other procs const std::array send_shape = {slave_indices_remote.size(), 3}; - std::vector coordinates_send(send_shape.front() * send_shape.back()); - std::vector normals_send(send_shape.front() * send_shape.back()); + std::vector coordinates_send(send_shape.front() * send_shape.back()); + std::vector normals_send(send_shape.front() * send_shape.back()); std::vector send_rems(slave_indices_remote.size()); for (std::size_t i = 0; i < slave_indices_remote.size(); ++i) { @@ -564,35 +564,33 @@ mpc_data create_contact_slip_condition( disp3.begin() + 1); // Send slave normal and coordinate to neighbors - std::vector recv_coords(disp.back() * 3); + std::vector recv_coords(disp.back() * 3); MPI_Neighbor_allgatherv(coordinates_send.data(), (int)coordinates_send.size(), - dolfinx::MPI::mpi_type(), recv_coords.data(), + dolfinx::MPI::mpi_type(), recv_coords.data(), num_slaves_recv3.data(), disp3.data(), - dolfinx::MPI::mpi_type(), - neighborhood_comms[0]); - std::vector slave_normals(disp.back() * 3); + dolfinx::MPI::mpi_type(), neighborhood_comms[0]); + std::vector slave_normals(disp.back() * 3); MPI_Neighbor_allgatherv(normals_send.data(), (int)normals_send.size(), - dolfinx::MPI::mpi_type(), - slave_normals.data(), num_slaves_recv3.data(), - disp3.data(), dolfinx::MPI::mpi_type(), - neighborhood_comms[0]); + dolfinx::MPI::mpi_type(), slave_normals.data(), + num_slaves_recv3.data(), disp3.data(), + dolfinx::MPI::mpi_type(), neighborhood_comms[0]); // Compute off-process contributions - mpc_data remote_data; + mpc_data remote_data; { std::vector remote_cell_collisions = dolfinx_mpc::find_local_collisions(*mesh, bb_tree, recv_coords, eps2); - auto [recv_basis_values, shape] = dolfinx_mpc::evaluate_basis_functions( + auto [recv_basis_values, shape] = dolfinx_mpc::evaluate_basis_functions( *V, recv_coords, remote_cell_collisions); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> basis_span(recv_basis_values.data(), shape[0], shape[1]); - std::experimental::mdspan< - double, std::experimental::extents< - std::size_t, std::experimental::dynamic_extent, 3>> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> normal_span(slave_normals.data(), disp.back(), 3); - remote_data = impl::compute_master_contributions( + remote_data = impl::compute_master_contributions( recv_rems, remote_cell_collisions, normal_span, V, basis_span); } @@ -673,12 +671,12 @@ mpc_data create_contact_slip_condition( remote_colliding_masters.data(), inc_num_collision_masters.data(), disp_inc_masters.data(), dolfinx::MPI::mpi_type(), neighborhood_comms[1], &requests[1]); - std::vector remote_colliding_coeffs(disp_inc_masters.back()); + std::vector remote_colliding_coeffs(disp_inc_masters.back()); MPI_Ineighbor_alltoallv( remote_data.coeffs.data(), num_collision_masters.data(), - send_disp_masters.data(), dolfinx::MPI::mpi_type(), + send_disp_masters.data(), dolfinx::MPI::mpi_type(), remote_colliding_coeffs.data(), inc_num_collision_masters.data(), - disp_inc_masters.data(), dolfinx::MPI::mpi_type(), + disp_inc_masters.data(), dolfinx::MPI::mpi_type(), neighborhood_comms[1], &requests[2]); std::vector remote_colliding_owners(disp_inc_masters.back()); MPI_Ineighbor_alltoallv( @@ -735,7 +733,7 @@ mpc_data create_contact_slip_condition( std::partial_sum(num_inc_masters.begin(), num_inc_masters.end(), offproc_offsets.begin() + 1); std::vector offproc_masters(offproc_offsets.back()); - std::vector offproc_coeffs(offproc_offsets.back()); + std::vector offproc_coeffs(offproc_offsets.back()); std::vector offproc_owners(offproc_offsets.back()); std::fill(slave_found.begin(), slave_found.end(), false); @@ -780,7 +778,7 @@ mpc_data create_contact_slip_condition( // First count number of local masters std::vector& masters_offsets = mpc_local.offsets; std::vector& masters_out = mpc_local.masters; - std::vector& coefficients_out = mpc_local.coeffs; + std::vector& coefficients_out = mpc_local.coeffs; std::vector& owners_out = mpc_local.owners; std::vector num_masters_per_slave(local_slaves.size(), 0); @@ -799,7 +797,7 @@ mpc_data create_contact_slip_condition( // Reuse num_masters_per_slave for input indices std::vector local_masters(local_offsets.back()); std::vector local_owners(local_offsets.back()); - std::vector local_coeffs(local_offsets.back()); + std::vector local_coeffs(local_offsets.back()); // Insert local contributions { @@ -842,10 +840,9 @@ mpc_data create_contact_slip_condition( } } // Distribute ghost data - dolfinx_mpc::mpc_data ghost_data - = dolfinx_mpc::distribute_ghost_data( - local_slaves, local_masters, local_coeffs, local_owners, - num_masters_per_slave, imap, block_size); + dolfinx_mpc::mpc_data ghost_data = dolfinx_mpc::distribute_ghost_data( + local_slaves, local_masters, local_coeffs, local_owners, + num_masters_per_slave, imap, block_size); // Add ghost data to existing arrays const std::vector& ghost_slaves = ghost_data.slaves; @@ -857,7 +854,7 @@ mpc_data create_contact_slip_condition( const std::vector& ghost_num = ghost_data.offsets; num_masters_per_slave.insert(std::end(num_masters_per_slave), std::cbegin(ghost_num), std::cend(ghost_num)); - const std::vector& ghost_coeffs = ghost_data.coeffs; + const std::vector& ghost_coeffs = ghost_data.coeffs; local_coeffs.insert(std::end(local_coeffs), std::cbegin(ghost_coeffs), std::cend(ghost_coeffs)); const std::vector& ghost_owner_ranks = ghost_data.owners; @@ -869,7 +866,7 @@ mpc_data create_contact_slip_condition( std::partial_sum(num_masters_per_slave.begin(), num_masters_per_slave.end(), offsets.begin() + 1); - dolfinx_mpc::mpc_data output; + dolfinx_mpc::mpc_data output; output.offsets = offsets; output.masters = local_masters; output.coeffs = local_coeffs; @@ -885,8 +882,8 @@ mpc_data create_contact_slip_condition( /// @param[in] master_marker Tag for the other interface /// @param[in] eps2 The tolerance for the squared distance to be considered a /// collision -template -mpc_data create_contact_inelastic_condition( +template +mpc_data create_contact_inelastic_condition( std::shared_ptr> V, dolfinx::mesh::MeshTags meshtags, std::int32_t slave_marker, std::int32_t master_marker, const double eps2 = 1e-20) @@ -973,33 +970,34 @@ mpc_data create_contact_inelastic_condition( // 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); - std::vector slave_coordinates; + std::vector slave_coordinates; { std::array c_shape; std::tie(slave_coordinates, c_shape) - = dolfinx_mpc::tabulate_dof_coordinates(*V, local_blocks, slave_cells); + = dolfinx_mpc::tabulate_dof_coordinates(*V, local_blocks, + slave_cells); } // Loop through all masters on current processor and check if they // collide with a local master facet std::map> local_owners; std::map> local_masters; - std::map> local_coeffs; + std::map> local_coeffs; std::vector blocks_wo_local_collision; std::vector collision_to_local; { std::vector colliding_cells = dolfinx_mpc::find_local_collisions(*V->mesh(), bb_tree, slave_coordinates, eps2); - auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( + auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( *V, slave_coordinates, colliding_cells); assert(basis_shape.back() == 1); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> basis(basis_values.data(), basis_shape[0], basis_shape[1]); std::vector master_block_global; std::vector l_master; - std::vector coeff; + std::vector coeff; for (std::size_t i = 0; i < local_blocks.size(); ++i) { @@ -1013,7 +1011,7 @@ mpc_data create_contact_inelastic_condition( // Store block and non-zero basis value for (std::size_t j = 0; j < cell_blocks.size(); ++j) { - if (const PetscScalar c = basis(i, j); std::abs(c) > 1e-6) + if (const T c = basis(i, j); std::abs(c) > 1e-6) { coeff.push_back(c); l_master.push_back(cell_blocks[j]); @@ -1060,14 +1058,13 @@ mpc_data create_contact_inelastic_condition( } } // Extract coordinates and normals to distribute - std::vector distribute_coordinates(blocks_wo_local_collision.size() - * 3); + std::vector distribute_coordinates(blocks_wo_local_collision.size() * 3); for (std::size_t i = 0; i < collision_to_local.size(); ++i) { std::copy_n(std::next(slave_coordinates.begin(), 3 * collision_to_local[i]), 3, std::next(distribute_coordinates.begin(), 3 * i)); } - dolfinx_mpc::mpc_data mpc; + dolfinx_mpc::mpc_data mpc; // If serial, we only have to gather slaves, masters, coeffs in 1D // arrays @@ -1084,7 +1081,7 @@ mpc_data create_contact_inelastic_condition( } std::vector masters_out; - std::vector coeffs_out; + std::vector coeffs_out; std::vector offsets_out = {0}; std::vector slaves_out; // Flatten the maps to 1D arrays (assuming all slaves are local @@ -1155,20 +1152,20 @@ mpc_data create_contact_inelastic_condition( coordinate_disp.begin() + 1); // Send slave coordinates to neighbors - std::vector recv_coords(disp.back() * 3); - MPI_Neighbor_allgatherv( - distribute_coordinates.data(), (int)distribute_coordinates.size(), - dolfinx::MPI::mpi_type(), recv_coords.data(), - num_block_coordinates.data(), coordinate_disp.data(), - dolfinx::MPI::mpi_type(), neighborhood_comms[0]); + std::vector recv_coords(disp.back() * 3); + MPI_Neighbor_allgatherv(distribute_coordinates.data(), + (int)distribute_coordinates.size(), + dolfinx::MPI::mpi_type(), recv_coords.data(), + num_block_coordinates.data(), coordinate_disp.data(), + dolfinx::MPI::mpi_type(), neighborhood_comms[0]); // Vector for processes with slaves, mapping slaves with // collision on this process std::vector> collision_slaves(indegree); std::vector>> collision_masters(indegree); - std::vector>> - collision_coeffs(indegree); + std::vector>> collision_coeffs( + indegree); std::vector>> collision_owners(indegree); std::vector>> @@ -1177,17 +1174,17 @@ mpc_data create_contact_inelastic_condition( std::vector remote_cell_collisions = dolfinx_mpc::find_local_collisions(*V->mesh(), bb_tree, recv_coords, eps2); - auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( + auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( *V, recv_coords, remote_cell_collisions); assert(basis_shape.back() == 1); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> basis_span(basis.data(), basis_shape[0], basis_shape[1]); // TODO: Rework this so it is the same as the code on owning process. // Preferably get rid of all the std::map's // Work arrays for loop std::vector r_master; - std::vector r_coeff; + std::vector r_coeff; std::vector remote_master_global; for (std::int32_t i = 0; i < indegree; ++i) { @@ -1207,7 +1204,7 @@ mpc_data create_contact_inelastic_condition( // Store block and non-zero basis values for (std::size_t k = 0; k < cell_blocks.size(); ++k) { - if (const PetscScalar c = basis_span(j, k); std::abs(c) > 1e-6) + if (const T c = basis_span(j, k); std::abs(c) > 1e-6) { r_coeff.push_back(c); r_master.push_back(cell_blocks[k]); @@ -1268,7 +1265,7 @@ mpc_data create_contact_inelastic_condition( std::vector offset_for_blocks; std::vector offsets_in_blocks; std::vector found_owners; - std::vector found_coefficients; + std::vector found_coefficients; for (std::int32_t i = 0; i < indegree; ++i) { std::int32_t master_offset = 0; @@ -1282,7 +1279,7 @@ mpc_data create_contact_inelastic_condition( found_masters.insert(found_masters.end(), masters_ij.begin(), masters_ij.end()); - std::vector& coeffs_ij = collision_coeffs[i][slave]; + std::vector& coeffs_ij = collision_coeffs[i][slave]; found_coefficients.insert(found_coefficients.end(), coeffs_ij.begin(), coeffs_ij.end()); @@ -1348,13 +1345,13 @@ mpc_data create_contact_inelastic_condition( remote_colliding_masters.data(), num_inc_masters.data(), disp_inc_masters.data(), dolfinx::MPI::mpi_type(), neighborhood_comms[1]); - std::vector remote_colliding_coeffs(disp_inc_masters.back()); - MPI_Neighbor_alltoallv( - found_coefficients.data(), num_collision_masters.data(), - send_disp_masters.data(), dolfinx::MPI::mpi_type(), - remote_colliding_coeffs.data(), num_inc_masters.data(), - disp_inc_masters.data(), dolfinx::MPI::mpi_type(), - neighborhood_comms[1]); + std::vector remote_colliding_coeffs(disp_inc_masters.back()); + MPI_Neighbor_alltoallv(found_coefficients.data(), + num_collision_masters.data(), send_disp_masters.data(), + dolfinx::MPI::mpi_type(), + remote_colliding_coeffs.data(), num_inc_masters.data(), + disp_inc_masters.data(), dolfinx::MPI::mpi_type(), + neighborhood_comms[1]); std::vector remote_colliding_owners(disp_inc_masters.back()); MPI_Neighbor_alltoallv( found_owners.data(), num_collision_masters.data(), @@ -1427,7 +1424,7 @@ mpc_data create_contact_inelastic_condition( local_masters[local_slave].insert(local_masters[local_slave].end(), masters_.begin(), masters_.end()); - std::vector coeffs_( + std::vector coeffs_( remote_colliding_coeffs.begin() + disp_inc_masters[i] + master_offsets[j] + block_offsets[k], remote_colliding_coeffs.begin() + disp_inc_masters[i] @@ -1494,7 +1491,7 @@ mpc_data create_contact_inelastic_condition( // (can include repeats of data) std::map> proc_to_ghost; std::map> proc_to_ghost_masters; - std::map> proc_to_ghost_coeffs; + std::map> proc_to_ghost_coeffs; std::map> proc_to_ghost_owners; std::map> proc_to_ghost_offsets; std::vector loc_block(1); @@ -1512,7 +1509,7 @@ mpc_data create_contact_inelastic_condition( { const std::int32_t slave = block * block_size + j; const std::vector& masters_i = local_masters[slave]; - const std::vector& coeffs_i = local_coeffs[slave]; + const std::vector& coeffs_i = local_coeffs[slave]; const std::vector& owners_i = local_owners[slave]; const auto num_masters = (std::int32_t)masters_i.size(); for (auto proc : shared_indices.links(slave / block_size)) @@ -1546,7 +1543,7 @@ mpc_data create_contact_inelastic_condition( // Flatten map of global slave ghost dofs to use alltoallv std::vector out_ghost_slaves; std::vector out_ghost_masters; - std::vector out_ghost_coeffs; + std::vector out_ghost_coeffs; std::vector out_ghost_owners; std::vector out_ghost_offsets; std::vector num_send_slaves(dest_ranks_ghost.size()); @@ -1617,13 +1614,12 @@ mpc_data create_contact_inelastic_condition( in_ghost_masters.data(), inc_num_masters.data(), disp_recv_ghost_masters.data(), dolfinx::MPI::mpi_type(), slave_to_ghost); - std::vector in_ghost_coeffs(disp_recv_ghost_masters.back()); + std::vector in_ghost_coeffs(disp_recv_ghost_masters.back()); MPI_Neighbor_alltoallv(out_ghost_coeffs.data(), num_send_masters.data(), disp_send_ghost_masters.data(), - dolfinx::MPI::mpi_type(), - in_ghost_coeffs.data(), inc_num_masters.data(), - disp_recv_ghost_masters.data(), - dolfinx::MPI::mpi_type(), slave_to_ghost); + dolfinx::MPI::mpi_type(), in_ghost_coeffs.data(), + inc_num_masters.data(), disp_recv_ghost_masters.data(), + dolfinx::MPI::mpi_type(), slave_to_ghost); std::vector in_ghost_owners(disp_recv_ghost_masters.back()); MPI_Neighbor_alltoallv( out_ghost_owners.data(), num_send_masters.data(), @@ -1652,7 +1648,7 @@ mpc_data create_contact_inelastic_condition( slaves.reserve(tdim * local_blocks.size()); std::vector masters; masters.reserve(slaves.size()); - std::vector coeffs_out; + std::vector coeffs_out; coeffs_out.reserve(slaves.size()); std::vector owners_out; owners_out.reserve(slaves.size()); diff --git a/cpp/PeriodicConstraint.h b/cpp/PeriodicConstraint.h index 36625725..32c866f4 100644 --- a/cpp/PeriodicConstraint.h +++ b/cpp/PeriodicConstraint.h @@ -25,8 +25,7 @@ template dolfinx_mpc::mpc_data _create_periodic_condition( const dolfinx::fem::FunctionSpace& V, std::span slave_blocks, - const std::function(std::span)>& relation, - T scale, + const std::function(std::span)>& relation, T scale, const std::function& parent_map, const dolfinx::fem::FunctionSpace& parent_space) { @@ -74,7 +73,7 @@ dolfinx_mpc::mpc_data _create_periodic_condition( // Tolerance for adding scaled basis values to MPC. Any scaled basis // value with lower absolute value than the tolerance is ignored - const double tol = 1e-13; + const U tol = 1e-13; auto mesh = V.mesh(); auto dofmap = V.dofmap(); @@ -100,23 +99,23 @@ dolfinx_mpc::mpc_data _create_periodic_condition( *V.mesh()->topology(), *V.dofmap(), local_blocks); // Compute relation(slave_blocks) - std::vector mapped_T_b(local_blocks.size() * 3); + std::vector mapped_T_b(local_blocks.size() * 3); const std::array T_shape = {local_blocks.size(), 3}; { - std::experimental::mdspan< - double, std::experimental::extents< - std::size_t, std::experimental::dynamic_extent, 3>> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> mapped_T(mapped_T_b.data(), T_shape); // Tabulate dof coordinates for each dof auto [x, x_shape] = dolfinx_mpc::tabulate_dof_coordinates( V, local_blocks, slave_cells, true); // Map all slave coordinates using the relation - std::vector mapped_x = relation(x); - std::experimental::mdspan< - double, std::experimental::extents> + std::vector mapped_x = relation(x); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>> x_(mapped_x.data(), x_shape); for (std::size_t i = 0; i < mapped_T.extent(0); ++i) for (std::size_t j = 0; j < mapped_T.extent(1); ++j) @@ -138,10 +137,10 @@ dolfinx_mpc::mpc_data _create_periodic_condition( std::vector local_cell_collisions = dolfinx_mpc::find_local_collisions(*mesh, tree, mapped_T_b, 1e-20); dolfinx::common::Timer t0("~~Periodic: Local cell and eval basis"); - auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( + auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( V, mapped_T_b, local_cell_collisions); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> tabulated_basis_values(basis_values.data(), basis_shape); t0.stop(); // Create output arrays @@ -285,7 +284,7 @@ dolfinx_mpc::mpc_data _create_periodic_condition( // Communicate coordinates of slave blocks std::vector out_placement(outdegree, 0); - std::vector coords_out(disp_out.back() * 3); + std::vector coords_out(disp_out.back() * 3); for (std::size_t i = 0; i < local_cell_collisions.size(); i++) { @@ -316,10 +315,11 @@ dolfinx_mpc::mpc_data _create_periodic_condition( // Communciate coordinates with other process const std::array cr_shape = {(std::size_t)disp_in.back(), 3}; - std::vector coords_recvb(cr_shape.front() * cr_shape.back()); - std::experimental::mdspan< - const double, std::experimental::extents< - std::size_t, std::experimental::dynamic_extent, 3>> + std::vector coords_recvb(cr_shape.front() * cr_shape.back()); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, + MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent, 3>> coords_recv(coords_recvb.data(), cr_shape); // Take into account that we send three values per slave @@ -330,11 +330,10 @@ dolfinx_mpc::mpc_data _create_periodic_condition( std::for_each(num_recv_slaves.begin(), num_recv_slaves.end(), m_3); // Communicate coordinates - MPI_Neighbor_alltoallv(coords_out.data(), num_out_slaves.data(), - disp_out.data(), dolfinx::MPI::mpi_type(), - coords_recvb.data(), num_recv_slaves.data(), - disp_in.data(), dolfinx::MPI::mpi_type(), - slave_to_master); + MPI_Neighbor_alltoallv( + coords_out.data(), num_out_slaves.data(), disp_out.data(), + dolfinx::MPI::mpi_type(), coords_recvb.data(), num_recv_slaves.data(), + disp_in.data(), dolfinx::MPI::mpi_type(), slave_to_master); // Reset in_displacements to be per block for later usage auto d_3 = [](auto& num) { num /= 3; }; @@ -358,10 +357,10 @@ dolfinx_mpc::mpc_data _create_periodic_condition( std::vector remote_cell_collisions = dolfinx_mpc::find_local_collisions(*mesh, tree, coords_recvb, 1e-20); auto [remote_basis_valuesb, r_basis_shape] - = dolfinx_mpc::evaluate_basis_functions(V, coords_recvb, - remote_cell_collisions); - std::experimental::mdspan> + = dolfinx_mpc::evaluate_basis_functions(V, coords_recvb, + remote_cell_collisions); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> remote_basis_values(remote_basis_valuesb.data(), r_basis_shape); // Find remote masters and count how many to send to each process @@ -496,12 +495,12 @@ template dolfinx_mpc::mpc_data geometrical_condition( const std::shared_ptr> V, const std::function( - std::experimental::mdspan< - const double, - std::experimental::extents>)>& + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)>& indicator, - const std::function(std::span)>& relation, + const std::function(std::span)>& relation, const std::vector>>& bcs, T scale, bool collapse) { @@ -564,7 +563,7 @@ dolfinx_mpc::mpc_data topological_condition( const std::shared_ptr> V, const std::shared_ptr> meshtag, const std::int32_t tag, - const std::function(std::span)>& relation, + const std::function(std::span)>& relation, const std::vector>>& bcs, T scale, bool collapse) { @@ -611,8 +610,8 @@ dolfinx_mpc::mpc_data topological_condition( // Identity map const auto sub_map = [](const std::int32_t& dof) { return dof; }; - return _create_periodic_condition(*V, std::span(reduced_blocks), - relation, scale, sub_map, *V); + return _create_periodic_condition(*V, std::span(reduced_blocks), + relation, scale, sub_map, *V); } }; @@ -624,35 +623,35 @@ namespace dolfinx_mpc mpc_data create_periodic_condition_geometrical( const std::shared_ptr> V, const std::function( - std::experimental::mdspan< + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const double, - std::experimental::extents>)>& - indicator, + MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)>& indicator, const std::function(std::span)>& relation, const std::vector>>& bcs, double scale, bool collapse) { - return impl::geometrical_condition(V, indicator, relation, bcs, scale, - collapse); + return impl::geometrical_condition(V, indicator, relation, + bcs, scale, collapse); } mpc_data> create_periodic_condition_geometrical( const std::shared_ptr> V, const std::function( - std::experimental::mdspan< + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const double, - std::experimental::extents>)>& - indicator, + MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)>& indicator, const std::function(std::span)>& relation, const std::vector< std::shared_ptr>>>& bcs, std::complex scale, bool collapse) { - return impl::geometrical_condition>( + return impl::geometrical_condition, double>( V, indicator, relation, bcs, scale, collapse); } @@ -665,8 +664,8 @@ mpc_data create_periodic_condition_topological( bcs, double scale, bool collapse) { - return impl::topological_condition(V, meshtag, tag, relation, bcs, - scale, collapse); + return impl::topological_condition(V, meshtag, tag, relation, + bcs, scale, collapse); } mpc_data> create_periodic_condition_topological( @@ -679,7 +678,70 @@ mpc_data> create_periodic_condition_topological( bcs, std::complex scale, bool collapse) { - return impl::topological_condition>( + return impl::topological_condition, double>( V, meshtag, tag, relation, bcs, scale, collapse); } + +mpc_data create_periodic_condition_geometrical( + const std::shared_ptr> V, + const std::function( + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const float, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)>& + indicator, + const std::function(std::span)>& relation, + const std::vector>>& + bcs, + float scale, bool collapse) +{ + return impl::geometrical_condition(V, indicator, relation, bcs, + scale, collapse); +} + +mpc_data> create_periodic_condition_geometrical( + const std::shared_ptr> V, + const std::function( + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const float, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)>& + indicator, + const std::function(std::span)>& relation, + const std::vector< + std::shared_ptr>>>& + bcs, + std::complex scale, bool collapse) +{ + return impl::geometrical_condition, float>( + V, indicator, relation, bcs, scale, collapse); +} + +mpc_data create_periodic_condition_topological( + const std::shared_ptr> V, + const std::shared_ptr> meshtag, + const std::int32_t tag, + const std::function(std::span)>& relation, + const std::vector>>& + bcs, + float scale, bool collapse) +{ + return impl::topological_condition(V, meshtag, tag, relation, + bcs, scale, collapse); +} + +mpc_data> create_periodic_condition_topological( + const std::shared_ptr> V, + const std::shared_ptr> meshtag, + const std::int32_t tag, + const std::function(std::span)>& relation, + const std::vector< + std::shared_ptr>>>& + bcs, + std::complex scale, bool collapse) +{ + return impl::topological_condition, float>( + V, meshtag, tag, relation, bcs, scale, collapse); +} + } // namespace dolfinx_mpc \ No newline at end of file diff --git a/cpp/SlipConstraint.h b/cpp/SlipConstraint.h index 0199bd4d..e6f7eb3e 100644 --- a/cpp/SlipConstraint.h +++ b/cpp/SlipConstraint.h @@ -12,27 +12,26 @@ namespace dolfinx_mpc { // -template -mpc_data create_slip_condition( +template +mpc_data create_slip_condition( std::shared_ptr>& space, const dolfinx::mesh::MeshTags& meshtags, std::int32_t marker, - const dolfinx::fem::Function& v, - std::vector>> - bcs, + const dolfinx::fem::Function& v, + std::vector>> bcs, const bool sub_space) { // Map from collapsed sub space to parent space std::function parent_map; // Interpolate input function into suitable space - std::shared_ptr> n; + std::shared_ptr> n; if (sub_space) { std::pair, std::vector> collapsed_space = space->collapse(); auto V_ptr = std::make_shared>( std::move(collapsed_space.first)); - n = std::make_shared>(V_ptr); + n = std::make_shared>(V_ptr); n->interpolate(v); parent_map = [sub_map = collapsed_space.second](const std::int32_t dof) { return sub_map[dof]; }; @@ -40,7 +39,7 @@ mpc_data create_slip_condition( else { parent_map = [](const std::int32_t dof) { return dof; }; - n = std::make_shared>(space); + n = std::make_shared>(space); n->interpolate(v); } @@ -74,7 +73,7 @@ mpc_data create_slip_condition( // Remove Dirichlet BC dofs const std::vector bc_marker - = dolfinx_mpc::is_bc(*space, entity_dofs[0], bcs); + = dolfinx_mpc::is_bc(*space, entity_dofs[0], bcs); num_normal_components = n->function_space()->dofmap()->index_map_bs(); slave_blocks.reserve(entity_dofs[0].size()); for (std::size_t i = 0; i < bc_marker.size(); i++) @@ -94,32 +93,32 @@ mpc_data create_slip_condition( meshtags.dim(), slave_facets); // Remove Dirichlet BC dofs const std::vector bc_marker - = dolfinx_mpc::is_bc(*space, all_slave_blocks, bcs); + = dolfinx_mpc::is_bc(*space, all_slave_blocks, bcs); for (std::size_t i = 0; i < bc_marker.size(); i++) if (!bc_marker[i]) slave_blocks.push_back(all_slave_blocks[i]); } - const std::span& n_vec = n->x()->array(); + const std::span& n_vec = n->x()->array(); // Arrays holding MPC data std::vector slaves; std::vector masters; - std::vector coeffs; + std::vector coeffs; std::vector owners; std::vector offsets(1, 0); // Temporary arrays used to hold information about masters std::vector pair_m; for (auto block : slave_blocks) { - std::span normal( + std::span normal( std::next(n_vec.begin(), block * num_normal_components), num_normal_components); // Determine slave dof by component with biggest normal vector (to avoid // issues with grids aligned with coordiante system) auto max_el = std::max_element(normal.begin(), normal.end(), - [](PetscScalar a, PetscScalar b) + [](T a, T b) { return std::norm(a) < std::norm(b); }); auto slave_index = std::distance(normal.begin(), max_el); assert(slave_index < num_normal_components); @@ -128,7 +127,7 @@ mpc_data create_slip_condition( slaves.push_back(parent_slave); std::vector parent_masters; - std::vector pair_c; + std::vector pair_c; std::vector pair_o; for (std::int32_t i = 0; i < num_normal_components; ++i) { @@ -136,7 +135,7 @@ mpc_data create_slip_condition( { const std::int32_t parent_dof = parent_map(block * num_normal_components + i); - PetscScalar coeff = -normal[i] / normal[slave_index]; + T coeff = -normal[i] / normal[slave_index]; parent_masters.push_back(parent_dof); pair_c.push_back(coeff); const std::int32_t m_rank @@ -164,7 +163,7 @@ mpc_data create_slip_condition( offsets.push_back((std::int32_t)masters.size()); owners.insert(owners.end(), pair_o.begin(), pair_o.end()); } - mpc_data data; + mpc_data data; data.slaves = slaves; data.masters = masters; diff --git a/cpp/assemble_matrix.cpp b/cpp/assemble_matrix.cpp index 093fb7f7..a336aba1 100644 --- a/cpp/assemble_matrix.cpp +++ b/cpp/assemble_matrix.cpp @@ -11,9 +11,9 @@ #include #include -namespace stdex = std::experimental; -using mdspan2_t - = stdex::mdspan>; +using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; namespace { @@ -32,9 +32,11 @@ namespace /// @returns The matrix stripped of slave contributions template void fill_stripped_matrix( - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> Ae_stripped, - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> Ae, const std::array& num_dofs, const std::array& bs, @@ -100,7 +102,8 @@ void modify_mpc_cell( std::span, std::span)>& mat_set, const std::array& num_dofs, - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> Ae, const std::array, 2>& dofs, const std::array& bs, @@ -141,13 +144,15 @@ void modify_mpc_cell( >= std::size_t(2 * ndim0 * ndim1 + ndim0 + ndim1)); std::fill(scratch_memory.begin(), scratch_memory.end(), T(0)); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> Ae_original(scratch_memory.data(), ndim0, ndim1); // Copy Ae into new matrix for distirbution of master dofs std::copy_n(Ae.data_handle(), ndim0 * ndim1, Ae_original.data_handle()); - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> Ae_stripped(std::next(scratch_memory.data(), ndim0 * ndim1), ndim0, ndim1); // Build matrix where all slave-slave entries are 0 for usage to row and @@ -267,8 +272,7 @@ void assemble_exterior_facets( const std::function, std::span, const std::span)>& mat_add_values, - const dolfinx::mesh::Mesh& mesh, - std::span facets, + const dolfinx::mesh::Mesh& mesh, std::span facets, const std::function&, const std::span&, std::int32_t, int)>& apply_dof_transformation, @@ -278,7 +282,7 @@ void assemble_exterior_facets( std::int32_t, int)>& apply_dof_transformation_to_transpose, const dolfinx::fem::DofMap& dofmap1, const std::vector& bc0, const std::vector& bc1, - const std::function& kernel, const std::span coeffs, int cstride, const std::vector& constants, @@ -300,8 +304,9 @@ void assemble_exterior_facets( cell_to_slaves = {mpc0->cell_to_slaves(), mpc1->cell_to_slaves()}; // Get mesh data - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> x_dofmap = mesh.geometry().dofmap(); const int num_dofs_g = x_dofmap.extent(1); @@ -319,8 +324,9 @@ void assemble_exterior_facets( const std::array bs = {bs0, bs1}; std::vector Aeb(ndim0 * ndim1); - std::experimental::mdspan> Ae( - Aeb.data(), ndim0, ndim1); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + Ae(Aeb.data(), ndim0, ndim1); const std::span _Ae(Aeb); std::vector scratch_memory(2 * ndim0 * ndim1 + ndim0 + ndim1); @@ -331,7 +337,9 @@ void assemble_exterior_facets( const int local_facet = facets[l + 1]; // Get cell vertex coordinates - auto x_dofs = stdex::submdspan(x_dofmap, cell, stdex::full_extent); + namespace stdex = std::experimental; + auto x_dofs = stdex::submdspan(x_dofmap, cell, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) { std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3, @@ -402,7 +410,7 @@ void assemble_cells_impl( const std::function, std::span, const std::span)>& mat_add_values, - const dolfinx::mesh::Geometry& geometry, + const dolfinx::mesh::Geometry& geometry, std::span active_cells, std::function, const std::span, const std::int32_t, const int)> @@ -413,7 +421,7 @@ void assemble_cells_impl( apply_dof_transformation_to_transpose, const dolfinx::fem::DofMap& dofmap1, const std::vector& bc0, const std::vector& bc1, - const std::function& kernel, const std::span& coeffs, int cstride, const std::vector& constants, @@ -436,9 +444,9 @@ void assemble_cells_impl( cell_to_slaves = {mpc0->cell_to_slaves(), mpc1->cell_to_slaves()}; // Prepare cell geometry - namespace stdex = std::experimental; - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = x_dofmap.extent(1); std::span x_g = geometry.x(); @@ -453,16 +461,18 @@ void assemble_cells_impl( const std::array num_dofs = {num_dofs0, num_dofs1}; std::vector Aeb(ndim0 * ndim1); - std::experimental::mdspan> Ae( - Aeb.data(), ndim0, ndim1); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + Ae(Aeb.data(), ndim0, ndim1); const std::span _Ae(Aeb); std::vector scratch_memory(2 * ndim0 * ndim1 + ndim0 + ndim1); - + namespace stdex = std::experimental; for (std::size_t c = 0; c < active_cells.size(); c++) { const std::int32_t cell = active_cells[c]; // Get cell coordinates/geometry - auto x_dofs = stdex::submdspan(x_dofmap, cell, stdex::full_extent); + auto x_dofs = stdex::submdspan(x_dofmap, cell, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) { std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3, @@ -721,3 +731,44 @@ void dolfinx_mpc::assemble_matrix( { _assemble_matrix(mat_add_block, mat_add, a, mpc0, mpc1, bcs, diagval); } +//----------------------------------------------------------------------------- +void dolfinx_mpc::assemble_matrix( + const std::function, + std::span, + const std::span&)>& mat_add_block, + const std::function, + std::span, + const std::span&)>& mat_add, + const dolfinx::fem::Form& a, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc0, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc1, + const std::vector>>& + bcs, + const float diagval) +{ + _assemble_matrix(mat_add_block, mat_add, a, mpc0, mpc1, bcs, diagval); +} +//----------------------------------------------------------------------------- +void dolfinx_mpc::assemble_matrix( + const std::function< + int(std::span, std::span, + const std::span>&)>& mat_add_block, + const std::function< + int(std::span, std::span, + const std::span>&)>& mat_add, + const dolfinx::fem::Form>& a, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc0, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc1, + const std::vector< + std::shared_ptr>>>& + bcs, + const std::complex diagval) +{ + _assemble_matrix(mat_add_block, mat_add, a, mpc0, mpc1, bcs, diagval); +} diff --git a/cpp/assemble_matrix.h b/cpp/assemble_matrix.h index a68d1317..baf90212 100644 --- a/cpp/assemble_matrix.h +++ b/cpp/assemble_matrix.h @@ -69,4 +69,58 @@ void assemble_matrix( bcs, const std::complex diagval = 1.0); +/// Assemble bilinear form into a matrix +/// @param[in] mat_add_block The function for adding block values into the +/// matrix +/// @param[in] mat_add The function for adding values into the matrix +/// @param[in] a The bilinear from to assemble +/// @param[in] bcs Boundary conditions to apply. For boundary condition +/// dofs the row and column are zeroed. The diagonal entry is not set. +/// @param[in] diagval Value to set on diagonal of matrix for slave dofs and +/// Dirichlet BC (default=1) +void assemble_matrix( + const std::function, + std::span, + const std::span&)>& mat_add_block, + const std::function, + std::span, + const std::span&)>& mat_add, + const dolfinx::fem::Form& a, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc0, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc1, + const std::vector>>& + bcs, + const float diagval = 1.0); + +//----------------------------------------------------------------------------- +/// Assemble bilinear form into a matrix +/// @param[in] mat_add_block The function for adding block values into the +/// matrix +/// @param[in] mat_add The function for adding values into the matrix +/// @param[in] a The bilinear from to assemble +/// @param[in] bcs Boundary conditions to apply. For boundary condition +/// dofs the row and column are zeroed. The diagonal entry is not set. +/// @param[in] diagval Value to set on diagonal of matrix for slave dofs and +/// Dirichlet BC (default=1) +void assemble_matrix( + const std::function< + int(std::span, std::span, + const std::span>&)>& mat_add_block, + const std::function< + int(std::span, std::span, + const std::span>&)>& mat_add, + const dolfinx::fem::Form>& a, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc0, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc1, + const std::vector< + std::shared_ptr>>>& + bcs, + const std::complex diagval = 1.0); + } // namespace dolfinx_mpc \ No newline at end of file diff --git a/cpp/assemble_vector.cpp b/cpp/assemble_vector.cpp index 32016e8c..a762272b 100644 --- a/cpp/assemble_vector.cpp +++ b/cpp/assemble_vector.cpp @@ -11,9 +11,9 @@ #include #include -namespace stdex = std::experimental; -using mdspan2_t - = stdex::mdspan>; +using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; namespace { @@ -30,12 +30,11 @@ namespace /// assembles into a local element matrix for a given entity /// @tparam T Scalar type for vector /// @tparam e stride Stride for each entity in active_entities -template +template void _assemble_entities_impl( std::span b, std::span active_entities, const dolfinx::fem::DofMap& dofmap, - const std::shared_ptr>& - mpc, + const std::shared_ptr>& mpc, const std::function)> fetch_cells, const std::function, std::span, @@ -93,7 +92,7 @@ void _assemble_vector( std::span b, const dolfinx::fem::Form& L, const std::shared_ptr>& mpc) { - + namespace stdex = std::experimental; const auto mesh = L.mesh(); assert(mesh); @@ -109,11 +108,11 @@ void _assemble_vector( auto coefficients = dolfinx::fem::make_coefficients_span(coeff_vec); // Prepare cell geometry - namespace stdex = std::experimental; - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> x_dofmap = mesh->geometry().dofmap(); - std::span x_g = mesh->geometry().x(); + std::span x_g = mesh->geometry().x(); // Prepare dof tranformation data auto element = L.function_spaces().at(0)->element(); @@ -130,7 +129,7 @@ void _assemble_vector( cell_info = std::span(mesh->topology()->get_cell_permutation_info()); } const std::size_t num_dofs_g = x_dofmap.extent(1); - std::vector coordinate_dofs(3 * num_dofs_g); + std::vector coordinate_dofs(3 * num_dofs_g); if (L.num_integrals(dolfinx::fem::IntegralType::cell) > 0) { @@ -153,7 +152,8 @@ void _assemble_vector( auto cell = entity.front(); // Fetch the coordinates of the cell - auto x_dofs = stdex::submdspan(x_dofmap, cell, stdex::full_extent); + auto x_dofs = stdex::submdspan( + x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) { std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3, @@ -172,8 +172,8 @@ void _assemble_vector( // Assemble over all active cells std::span active_cells = L.domain(dolfinx::fem::IntegralType::cell, i); - _assemble_entities_impl(b, active_cells, *dofmap, mpc, fetch_cell, - assemble_local_cell_vector); + _assemble_entities_impl(b, active_cells, *dofmap, mpc, + fetch_cell, assemble_local_cell_vector); } } // Prepare permutations for exterior and interior facet integrals @@ -201,7 +201,8 @@ void _assemble_vector( // Fetch the coordinates of the cell const std::int32_t cell = entity[0]; const int local_facet = entity[1]; - auto x_dofs = stdex::submdspan(x_dofmap, cell, stdex::full_extent); + auto x_dofs = stdex::submdspan( + x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) { std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3, @@ -220,8 +221,9 @@ void _assemble_vector( // Assemble over all active cells std::span active_facets = L.domain(dolfinx::fem::IntegralType::exterior_facet, i); - _assemble_entities_impl(b, active_facets, *dofmap, mpc, fetch_cell, - assemble_local_exterior_facet_vector); + _assemble_entities_impl(b, active_facets, *dofmap, mpc, + fetch_cell, + assemble_local_exterior_facet_vector); } } @@ -263,4 +265,22 @@ void dolfinx_mpc::assemble_vector( { _assemble_vector>(b, L, mpc); } + +void dolfinx_mpc::assemble_vector( + std::span b, const dolfinx::fem::Form& L, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc) +{ + _assemble_vector(b, L, mpc); +} + +void dolfinx_mpc::assemble_vector( + std::span> b, + const dolfinx::fem::Form>& L, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc) +{ + _assemble_vector>(b, L, mpc); +} //----------------------------------------------------------------------------- diff --git a/cpp/assemble_vector.h b/cpp/assemble_vector.h index 6a411a0e..6e906a3b 100644 --- a/cpp/assemble_vector.h +++ b/cpp/assemble_vector.h @@ -11,9 +11,9 @@ #include #include #include -namespace stdex = std::experimental; -using mdspan2_t - = stdex::mdspan>; +using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; namespace dolfinx_mpc { @@ -90,4 +90,26 @@ void assemble_vector( const dolfinx_mpc::MultiPointConstraint, double>>& mpc); +/// Assemble a linear form into a vector +/// @param[in] b The vector to be assembled. It will not be zeroed before +/// assembly. +/// @param[in] L The linear forms to assemble into b +/// @param[in] mpc The multi-point constraint +void assemble_vector( + std::span b, const dolfinx::fem::Form& L, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc); + +/// Assemble a linear form into a vector +/// @param[in] b The vector to be assembled. It will not be zeroed before +/// assembly. +/// @param[in] L The linear forms to assemble into b +/// @param[in] mpc The multi-point constraint +void assemble_vector( + std::span> b, + const dolfinx::fem::Form>& L, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc); + } // namespace dolfinx_mpc diff --git a/cpp/lifting.h b/cpp/lifting.h index 8257574d..767a661d 100644 --- a/cpp/lifting.h +++ b/cpp/lifting.h @@ -140,7 +140,7 @@ template void apply_lifting( std::span b, const std::shared_ptr> a, const std::vector>>& bcs, - const std::span& x0, double scale, + const std::span& x0, T scale, const std::shared_ptr>& mpc1) { const std::vector constants = pack_constants(*a); @@ -183,13 +183,14 @@ void apply_lifting( // Prepare cell geometry namespace stdex = std::experimental; - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> x_dofmap = mesh->geometry().dofmap(); - std::span x_g = mesh->geometry().x(); + std::span x_g = mesh->geometry().x(); const int tdim = mesh->topology()->dim(); const std::size_t num_dofs_g = x_dofmap.extent(1); - std::vector coordinate_dofs(3 * num_dofs_g); + std::vector coordinate_dofs(3 * num_dofs_g); std::span cell_info; if (needs_transformation_data) @@ -229,7 +230,8 @@ void apply_lifting( auto cell = entity.front(); // Fetch the coordinates of the cell - auto x_dofs = stdex::submdspan(x_dofmap, cell, stdex::full_extent); + auto x_dofs = stdex::submdspan( + x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) { std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3, @@ -312,7 +314,8 @@ void apply_lifting( const int local_facet = entity[1]; // Fetch the coordinates of the cell - auto x_dofs = stdex::submdspan(x_dofmap, cell, stdex::full_extent); + auto x_dofs = stdex::submdspan( + x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) { std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3, @@ -498,4 +501,117 @@ void apply_lifting( } } } + +/// Modify b such that: +/// +/// b <- b - scale * K^T (A_j (g_j 0 x0_j)) +/// +/// where j is a block (nest) row index and K^T is the reduction matrix stemming +/// from the multi point constraint. For non - blocked problems j = 0. +/// The boundary conditions bcs1 are on the trial spaces V_j. +/// The forms in [a] must have the same test space as L (from +/// which b was built), but the trial space may differ. If x0 is not +/// supplied, then it is treated as 0. +/// @param[in,out] b The vector to be modified +/// @param[in] a The bilinear formss, where a[j] is the form that +/// generates A[j] +/// @param[in] bcs List of boundary conditions for each block, i.e. bcs1[2] +/// are the boundary conditions applied to the columns of a[2] / x0[2] +/// block. +/// @param[in] x0 The vectors used in the lifitng. +/// @param[in] scale Scaling to apply +/// @param[in] mpc The multi point constraints +void apply_lifting( + std::span b, + const std::vector>> a, + const std::vector< + std::vector>>>& + bcs1, + const std::vector>& x0, float scale, + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint>& mpc) +{ + if (!x0.empty() and x0.size() != a.size()) + { + throw std::runtime_error( + "Mismatch in size between x0 and bilinear form in assembler."); + } + + if (a.size() != bcs1.size()) + { + throw std::runtime_error( + "Mismatch in size between a and bcs in assembler."); + } + for (std::size_t j = 0; j < a.size(); ++j) + { + if (x0.empty()) + { + impl::apply_lifting(b, a[j], bcs1[j], std::span(), + scale, mpc); + } + else + { + impl::apply_lifting(b, a[j], bcs1[j], x0[j], scale, mpc); + } + } +} +/// Modify b such that: +/// +/// b <- b - scale * K^T (A_j (g_j 0 x0_j)) +/// +/// where j is a block (nest) row index and K^T is the reduction matrix stemming +/// from the multi point constraint. For non - blocked problems j = 0. +/// The boundary conditions bcs1 are on the trial spaces V_j. +/// The forms in [a] must have the same test space as L (from +/// which b was built), but the trial space may differ. If x0 is not +/// supplied, then it is treated as 0. +/// @param[in,out] b The vector to be modified +/// @param[in] a The bilinear formss, where a[j] is the form that +/// generates A[j] +/// @param[in] bcs List of boundary conditions for each block, i.e. bcs1[2] +/// are the boundary conditions applied to the columns of a[2] / x0[2] +/// block. +/// @param[in] x0 The vectors used in the lifitng. +/// @param[in] scale Scaling to apply +/// @param[in] mpc The multi point constraints +void apply_lifting( + std::span> b, + const std::vector< + std::shared_ptr>>> + a, + const std::vector>>>>& + bcs1, + const std::vector>>& x0, float scale, + + const std::shared_ptr< + const dolfinx_mpc::MultiPointConstraint, float>>& + mpc) +{ + if (!x0.empty() and x0.size() != a.size()) + { + throw std::runtime_error( + "Mismatch in size between x0 and bilinear form in assembler."); + } + + if (a.size() != bcs1.size()) + { + throw std::runtime_error( + "Mismatch in size between a and bcs in assembler."); + } + for (std::size_t j = 0; j < a.size(); ++j) + { + if (x0.empty()) + { + impl::apply_lifting>( + b, a[j], bcs1[j], std::span>(), scale, mpc); + } + else + { + impl::apply_lifting>(b, a[j], bcs1[j], x0[j], scale, + mpc); + } + } +} + } // namespace dolfinx_mpc \ No newline at end of file diff --git a/cpp/mpc_helpers.h b/cpp/mpc_helpers.h index a5adb10f..afbc29ba 100644 --- a/cpp/mpc_helpers.h +++ b/cpp/mpc_helpers.h @@ -209,9 +209,10 @@ dolfinx::fem::FunctionSpace create_extended_functionspace( } // Extract information from the old dofmap to create a new one - namespace stdex = std::experimental; - stdex::mdspan> dofmap_adj - = old_dofmap.map(); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + dofmap_adj = old_dofmap.map(); // Copy dofmap std::vector flattened_dofmap; flattened_dofmap.reserve(dofmap_adj.extent(0) * dofmap_adj.extent(1)); diff --git a/cpp/utils.cpp b/cpp/utils.cpp index c46986f0..290d1d97 100644 --- a/cpp/utils.cpp +++ b/cpp/utils.cpp @@ -14,42 +14,6 @@ using namespace dolfinx_mpc; -//----------------------------------------------------------------------------- -dolfinx::la::petsc::Matrix dolfinx_mpc::create_matrix( - const dolfinx::fem::Form& a, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc0, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc1, - const std::string& type) -{ - dolfinx::common::Timer timer("~MPC: Create Matrix"); - - // Build sparsitypattern - dolfinx::la::SparsityPattern pattern = create_sparsity_pattern(a, mpc0, mpc1); - - // Finalise communication - dolfinx::common::Timer timer_s("~MPC: Assemble sparsity pattern"); - pattern.finalize(); - timer_s.stop(); - - // Initialize matrix - dolfinx::la::petsc::Matrix A(a.mesh()->comm(), pattern, type); - - return A; -} -//----------------------------------------------------------------------------- -dolfinx::la::petsc::Matrix dolfinx_mpc::create_matrix( - const dolfinx::fem::Form& a, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc, - const std::string& type) -{ - return dolfinx_mpc::create_matrix(a, mpc, mpc, type); -} //----------------------------------------------------------------------------- std::array dolfinx_mpc::create_neighborhood_comms( MPI_Comm comm, const dolfinx::mesh::MeshTags& meshtags, diff --git a/cpp/utils.h b/cpp/utils.h index 1aacfe4a..b5101624 100644 --- a/cpp/utils.h +++ b/cpp/utils.h @@ -133,22 +133,41 @@ compute_shared_indices(std::shared_ptr> V) return V->dofmap()->index_map->index_to_dest_ranks(); } -dolfinx::la::petsc::Matrix -create_matrix(const dolfinx::fem::Form& a, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc0, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc1, - const std::string& type = std::string()); - -dolfinx::la::petsc::Matrix -create_matrix(const dolfinx::fem::Form& a, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc, - const std::string& type = std::string()); +template +dolfinx::la::petsc::Matrix create_matrix( + const dolfinx::fem::Form& a, + const std::shared_ptr> + mpc0, + const std::shared_ptr> + mpc1, + const std::string& type = std::string()) +{ + dolfinx::common::Timer timer("~MPC: Create Matrix"); + + // Build sparsitypattern + dolfinx::la::SparsityPattern pattern = create_sparsity_pattern(a, mpc0, mpc1); + + // Finalise communication + dolfinx::common::Timer timer_s("~MPC: Assemble sparsity pattern"); + pattern.finalize(); + timer_s.stop(); + + // Initialize matrix + dolfinx::la::petsc::Matrix A(a.mesh()->comm(), pattern, type); + + return A; +} + +template +dolfinx::la::petsc::Matrix create_matrix( + const dolfinx::fem::Form& a, + const std::shared_ptr> + mpc, + const std::string& type = std::string()) +{ + return dolfinx_mpc::create_matrix(a, mpc, mpc, type); +} + /// Create neighborhood communicators from every processor with a slave dof on /// it, to the processors with a set of master facets. /// @param[in] comm The MPI communicator to base communications on @@ -172,6 +191,7 @@ MPI_Comm create_owner_to_ghost_comm( /// Creates a normal approximation for the dofs in the closure of the attached /// facets, where the normal is an average if a dof belongs to multiple facets +/// FIXME: Remove petsc dependency here template dolfinx::fem::Function create_normal_approximation(std::shared_ptr> V, @@ -195,8 +215,8 @@ create_normal_approximation(std::shared_ptr> V, std::span _n(array, n); const std::int32_t bs = V->dofmap()->index_map_bs(); - std::array normal; - std::array n_0; + std::array normal; + std::array n_0; for (std::int32_t i = 0; i < block_to_entities.num_nodes(); i++) { auto ents = block_to_entities.links(i); @@ -209,7 +229,7 @@ create_normal_approximation(std::shared_ptr> V, for (std::size_t j = 1; j < normals.size() / 3; ++j) { // Align direction of normal vectors n_0 and n_j - double n_nj = std::transform_reduce( + U n_nj = std::transform_reduce( n_0.begin(), n_0.end(), std::next(normals.begin(), 3 * j), 0., std::plus{}, [](auto x, auto y) { return x * y; }); auto sign = n_nj / std::abs(n_nj); @@ -230,7 +250,7 @@ create_normal_approximation(std::shared_ptr> V, PetscScalar acc = 0; for (std::int32_t j = 0; j < bs; j++) acc += _n[i * bs + j] * _n[i * bs + j]; - if (double abs = std::abs(acc); abs > 1e-10) + if (U abs = std::abs(acc); abs > 1e-10) { for (std::int32_t j = 0; j < bs; j++) _n[i * bs + j] /= abs; @@ -341,13 +361,11 @@ create_block_to_cell_map(const dolfinx::mesh::Topology& topology, /// matrix. /// @param[in] mpc1 The multi point constraint to apply to the columns of the /// matrix. -template +template dolfinx::la::SparsityPattern create_sparsity_pattern( - const dolfinx::fem::Form& a, - const std::shared_ptr> - mpc0, - const std::shared_ptr> - mpc1) + const dolfinx::fem::Form& a, + const std::shared_ptr> mpc0, + const std::shared_ptr> mpc1) { { LOG(INFO) << "Generating MPC sparsity pattern"; @@ -376,17 +394,14 @@ dolfinx::la::SparsityPattern create_sparsity_pattern( LOG(INFO) << "Build standard pattern\n"; /// Create and build sparsity pattern for original form. Should be /// equivalent to calling create_sparsity_pattern(Form a) - build_standard_pattern(pattern, a); + build_standard_pattern(pattern, a); LOG(INFO) << "Build new pattern\n"; // Arrays replacing slave dof with master dof in sparsity pattern auto pattern_populator = [](dolfinx::la::SparsityPattern& pattern, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> - mpc, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint> + const std::shared_ptr> mpc, + const std::shared_ptr> mpc_off_axis, const auto& pattern_inserter, const auto& master_inserter) { @@ -926,9 +941,9 @@ dolfinx_mpc::mpc_data distribute_ghost_data( /// @returns basis values (not unrolled for block size) for each point. shape /// (num_points, number_of_dofs, value_size). Flattened row major template -std::pair, std::array> +std::pair, std::array> evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, - std::span x, + std::span x, std::span cells) { assert(x.size() % 3 == 0); @@ -948,8 +963,9 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, // Get geometry data namespace stdex = std::experimental; - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> x_dofmap = mesh->geometry().dofmap(); auto cmaps = mesh->geometry().cmaps(); @@ -987,7 +1003,7 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, assert(basis_shape[3] == std::size_t(element->value_size() / bs_element)); std::array reference_shape = {basis_shape[1], basis_shape[2], basis_shape[3]}; - std::vector output_basis(std::reduce( + std::vector output_basis(std::reduce( reference_shape.begin(), reference_shape.end(), 1, std::multiplies{})); if (num_points == 0) @@ -1000,10 +1016,12 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, cell_info = std::span(mesh->topology()->get_cell_permutation_info()); } - using cmdspan4_t - = stdex::mdspan>; - using mdspan2_t = stdex::mdspan>; - using mdspan3_t = stdex::mdspan>; + using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using mdspan3_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; // Create buffer for coordinate dofs and point in physical space std::vector coord_dofs_b(num_dofs_g * gdim); @@ -1019,16 +1037,16 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, cmdspan4_t phi0(phi0_b.data(), phi0_shape); cmaps[0].tabulate(1, std::vector(tdim, 0), {1, tdim}, phi0_b); auto dphi0 = stdex::submdspan(phi0, std::pair(1, tdim + 1), 0, - stdex::full_extent, 0); + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); // Data structure for evaluating geometry basis at specific points. // Used in non-affine case. std::array phi_shape = cmaps[0].tabulate_shape(1, 1); - std::vector phi_b( + std::vector phi_b( std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{})); cmdspan4_t phi(phi_b.data(), phi_shape); - auto dphi - = stdex::submdspan(phi, std::pair(1, tdim + 1), 0, stdex::full_extent, 0); + auto dphi = stdex::submdspan(phi, std::pair(1, tdim + 1), 0, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); // Reference coordinates for each point std::vector Xb(num_points * tdim); @@ -1037,10 +1055,10 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, // Geometry data at each point std::vector J_b(num_points * gdim * tdim); mdspan3_t J(J_b.data(), num_points, gdim, tdim); - std::vector K_b(num_points * tdim * gdim); + std::vector K_b(num_points * tdim * gdim); mdspan3_t K(K_b.data(), num_points, tdim, gdim); - std::vector detJ(num_points); - std::vector det_scratch(2 * gdim * tdim); + std::vector detJ(num_points); + std::vector det_scratch(2 * gdim * tdim); // Prepare geometry data in each cell for (std::size_t p = 0; p < cells.size(); ++p) @@ -1052,7 +1070,8 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, continue; // Get cell geometry (coordinate dofs) - auto x_dofs = stdex::submdspan(x_dofmap, cell_index, stdex::full_extent); + auto x_dofs = stdex::submdspan(x_dofmap, cell_index, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < num_dofs_g; ++i) { const int pos = 3 * x_dofs[i]; @@ -1063,11 +1082,17 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, for (std::size_t j = 0; j < gdim; ++j) xp(0, j) = x[3 * p + j]; - auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent); - auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent); - - std::array Xpb = {0, 0, 0}; - stdex::mdspan> + auto _J + = stdex::submdspan(J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _K + = stdex::submdspan(K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + + std::array Xpb = {0, 0, 0}; + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 1, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>> Xp(Xpb.data(), 1, tdim); // Compute reference coordinates X, and J, detJ and K @@ -1076,7 +1101,7 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, dolfinx::fem::CoordinateElement::compute_jacobian(dphi0, coord_dofs, _J); dolfinx::fem::CoordinateElement::compute_jacobian_inverse(_J, _K); - std::array x0 = {0, 0, 0}; + std::array x0 = {0, 0, 0}; for (std::size_t i = 0; i < coord_dofs.extent(1); ++i) x0[i] += coord_dofs(0, i); dolfinx::fem::CoordinateElement::pull_back_affine(Xp, _K, x0, xp); @@ -1112,10 +1137,14 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, std::vector basis_valuesb(num_basis_values); mdspan2_t basis_values(basis_valuesb.data(), basis_shape[2], basis_shape[3]); - using xu_t = stdex::mdspan>; - using xU_t = stdex::mdspan>; - using xJ_t = stdex::mdspan>; - using xK_t = stdex::mdspan>; + using xu_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using xU_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using xJ_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; + using xK_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; auto push_forward_fn = element->basix_element().template map_fn(); @@ -1136,10 +1165,15 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, apply_dof_transformation(basis_valuesb, cell_info, cell_index, (int)reference_value_size); - auto _U = stdex::submdspan(full_basis, p, stdex::full_extent, - stdex::full_extent); - auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent); - auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent); + auto _U = stdex::submdspan(full_basis, p, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _J + = stdex::submdspan(J, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); + auto _K + = stdex::submdspan(K, p, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); push_forward_fn(_U, basis_values, _J, detJ[p], _K); } return {output_basis, reference_shape}; @@ -1206,9 +1240,9 @@ std::pair, std::array> tabulate_dof_coordinates( } // Prepare cell geometry - namespace stdex = std::experimental; - std::experimental::mdspan> + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> x_dofmap = mesh->geometry().dofmap(); std::span x_g = mesh->geometry().x(); const std::size_t num_dofs_g = x_dofmap.extent(1); @@ -1221,7 +1255,8 @@ std::pair, std::array> tabulate_dof_coordinates( std::vector coordsb(std::reduce(coord_shape.cbegin(), coord_shape.cend(), 1, std::multiplies{})); - using mdspan2_t = stdex::mdspan>; + using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>; // Loop over cells and tabulate dofs assert(space_dimension == X_shape[0]); @@ -1247,10 +1282,12 @@ std::pair, std::array> tabulate_dof_coordinates( std::vector phi_b( std::reduce(bsize.begin(), bsize.end(), 1, std::multiplies{})); cmaps[0].tabulate(0, X_b, X_shape, phi_b); - stdex::mdspan> phi_full(phi_b.data(), - bsize); - auto phi = stdex::submdspan(phi_full, 0, stdex::full_extent, - stdex::full_extent, 0); + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + phi_full(phi_b.data(), bsize); + auto phi = stdex::submdspan(phi_full, 0, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0); // Create insertion function std::function inserter; @@ -1270,7 +1307,8 @@ std::pair, std::array> tabulate_dof_coordinates( for (std::size_t c = 0; c < cells.size(); ++c) { // Fetch the coordinates of the cell - auto x_dofs = stdex::submdspan(x_dofmap, cells[c], stdex::full_extent); + auto x_dofs = stdex::submdspan(x_dofmap, cells[c], + MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent); for (std::size_t i = 0; i < num_dofs_g; ++i) { const int pos = 3 * x_dofs[i]; diff --git a/python/benchmarks/bench_elasticity_edge.py b/python/benchmarks/bench_elasticity_edge.py index 5a580592..c617c022 100644 --- a/python/benchmarks/bench_elasticity_edge.py +++ b/python/benchmarks/bench_elasticity_edge.py @@ -67,7 +67,8 @@ def periodic_relation(x): periodic_mt = meshtags(mesh, edim, edges[arg_sort], np.full(len(edges), 2, dtype=np.int32)) mpc = MultiPointConstraint(V) - mpc.create_periodic_constraint_topological(V, periodic_mt, 2, periodic_relation, bcs, scale=0.5) + mpc.create_periodic_constraint_topological(V, periodic_mt, 2, periodic_relation, bcs, + scale=PETSc.ScalarType(0.5)) mpc.finalize() # Create traction meshtag diff --git a/python/demos/demo_elasticity_disconnect_2D.py b/python/demos/demo_elasticity_disconnect_2D.py index d274d1cf..8b8732fd 100644 --- a/python/demos/demo_elasticity_disconnect_2D.py +++ b/python/demos/demo_elasticity_disconnect_2D.py @@ -10,20 +10,19 @@ import gmsh import numpy as np -from dolfinx.fem import (Constant, Function, FunctionSpace, - VectorFunctionSpace, dirichletbc, - locate_dofs_geometrical, locate_dofs_topological) +from dolfinx.fem import (Constant, Function, VectorFunctionSpace, dirichletbc, + functionspace, locate_dofs_geometrical, + locate_dofs_topological, FunctionSpaceBase) from dolfinx.io import XDMFFile, gmshio from dolfinx.mesh import locate_entities_boundary +from dolfinx_mpc import LinearProblem, MultiPointConstraint +from dolfinx_mpc.utils import (create_point_to_point_constraint, + rigid_motions_nullspace) from mpi4py import MPI from petsc4py import PETSc from ufl import (Identity, Measure, SpatialCoordinate, TestFunction, TrialFunction, grad, inner, sym, tr) -from dolfinx_mpc import LinearProblem, MultiPointConstraint -from dolfinx_mpc.utils import (create_point_to_point_constraint, - rigid_motions_nullspace) - # Mesh parameters for creating a mesh consisting of two disjoint rectangles right_tag = 1 left_tag = 2 @@ -63,7 +62,7 @@ fdim = tdim - 1 # Locate cells with different elasticity parameters -DG0 = FunctionSpace(mesh, ("DG", 0)) +DG0 = functionspace(mesh, ("DG", 0)) left_dofs = locate_dofs_topological(DG0, tdim, ct.find(left_tag)) right_dofs = locate_dofs_topological(DG0, tdim, ct.find(right_tag)) @@ -108,7 +107,7 @@ def sigma(v): bcs = [bc_push, bc_fix] -def gather_dof_coordinates(V: FunctionSpace, dofs: np.ndarray): +def gather_dof_coordinates(V: FunctionSpaceBase, dofs: np.ndarray): """ Distributes the dof coordinates of this subset of dofs to all processors """ diff --git a/python/demos/demo_periodic3d_topological.py b/python/demos/demo_periodic3d_topological.py index 5025c8c8..298c4c4b 100644 --- a/python/demos/demo_periodic3d_topological.py +++ b/python/demos/demo_periodic3d_topological.py @@ -77,7 +77,7 @@ def periodic_relation(x): return out_x with Timer("~~Periodic: Compute mpc condition"): mpc = dolfinx_mpc.MultiPointConstraint(V) - mpc.create_periodic_constraint_topological(V.sub(0), mt, 2, periodic_relation, bcs, 1) + mpc.create_periodic_constraint_topological(V.sub(0), mt, 2, periodic_relation, bcs, PETSc.ScalarType(1)) mpc.finalize() # Define variational problem u = TrialFunction(V) diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index aada0771..8ffb3494 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -31,57 +31,54 @@ namespace py = pybind11; -namespace dolfinx_mpc_wrappers -{ -void mpc(py::module& m) +namespace + { - m.def("compute_shared_indices", &dolfinx_mpc::compute_shared_indices); +// Templating over mesh resolution +template +void declare_mpc(py::module& m, std::string type) +{ + std::string pyclass_name = "MultiPointConstraint_" + type; // dolfinx_mpc::MultiPointConstraint - py::class_< - dolfinx_mpc::MultiPointConstraint, - std::shared_ptr>> - multipointconstraint( - m, "MultiPointConstraint", - "Object for representing contact (non-penetrating) conditions"); + py::class_, + std::shared_ptr>> + multipointconstraint(m, pyclass_name.c_str(), + "Object for representing contact " + "(non-penetrating) conditions"); multipointconstraint - .def(py::init>, + .def(py::init>, std::vector, std::vector, - std::vector, std::vector, + std::vector, std::vector, std::vector>()) - .def_property_readonly( - "masters", - &dolfinx_mpc::MultiPointConstraint::masters) + .def_property_readonly("masters", + &dolfinx_mpc::MultiPointConstraint::masters) .def("coefficients", - [](dolfinx_mpc::MultiPointConstraint& self) + [](dolfinx_mpc::MultiPointConstraint& self) { - std::shared_ptr> - adj = self.coefficients(); + std::shared_ptr> adj + = self.coefficients(); const std::vector& offsets = adj->offsets(); - const std::vector& data = adj->array(); + const std::vector& data = adj->array(); - return std::pair, - py::array_t>( - py::array_t(data.size(), data.data(), - py::cast(self)), + 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))); }) .def_property_readonly( "constants", - [](dolfinx_mpc::MultiPointConstraint& self) + [](dolfinx_mpc::MultiPointConstraint& self) { - const std::vector& consts = self.constant_values(); - return py::array_t(consts.size(), consts.data(), - py::cast(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("owners", + &dolfinx_mpc::MultiPointConstraint::owners) .def_property_readonly( "slaves", - [](dolfinx_mpc::MultiPointConstraint& self) + [](dolfinx_mpc::MultiPointConstraint& self) { const std::vector& slaves = self.slaves(); return py::array_t(slaves.size(), slaves.data(), @@ -89,7 +86,7 @@ void mpc(py::module& m) }) .def_property_readonly( "is_slave", - [](dolfinx_mpc::MultiPointConstraint& self) + [](dolfinx_mpc::MultiPointConstraint& self) { std::span slaves = self.is_slave(); return py::array_t(slaves.size(), slaves.data(), @@ -98,42 +95,118 @@ void mpc(py::module& m) .def_property_readonly( "cell_to_slaves", - &dolfinx_mpc::MultiPointConstraint::cell_to_slaves) + &dolfinx_mpc::MultiPointConstraint::cell_to_slaves) .def_property_readonly( "num_local_slaves", - &dolfinx_mpc::MultiPointConstraint::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) + &dolfinx_mpc::MultiPointConstraint::function_space) + .def_property_readonly("owners", + &dolfinx_mpc::MultiPointConstraint::owners) .def( "backsubstitution", - [](dolfinx_mpc::MultiPointConstraint& self, - py::array_t u) { - self.backsubstitution( - std::span(u.mutable_data(), u.size())); - }, + [](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") .def( "homogenize", - [](dolfinx_mpc::MultiPointConstraint& self, - py::array_t u) { - self.homogenize(std::span(u.mutable_data(), u.size())); - }, + [](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"); - py::class_, - std::shared_ptr>> - mpc_data(m, "mpc_data", "Object with data arrays for mpc"); + // .def("ghost_masters", &dolfinx_mpc::mpc_data::ghost_masters); +} + +template +void declare_functions(py::module& m) +{ + m.def("compute_shared_indices", &dolfinx_mpc::compute_shared_indices); + + m.def("create_sparsity_pattern", &dolfinx_mpc::create_sparsity_pattern); + + m.def("create_contact_slip_condition", + &dolfinx_mpc::create_contact_slip_condition); + 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, + const std::vector>>& + bcs, + T scale, bool collapse) + { + auto _indicator + = [&indicator](MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const U, + MDSPAN_IMPL_STANDARD_NAMESPACE::extents< + std::size_t, 3, + MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>> + 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); + std::vector s(m.data(), m.data() + m.size()); + return s; + }; + + 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); + std::vector output(v.data(), v.data() + v.size()); + return output; + }; + return dolfinx_mpc::create_periodic_condition_geometrical( + V, _indicator, _relation, bcs, scale, collapse); + }); + + m.def( + "create_periodic_constraint_topological", + [](const std::shared_ptr>& V, + const std::shared_ptr>& + meshtags, + const int dim, + const std::function(const py::array_t&)>& 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); + 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); + }); +} + +template +void declare_mpc_data(py::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( "slaves", - [](dolfinx_mpc::mpc_data& self) + [](dolfinx_mpc::mpc_data& self) { const std::vector& slaves = self.slaves; return py::array_t(slaves.size(), slaves.data(), @@ -141,7 +214,7 @@ void mpc(py::module& m) }) .def_property_readonly( "masters", - [](dolfinx_mpc::mpc_data& self) + [](dolfinx_mpc::mpc_data& self) { const std::vector& masters = self.masters; return py::array_t(masters.size(), masters.data(), @@ -149,15 +222,14 @@ void mpc(py::module& m) }) .def_property_readonly( "coeffs", - [](dolfinx_mpc::mpc_data& self) + [](dolfinx_mpc::mpc_data& self) { - const std::vector& coeffs = self.coeffs; - return py::array_t(coeffs.size(), coeffs.data(), - py::cast(self)); + const std::vector& coeffs = self.coeffs; + return py::array_t(coeffs.size(), coeffs.data(), py::cast(self)); }) .def_property_readonly( "owners", - [](dolfinx_mpc::mpc_data& self) + [](dolfinx_mpc::mpc_data& self) { const std::vector& owners = self.owners; return py::array_t(owners.size(), owners.data(), @@ -165,38 +237,46 @@ void mpc(py::module& m) }) .def_property_readonly( "offsets", - [](dolfinx_mpc::mpc_data& self) + [](dolfinx_mpc::mpc_data& self) { const std::vector& offsets = self.offsets; return py::array_t(offsets.size(), offsets.data(), py::cast(self)); }); +} - // .def("ghost_masters", &dolfinx_mpc::mpc_data::ghost_masters); - m.def("create_sparsity_pattern", - &dolfinx_mpc::create_sparsity_pattern); - - m.def("assemble_matrix", - [](Mat A, const dolfinx::fem::Form& a, - const std::shared_ptr>& mpc0, - const std::shared_ptr>& mpc1, - const std::vector>>& bcs, - const PetscScalar diagval) +template +void declare_petsc_functions(py::module& m) +{ + m.def("create_normal_approximation", + [](std::shared_ptr> V, std::int32_t dim, + const py::array_t& entities) { - dolfinx_mpc::assemble_matrix( - dolfinx::la::petsc::Matrix::set_block_fn(A, ADD_VALUES), - dolfinx::la::petsc::Matrix::set_fn(A, ADD_VALUES), a, mpc0, mpc1, - bcs, diagval); + return dolfinx_mpc::create_normal_approximation( + V, dim, + std::span(entities.data(), entities.size())); }); + m.def( + "assemble_matrix", + [](Mat A, const dolfinx::fem::Form& a, + const std::shared_ptr>& + mpc0, + const std::shared_ptr>& + mpc1, + const std::vector>>& + bcs, + const T diagval) + { + dolfinx_mpc::assemble_matrix( + dolfinx::la::petsc::Matrix::set_block_fn(A, ADD_VALUES), + dolfinx::la::petsc::Matrix::set_fn(A, ADD_VALUES), a, mpc0, mpc1, + bcs, diagval); + }); m.def( "assemble_vector", - [](py::array_t b, - const dolfinx::fem::Form& L, - const std::shared_ptr< - const dolfinx_mpc::MultiPointConstraint>& mpc) + [](py::array_t b, const dolfinx::fem::Form& L, + const std::shared_ptr>& + mpc) { dolfinx_mpc::assemble_vector(std::span(b.mutable_data(), b.size()), L, mpc); @@ -206,16 +286,14 @@ void mpc(py::module& m) m.def( "apply_lifting", - [](py::array_t b, - std::vector>>& a, - const std::vector>>>& bcs1, - const std::vector>& x0, - double scale, - std::shared_ptr< - const dolfinx_mpc::MultiPointConstraint>& mpc) + [](py::array_t b, + std::vector>>& a, + const std::vector>>>& bcs1, + const std::vector>& x0, U scale, + std::shared_ptr>& mpc) { - std::vector> _x0; + std::vector> _x0; for (const auto& x : x0) _x0.emplace_back(x.data(), x.size()); @@ -228,9 +306,8 @@ void mpc(py::module& m) m.def( "create_matrix", - [](const dolfinx::fem::Form& a, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint>& mpc) + [](const dolfinx::fem::Form& a, + const std::shared_ptr>& mpc) { auto A = dolfinx_mpc::create_matrix(a, mpc); Mat _A = A.mat(); @@ -241,11 +318,9 @@ void mpc(py::module& m) "Create a PETSc Mat for bilinear form."); m.def( "create_matrix", - [](const dolfinx::fem::Form& a, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint>& mpc0, - const std::shared_ptr< - dolfinx_mpc::MultiPointConstraint>& mpc1) + [](const dolfinx::fem::Form& a, + const std::shared_ptr>& mpc0, + const std::shared_ptr>& mpc1) { auto A = dolfinx_mpc::create_matrix(a, mpc0, mpc1); Mat _A = A.mat(); @@ -254,83 +329,31 @@ void mpc(py::module& m) }, py::return_value_policy::take_ownership, "Create a PETSc Mat for bilinear form."); - m.def("create_contact_slip_condition", - &dolfinx_mpc::create_contact_slip_condition); - 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_normal_approximation", - [](std::shared_ptr> V, - std::int32_t dim, - const py::array_t& entities) - { - return dolfinx_mpc::create_normal_approximation( - V, dim, - std::span(entities.data(), entities.size())); - }); +} - 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::vector>>& bcs, - PetscScalar scale, bool collapse) - { - auto _indicator - = [&indicator]( - std::experimental::mdspan< - const double, - std::experimental::extents< - std::size_t, 3, std::experimental::dynamic_extent>> - 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); - std::vector s(m.data(), m.data() + m.size()); - return s; - }; +} // namespace - 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); - std::vector output(v.data(), v.data() + v.size()); - return output; - }; - return dolfinx_mpc::create_periodic_condition_geometrical( - V, _indicator, _relation, bcs, scale, collapse); - }); +namespace dolfinx_mpc_wrappers +{ - m.def("create_periodic_constraint_topological", - [](const std::shared_ptr>& V, - const std::shared_ptr>& - meshtags, - const int dim, - const std::function(const py::array_t&)>& - relation, - const std::vector>>& bcs, - PetscScalar 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); - 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); - }); +void mpc(py::module& m) +{ + + declare_mpc(m, "float"); + declare_mpc, float>(m, "complex_float"); + declare_mpc(m, "double"); + declare_mpc, double>(m, "complex_double"); + + declare_functions(m); + declare_functions, float>(m); + declare_functions(m); + declare_functions, double>(m); + + declare_mpc_data(m, "float"); + declare_mpc_data, float>(m, "complex_float"); + declare_mpc_data(m, "double"); + declare_mpc_data, double>(m, "complex_double"); + + declare_petsc_functions(m); } } // namespace dolfinx_mpc_wrappers diff --git a/python/dolfinx_mpc/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index e132d128..58145f4c 100644 --- a/python/dolfinx_mpc/multipointconstraint.py +++ b/python/dolfinx_mpc/multipointconstraint.py @@ -1,4 +1,4 @@ -# Copyright (C) 2020-2021 Jørgen S. Dokken +# Copyright (C) 2020-2023 Jørgen S. Dokken # # This file is part of DOLFINX_MPC # @@ -11,12 +11,41 @@ import dolfinx.mesh as _mesh import numpy import numpy.typing as npt +from dolfinx import default_scalar_type from petsc4py import PETSc as _PETSc import dolfinx_mpc.cpp from .dictcondition import create_dictionary_constraint +_mpc_classes = Union[dolfinx_mpc.cpp.mpc.MultiPointConstraint_double, dolfinx_mpc.cpp.mpc.MultiPointConstraint_float, + dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double, + dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float] +_float_classes = Union[numpy.float32, numpy.float64, numpy.complex128, numpy.complex64] +_float_array_types = Union[npt.NDArray[numpy.float32], npt.NDArray[numpy.float64], npt.NDArray[numpy.complex64], + npt.NDArray[numpy.complex128]] +_mpc_data_classes = Union[dolfinx_mpc.cpp.mpc.mpc_data_double, dolfinx_mpc.cpp.mpc.mpc_data_float, + dolfinx_mpc.cpp.mpc.mpc_data_complex_double, + dolfinx_mpc.cpp.mpc.mpc_data_complex_float] + + +class MPCData(): + _cpp_object: _mpc_data_classes + + def __init__(self, slaves: npt.NDArray[numpy.int32], + masters: npt.NDArray[numpy.int64], coeffs: _float_array_types, + owners: npt.NDArray[numpy.int32], offsets: npt.NDArray[numpy.int32]): + if coeffs.dtype.type == numpy.float32: + self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_float(slaves, masters, coeffs, owners, offsets) + elif coeffs.dtype.type == numpy.float64: + self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_double(slaves, masters, coeffs, owners, offsets) + elif coeffs.dtype.type == numpy.complex64: + self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_float(slaves, masters, coeffs, owners, offsets) + elif coeffs.dtype.type == numpy.complex128: + self._cpp_object = dolfinx_mpc.cpp.mpc.mpc_data_complex_double(slaves, masters, coeffs, owners, offsets) + else: + raise ValueError("Unsupported dtype {coeffs.dtype.type} for coefficients") + class MultiPointConstraint(): """ @@ -25,28 +54,31 @@ class MultiPointConstraint(): Args: V: The function space + dtype: The dtype of the underlying functions """ _slaves: npt.NDArray[numpy.int32] _masters: npt.NDArray[numpy.int64] - _coeffs: npt.NDArray[_PETSc.ScalarType] + _coeffs: _float_array_types _owners: npt.NDArray[numpy.int32] _offsets: npt.NDArray[numpy.int32] - V: _fem.FunctionSpace + V: _fem.FunctionSpaceBase finalized: bool - _cpp_object: dolfinx_mpc.cpp.mpc.MultiPointConstraint + _cpp_object: _mpc_classes + _dtype: _float_classes __slots__ = tuple(__annotations__) - def __init__(self, V: _fem.FunctionSpace): + def __init__(self, V: _fem.FunctionSpaceBase, dtype: _float_classes = default_scalar_type): self._slaves = numpy.array([], dtype=numpy.int32) self._masters = numpy.array([], dtype=numpy.int64) - self._coeffs = numpy.array([], dtype=_PETSc.ScalarType) + self._coeffs = numpy.array([], dtype=dtype) self._owners = numpy.array([], dtype=numpy.int32) self._offsets = numpy.array([0], dtype=numpy.int32) self.V = V self.finalized = False + self._dtype = dtype - def add_constraint(self, V: _fem.FunctionSpace, slaves: npt.NDArray[numpy.int32], - masters: npt.NDArray[numpy.int64], coeffs: npt.NDArray[_PETSc.ScalarType], + def add_constraint(self, V: _fem.FunctionSpaceBase, slaves: npt.NDArray[numpy.int32], + masters: npt.NDArray[numpy.int64], coeffs: _float_array_types, owners: npt.NDArray[numpy.int32], offsets: npt.NDArray[numpy.int32]): """ Add new constraint given by numpy arrays. @@ -76,12 +108,18 @@ def add_constraint(self, V: _fem.FunctionSpace, slaves: npt.NDArray[numpy.int32] self._coeffs = numpy.append(self._coeffs, coeffs) self._owners = numpy.append(self._owners, owners) - def add_constraint_from_mpc_data(self, V: _fem.FunctionSpace, mpc_data: dolfinx_mpc.cpp.mpc.mpc_data): + def add_constraint_from_mpc_data(self, V: _fem.FunctionSpaceBase, + mpc_data: Union[_mpc_data_classes, MPCData]): """ Add new constraint given by an `dolfinc_mpc.cpp.mpc.mpc_data`-object """ self._already_finalized() - self.add_constraint(V, mpc_data.slaves, mpc_data.masters, mpc_data.coeffs, mpc_data.owners, mpc_data.offsets) + try: + self.add_constraint(V, mpc_data.slaves, mpc_data.masters, mpc_data.coeffs, # type: ignore + mpc_data.owners, mpc_data.offsets) # type: ignore + except AttributeError: + self.add_constraint(V, mpc_data._cpp_object.slaves, mpc_data._cpp_object.masters, + mpc_data._cpp_object.coeffs, mpc_data._cpp_object.owners, mpc_data._cpp_object.offsets) def finalize(self) -> None: """ @@ -90,20 +128,39 @@ def finalize(self) -> None: freedom and builds a new index map and function space where unghosted master dofs are added as ghosts. """ self._already_finalized() - + self._coeffs.astype(numpy.dtype(self._dtype)) # Initialize C++ object and create slave->cell maps - self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint( - self.V._cpp_object, self._slaves, self._masters, self._coeffs, self._owners, self._offsets) + if self._dtype == numpy.float32: + self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_float( + self.V._cpp_object, self._slaves, self._masters, self._coeffs.astype(self._dtype), + self._owners, self._offsets) + elif self._dtype == numpy.float64: + self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_double( + self.V._cpp_object, self._slaves, self._masters, self._coeffs.astype(self._dtype), + self._owners, self._offsets) + elif self._dtype == numpy.complex64: + self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_float( + self.V._cpp_object, self._slaves, self._masters, self._coeffs.astype(self._dtype), + self._owners, self._offsets) + + elif self._dtype == numpy.complex128: + self._cpp_object = dolfinx_mpc.cpp.mpc.MultiPointConstraint_complex_double( + self.V._cpp_object, self._slaves, self._masters, self._coeffs.astype(self._dtype), + self._owners, self._offsets) + else: + raise ValueError("Unsupported dtype {coeffs.dtype.type} for coefficients") + # Replace function space - self.V = _fem.FunctionSpace(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space) + self.V = _fem.FunctionSpaceBase(self.V.mesh, self.V.ufl_element(), self._cpp_object.function_space) self.finalized = True # Delete variables that are no longer required del (self._slaves, self._masters, self._coeffs, self._owners, self._offsets) - def create_periodic_constraint_topological(self, V: _fem.FunctionSpace, meshtag: _mesh.MeshTags, tag: int, + def create_periodic_constraint_topological(self, V: _fem.FunctionSpaceBase, meshtag: _mesh.MeshTags, tag: int, relation: Callable[[numpy.ndarray], numpy.ndarray], - bcs: List[_fem.DirichletBC], scale: _PETSc.ScalarType = 1): + bcs: List[_fem.DirichletBC], + scale: _float_classes = default_scalar_type(1.)): """ Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`. :math:`u(x_i) = scale * u(relation(x_i))` for all of :math:`x_i` on marked entities. @@ -127,10 +184,11 @@ def create_periodic_constraint_topological(self, V: _fem.FunctionSpace, meshtag: raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC") self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data) - def create_periodic_constraint_geometrical(self, V: _fem.FunctionSpace, + def create_periodic_constraint_geometrical(self, V: _fem.FunctionSpaceBase, indicator: Callable[[numpy.ndarray], numpy.ndarray], relation: Callable[[numpy.ndarray], numpy.ndarray], - bcs: List[_fem.DirichletBC], scale: _PETSc.ScalarType = 1): + bcs: List[_fem.DirichletBC], + scale: _float_classes = default_scalar_type(1.)): """ Create a periodic condition for all degrees of freedom whose physical location satisfies :math:`indicator(x_i)==True`, i.e. @@ -155,7 +213,7 @@ def create_periodic_constraint_geometrical(self, V: _fem.FunctionSpace, raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC") self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data) - def create_slip_constraint(self, space: _fem.FunctionSpace, facet_marker: Tuple[_mesh.MeshTags, int], + def create_slip_constraint(self, space: _fem.FunctionSpaceBase, facet_marker: Tuple[_mesh.MeshTags, int], v: _fem.Function, bcs: List[_fem.DirichletBC] = []): """ Create a slip constraint :math:`u \\cdot v=0` over the entities defined in `facet_marker` with the given index. @@ -186,7 +244,7 @@ def create_slip_constraint(self, space: _fem.FunctionSpace, facet_marker: Tuple[ Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,)) Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1) me = basix.ufl.mixed_element([Ve, Qe], gdim=mesh.geometry.dim) - W = dolfinx.fem.FunctionSpace(mesh, me) + W = dolfinx.fem.functionspace(mesh, me) mpc = MultiPointConstraint(W) n_space, _ = W.sub(0).collapse() normal = dolfinx.fem.Function(n_space) @@ -202,7 +260,7 @@ def create_slip_constraint(self, space: _fem.FunctionSpace, facet_marker: Tuple[ Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,)) Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1) me = basix.ufl.mixed_element([Ve, Qe], gdim=mesh.geometry.dim) - W = dolfinx.fem.FunctionSpace(mesh, me) + W = dolfinx.fem.functionspace(mesh, me) mpc = MultiPointConstraint(W) n_space, _ = W.sub(0).collapse() normal = Function(n_space) @@ -303,9 +361,9 @@ def slaves(self): return self._cpp_object.slaves @property - def masters(self): + def masters(self) -> _cpp.graph.AdjacencyList_int32: """ - Returns a `dolfinx.cpp.graph.AdjacencyList_int32` whose ith node corresponds to + Returns an adjacency-list whose ith node corresponds to a degree of freedom (local to process), and links the corresponding master dofs (local to process). Examples: @@ -319,7 +377,7 @@ def masters(self): self._not_finalized() return self._cpp_object.masters - def coefficients(self): + def coefficients(self) -> _float_array_types: """ Returns a vector containing the coefficients for the constraint, and the corresponding offsets for the ith degree of freedom. diff --git a/python/dolfinx_mpc/utils/mpc_utils.py b/python/dolfinx_mpc/utils/mpc_utils.py index 408ad20b..3d727543 100644 --- a/python/dolfinx_mpc/utils/mpc_utils.py +++ b/python/dolfinx_mpc/utils/mpc_utils.py @@ -147,7 +147,7 @@ def log_info(message): _log.set_log_level(old_level) -def rigid_motions_nullspace(V: _fem.FunctionSpace): +def rigid_motions_nullspace(V: _fem.FunctionSpaceBase): """ Function to build nullspace for 2D/3D elasticity. @@ -397,7 +397,7 @@ def create_point_to_point_constraint(V, slave_point, master_point, vector=None): return slaves, masters, coeffs, owners, offsets -def create_normal_approximation(V: _fem.FunctionSpace, mt: _cpp.mesh.MeshTags_int32, value: int): +def create_normal_approximation(V: _fem.FunctionSpaceBase, mt: _cpp.mesh.MeshTags_int32, value: int): """ Creates a normal approximation for the dofs in the closure of the attached entities. Where a dof is attached to entities facets, an average is computed diff --git a/python/tests/test_cube_contact.py b/python/tests/test_cube_contact.py index 683580ca..e5e7cec2 100644 --- a/python/tests/test_cube_contact.py +++ b/python/tests/test_cube_contact.py @@ -174,7 +174,7 @@ def test_cube_contact(generate_hex_boxes, nonslip, get_assemblers): # noqa: F81 mesh, mt = mesh_data fdim = mesh.topology.dim - 1 # Create functionspaces - V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) # Helper for orienting traction diff --git a/python/tests/test_rectangular_assembly.py b/python/tests/test_rectangular_assembly.py index 392d79fa..57a8fc70 100644 --- a/python/tests/test_rectangular_assembly.py +++ b/python/tests/test_rectangular_assembly.py @@ -50,9 +50,9 @@ def test_mixed_element(cell_type, ghost_mode): Ve = basix.ufl.element(basix.ElementFamily.P, cellname, 2, shape=(mesh.geometry.dim,)) Qe = basix.ufl.element(basix.ElementFamily.P, cellname, 1) - V = dolfinx.fem.FunctionSpace(mesh, Ve) - Q = dolfinx.fem.FunctionSpace(mesh, Qe) - W = dolfinx.fem.FunctionSpace(mesh, Ve * Qe) + V = dolfinx.fem.functionspace(mesh, Ve) + Q = dolfinx.fem.functionspace(mesh, Qe) + W = dolfinx.fem.functionspace(mesh, Ve * Qe) inlet_velocity = dolfinx.fem.Function(V) inlet_velocity.interpolate( diff --git a/python/tests/test_surface_integral.py b/python/tests/test_surface_integral.py index 096599fc..f338e602 100644 --- a/python/tests/test_surface_integral.py +++ b/python/tests/test_surface_integral.py @@ -26,7 +26,7 @@ def test_surface_integrals(get_assemblers): # noqa: F811 N = 4 mesh = create_unit_square(MPI.COMM_WORLD, N, N) - V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, ))) # Fixed Dirichlet BC on the left wall def left_wall(x): @@ -144,7 +144,7 @@ def test_surface_integral_dependency(get_assemblers): # noqa: F811 assemble_matrix, assemble_vector = get_assemblers N = 10 mesh = create_unit_square(MPI.COMM_WORLD, N, N) - V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,))) def top(x): return np.isclose(x[1], 1) diff --git a/python/tests/test_vector_poisson.py b/python/tests/test_vector_poisson.py index 859a8afc..0d9827c7 100644 --- a/python/tests/test_vector_poisson.py +++ b/python/tests/test_vector_poisson.py @@ -29,7 +29,7 @@ def test_vector_possion(Nx, Ny, slave_space, master_space, get_assemblers): # n # Create mesh and function space mesh = create_unit_square(MPI.COMM_WORLD, Nx, Ny) - V = fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) + V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,))) def boundary(x): return np.isclose(x.T, [0, 0, 0]).all(axis=1)