diff --git a/cpp/dolfinx/fem/DofMap.cpp b/cpp/dolfinx/fem/DofMap.cpp index b39c0f6dcf5..01358c20a0a 100644 --- a/cpp/dolfinx/fem/DofMap.cpp +++ b/cpp/dolfinx/fem/DofMap.cpp @@ -225,11 +225,11 @@ std::pair> DofMap::collapse( // Create new element dof layout and reset parent ElementDofLayout collapsed_dof_layout = layout.copy(); - auto [_index_map, bs, dofmap] - = build_dofmap_data(comm, topology, collapsed_dof_layout, reorder_fn); + auto [_index_map, bs, dofmaps] = build_dofmap_data( + comm, topology, {collapsed_dof_layout}, reorder_fn); auto index_map = std::make_shared(std::move(_index_map)); - return DofMap(layout, index_map, bs, std::move(dofmap), bs); + return DofMap(layout, index_map, bs, std::move(dofmaps.front()), bs); } else { diff --git a/cpp/dolfinx/fem/dofmapbuilder.cpp b/cpp/dolfinx/fem/dofmapbuilder.cpp index 2dc73c2ada6..60ecf1cbf0a 100644 --- a/cpp/dolfinx/fem/dofmapbuilder.cpp +++ b/cpp/dolfinx/fem/dofmapbuilder.cpp @@ -541,10 +541,10 @@ std::pair, std::vector> get_global_indices( } // namespace //----------------------------------------------------------------------------- -std::tuple> +std::tuple>> fem::build_dofmap_data( MPI_Comm comm, const mesh::Topology& topology, - const ElementDofLayout& element_dof_layout, + const std::vector& element_dof_layouts, const std::function( const graph::AdjacencyList&)>& reorder_fn) { @@ -557,12 +557,12 @@ fem::build_dofmap_data( // pair {dimension, mesh entity index} giving the mesh entity that dof // i is associated with. const auto [node_graph0, local_to_global0, dof_entity0] - = build_basic_dofmap(topology, element_dof_layout); + = build_basic_dofmap(topology, element_dof_layouts.front()); std::vector> index_maps(D + 1); for (int d = 0; d <= D; ++d) { - if (element_dof_layout.num_entity_dofs(d) > 0) + if (element_dof_layouts.front().num_entity_dofs(d) > 0) { assert(topology.index_map(d)); index_maps[d] = topology.index_map(d); @@ -578,10 +578,10 @@ fem::build_dofmap_data( std::int64_t offset = 0; for (int d = 0; d <= D; ++d) { - if (element_dof_layout.num_entity_dofs(d) > 0) + if (element_dof_layouts.front().num_entity_dofs(d) > 0) { offset += index_maps[d]->local_range()[0] - * element_dof_layout.num_entity_dofs(d); + * element_dof_layouts.front().num_entity_dofs(d); } } @@ -596,11 +596,12 @@ fem::build_dofmap_data( local_to_global_owner); // Build re-ordered dofmap - std::vector dofmap(node_graph0.array.size()); - for (std::size_t i = 0; i < dofmap.size(); ++i) - dofmap[i] = old_to_new[node_graph0.array[i]]; + std::vector> dofmaps(1); + dofmaps.front().resize(node_graph0.array.size()); + for (std::size_t i = 0; i < dofmaps.front().size(); ++i) + dofmaps.front()[i] = old_to_new[node_graph0.array[i]]; - return {std::move(index_map), element_dof_layout.block_size(), - std::move(dofmap)}; + return {std::move(index_map), element_dof_layouts.front().block_size(), + std::move(dofmaps)}; } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/fem/dofmapbuilder.h b/cpp/dolfinx/fem/dofmapbuilder.h index ea78799e5bf..28460d8267b 100644 --- a/cpp/dolfinx/fem/dofmapbuilder.h +++ b/cpp/dolfinx/fem/dofmapbuilder.h @@ -35,13 +35,14 @@ class ElementDofLayout; /// Build dofmap data for elements on a mesh topology /// @param[in] comm MPI communicator /// @param[in] topology The mesh topology -/// @param[in] element_dof_layout The element dof layout +/// @param[in] element_dof_layouts The element dof layouts for each cell type in +/// @p topology /// @param[in] reorder_fn Graph reordering function that is applied to /// the dofmaps -/// @return The index map, block size, and dofmap for element type 0 -std::tuple> +/// @return The index map, block size, and dofmaps for each element type +std::tuple>> build_dofmap_data(MPI_Comm comm, const mesh::Topology& topology, - const ElementDofLayout& element_dof_layout, + const std::vector& element_dof_layouts, const std::function( const graph::AdjacencyList&)>& reorder_fn); diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index ed7fc2dab7e..ecd2f69ae7a 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -110,8 +110,8 @@ fem::DofMap fem::create_dofmap( topology.create_entities(d); } - auto [_index_map, bs, dofmap] - = build_dofmap_data(comm, topology, layout, reorder_fn); + auto [_index_map, bs, dofmaps] + = build_dofmap_data(comm, topology, {layout}, reorder_fn); auto index_map = std::make_shared(std::move(_index_map)); // If the element's DOF transformations are permutations, permute the @@ -126,12 +126,12 @@ fem::DofMap fem::create_dofmap( int dim = layout.num_dofs(); for (std::int32_t cell = 0; cell < num_cells; ++cell) { - std::span dofs(dofmap.data() + cell * dim, dim); + std::span dofs(dofmaps.front().data() + cell * dim, dim); unpermute_dofs(dofs, cell_info[cell]); } } - return DofMap(layout, index_map, bs, std::move(dofmap), bs); + return DofMap(layout, index_map, bs, std::move(dofmaps.front()), bs); } //----------------------------------------------------------------------------- std::vector fem::get_coefficient_names(const ufcx_form& ufcx_form) diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h index 217a3bae7f7..ceaebbae69c 100644 --- a/cpp/dolfinx/mesh/Geometry.h +++ b/cpp/dolfinx/mesh/Geometry.h @@ -58,8 +58,38 @@ class Geometry const fem::CoordinateElement< typename std::remove_reference_t>& element, V&& x, int dim, W&& input_global_indices) - : _dim(dim), _dofmap(std::forward(dofmap)), _index_map(index_map), - _cmap(element), _x(std::forward(x)), + : _dim(dim), _dofmaps({dofmap}), _index_map(index_map), _cmaps({element}), + _x(std::forward(x)), + _input_global_indices(std::forward(input_global_indices)) + { + assert(_x.size() % 3 == 0); + if (_x.size() / 3 != _input_global_indices.size()) + throw std::runtime_error("Geometry size mis-match"); + } + + /// @brief Constructor of object that holds mesh geometry data. + /// + /// @param[in] index_map Index map associated with the geometry dofmap + /// @param[in] dofmaps The geometry (point) dofmaps. For a cell, it + /// gives the position in the point array of each local geometry node + /// @param[in] elements Elements that describes the cell geometry maps. + /// @param[in] x The point coordinates. The shape is `(num_points, 3)` + /// and the storage is row-major. + /// @param[in] dim The geometric dimension (`0 < dim <= 3`). + /// @param[in] input_global_indices The 'global' input index of each + /// point, commonly from a mesh input file. + template + requires std::is_convertible_v, std::vector> + and std::is_convertible_v, + std::vector> + Geometry( + std::shared_ptr index_map, + const std::vector>& dofmaps, + const std::vector>>& elements, + V&& x, int dim, W&& input_global_indices) + : _dim(dim), _dofmaps(dofmaps), _index_map(index_map), _cmaps(elements), + _x(std::forward(x)), _input_global_indices(std::forward(input_global_indices)) { assert(_x.size() % 3 == 0); @@ -85,20 +115,47 @@ class Geometry /// Return Euclidean dimension of coordinate system int dim() const { return _dim; } - /// DOF map + /// @brief DofMap for the geometry + /// @return A 2D array with shape [num_cells, dofs_per_cell] MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> dofmap() const { - int ndofs = _cmap.dim(); + if (_dofmaps.size() != 1) + throw std::runtime_error("Multiple dofmaps"); + + int ndofs = _cmaps.front().dim(); return MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>( - _dofmap.data(), _dofmap.size() / ndofs, ndofs); + _dofmaps.front().data(), _dofmaps.front().size() / ndofs, ndofs); } - /// Index map + /// @brief The dofmap associated with the `i`th coordinate map in the + /// geometry. + /// @param i Index + /// @return A 2D array with shape [num_cells, dofs_per_cell] + MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> + dofmap(std::int32_t i) const + { + if (i < 0 or i >= (int)_dofmaps.size()) + { + throw std::out_of_range("Cannot get dofmap:" + std::to_string(i) + + " out of range"); + } + int ndofs = _cmaps[i].dim(); + + return MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< + const std::int32_t, + MDSPAN_IMPL_STANDARD_NAMESPACE::dextents>( + _dofmaps[i].data(), _dofmaps[i].size() / ndofs, ndofs); + } + + /// @brief Index map + /// @return The index map for the geometry dofs std::shared_ptr index_map() const { return _index_map; @@ -117,10 +174,28 @@ class Geometry /// (num_points, 3) std::span x() { return _x; } - /// @brief The elements that describes the geometry maps. + /// @brief The element that describes the geometry map. /// /// @return The coordinate/geometry element - const fem::CoordinateElement& cmap() const { return _cmap; } + const fem::CoordinateElement& cmap() const + { + if (_cmaps.size() != 1) + throw std::runtime_error("Multiple cmaps."); + return _cmaps.front(); + } + + /// @brief The element that describe the `i`th geometry map + /// @param i Index of the coordinate element + /// @return Coordinate element + const fem::CoordinateElement& cmap(std::int32_t i) const + { + if (i < 0 or i >= (int)_cmaps.size()) + { + throw std::out_of_range("Cannot get cmap:" + std::to_string(i) + + " out of range"); + } + return _cmaps[i]; + } /// Global user indices const std::vector& input_global_indices() const @@ -132,14 +207,14 @@ class Geometry // Geometric dimension int _dim; - // Map per cell for extracting coordinate data - std::vector _dofmap; + // Map per cell for extracting coordinate data for each cmap + std::vector> _dofmaps; // IndexMap for geometry 'dofmap' std::shared_ptr _index_map; // The coordinate elements - fem::CoordinateElement _cmap; + std::vector> _cmaps; // Coordinates for all points stored as a contiguous array (row-major, // column size = 3) @@ -153,12 +228,131 @@ class Geometry /// Template type deduction template Geometry(std::shared_ptr, U, - const fem::CoordinateElement< - typename std::remove_reference_t>&, + const std::vector>>&, V, int, W) -> Geometry>; /// @endcond +/// @brief Build Geometry from input data. +/// +/// This function should be called after the mesh topology is built and +/// 'node' coordinate data has been distributed to the processes where +/// it is required. +/// +/// @param[in] topology Mesh topology. +/// @param[in] elements List of elements that defines the geometry map for +/// each cell type. +/// @param[in] nodes Geometry node global indices for cells on this +/// process. @pre Must be sorted. +/// @param[in] xdofs Geometry degree-of-freedom map (using global +/// indices) for cells on this process. `nodes` is a sorted and unique +/// list of the indices in `xdofs`. +/// @param[in] x The node coordinates (row-major, with shape +/// `(num_nodes, dim)`. The global index of each node is `i + +/// rank_offset`, where `i` is the local row index in `x` and +/// `rank_offset` is the sum of `x` rows on all processed with a lower +/// rank than the caller. +/// @param[in] dim Geometric dimension (1, 2, or 3). +/// @param[in] reorder_fn Function for re-ordering the degree-of-freedom +/// map associated with the geometry data. +/// @note Experimental new interface for multiple cmap/dofmap +/// @return A mesh geometry. +template +Geometry> +create_geometry( + const Topology& topology, + const std::vector>>& elements, + std::span nodes, std::span xdofs, + const U& x, int dim, + std::function(const graph::AdjacencyList&)> + reorder_fn + = nullptr) +{ + LOG(INFO) << "Create Geometry (multiple)"; + + assert(std::is_sorted(nodes.begin(), nodes.end())); + using T = typename std::remove_reference_t; + + // Check elements match cell types in topology + const int tdim = topology.dim(); + const std::size_t num_cell_types = topology.entity_types(tdim).size(); + if (elements.size() != num_cell_types) + throw std::runtime_error("Mismatch between topology and geometry."); + + std::vector dof_layouts; + for (const auto& el : elements) + dof_layouts.push_back(el.create_dof_layout()); + + LOG(INFO) << "Got " << dof_layouts.size() << " dof layouts"; + + // Build 'geometry' dofmap on the topology + auto [_dof_index_map, bs, dofmaps] + = fem::build_dofmap_data(topology.index_map(topology.dim())->comm(), + topology, dof_layouts, reorder_fn); + auto dof_index_map + = std::make_shared(std::move(_dof_index_map)); + + // If the mesh has higher order geometry, permute the dofmap + if (elements.first().needs_dof_permutations()) + { + const std::int32_t num_cells + = topology.connectivity(topology.dim(), 0)->num_nodes(); + const std::vector& cell_info + = topology.get_cell_permutation_info(); + int d = elements.front().dim(); + for (std::int32_t cell = 0; cell < num_cells; ++cell) + { + std::span dofs(dofmaps.front().data() + cell * d, d); + elements.front().unpermute_dofs(dofs, cell_info[cell]); + } + } + + LOG(INFO) << "Calling compute_local_to_global"; + // Compute local-to-global map from local indices in dofmap to the + // corresponding global indices in cells, and pass to function to + // compute local (dof) to local (position in coords) map from (i) + // local-to-global for dofs and (ii) local-to-global for entries in + // coords + + LOG(INFO) << "xdofs.size = " << xdofs.size(); + std::vector all_dofmaps; + std::stringstream s; + for (auto q : dofmaps) + { + s << q.size() << " "; + all_dofmaps.insert(all_dofmaps.end(), q.begin(), q.end()); + } + LOG(INFO) << "dofmap sizes = " << s.str(); + LOG(INFO) << "all_dofmaps.size = " << all_dofmaps.size(); + LOG(INFO) << "nodes.size = " << nodes.size(); + + const std::vector l2l = graph::build::compute_local_to_local( + graph::build::compute_local_to_global(xdofs, all_dofmaps), nodes); + + // Allocate space for input global indices and copy data + std::vector igi(nodes.size()); + std::transform(l2l.cbegin(), l2l.cend(), igi.begin(), + [&nodes](auto index) { return nodes[index]; }); + + // Build coordinate dof array, copying coordinates to correct position + assert(x.size() % dim == 0); + const std::size_t shape0 = x.size() / dim; + const std::size_t shape1 = dim; + std::vector xg(3 * shape0, 0); + for (std::size_t i = 0; i < shape0; ++i) + { + std::copy_n(std::next(x.cbegin(), shape1 * l2l[i]), shape1, + std::next(xg.begin(), 3 * i)); + } + + LOG(INFO) << "Creating geometry with " << dofmaps.size() << " dofmaps"; + + return Geometry(dof_index_map, std::move(dofmaps), elements, std::move(xg), + dim, std::move(igi)); +} + /// @brief Build Geometry from input data. /// /// This function should be called after the mesh topology is built and @@ -188,9 +382,8 @@ create_geometry( const Topology& topology, const fem::CoordinateElement< std::remove_reference_t>& element, - std::span nodes, - - std::span xdofs, const U& x, int dim, + std::span nodes, std::span xdofs, + const U& x, int dim, std::function(const graph::AdjacencyList&)> reorder_fn = nullptr) @@ -198,12 +391,12 @@ create_geometry( assert(std::is_sorted(nodes.begin(), nodes.end())); using T = typename std::remove_reference_t; - fem::ElementDofLayout doflayout = element.create_dof_layout(); + fem::ElementDofLayout dof_layout = element.create_dof_layout(); // Build 'geometry' dofmap on the topology - auto [_dof_index_map, bs, dofmap] + auto [_dof_index_map, bs, dofmaps] = fem::build_dofmap_data(topology.index_map(topology.dim())->comm(), - topology, doflayout, reorder_fn); + topology, {dof_layout}, reorder_fn); auto dof_index_map = std::make_shared(std::move(_dof_index_map)); @@ -217,7 +410,7 @@ create_geometry( int d = element.dim(); for (std::int32_t cell = 0; cell < num_cells; ++cell) { - std::span dofs(dofmap.data() + cell * d, d); + std::span dofs(dofmaps.front().data() + cell * d, d); element.unpermute_dofs(dofs, cell_info[cell]); } } @@ -228,7 +421,7 @@ create_geometry( // local-to-global for dofs and (ii) local-to-global for entries in // coords const std::vector l2l = graph::build::compute_local_to_local( - graph::build::compute_local_to_global(xdofs, dofmap), nodes); + graph::build::compute_local_to_global(xdofs, dofmaps.front()), nodes); // Allocate space for input global indices and copy data std::vector igi(nodes.size()); @@ -246,8 +439,8 @@ create_geometry( std::next(xg.begin(), 3 * i)); } - return Geometry(dof_index_map, std::move(dofmap), element, std::move(xg), dim, - std::move(igi)); + return Geometry(dof_index_map, std::move(dofmaps.front()), {element}, + std::move(xg), dim, std::move(igi)); } /// @brief Create a sub-geometry for a subset of entities. @@ -361,7 +554,7 @@ create_subgeometry(const Topology& topology, const Geometry& geometry, [&igi](std::int32_t sub_x_dof) { return igi[sub_x_dof]; }); // Create geometry - return {Geometry(sub_x_dof_index_map, std::move(sub_x_dofmap), sub_cmap, + return {Geometry(sub_x_dof_index_map, std::move(sub_x_dofmap), {sub_cmap}, std::move(sub_x), geometry.dim(), std::move(sub_igi)), std::move(subx_to_x_dofmap)}; } diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 7532875367f..d9466356786 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -985,8 +985,9 @@ void fem(nb::module_& m) [](MPICommWrapper comm, const dolfinx::mesh::Topology& topology, const dolfinx::fem::ElementDofLayout& layout) { + assert(topology.entity_types(topology.dim()).size() == 1); auto [map, bs, dofmap] = dolfinx::fem::build_dofmap_data( - comm.get(), topology, layout, + comm.get(), topology, {layout}, [](const dolfinx::graph::AdjacencyList& g) { return dolfinx::graph::reorder_gps(g); }); return std::tuple(std::move(map), bs, std::move(dofmap)); diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index e1a50d5d5cb..142b078a8f3 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -177,8 +177,9 @@ void declare_mesh(nb::module_& m, std::string type) nb::rv_policy::reference_internal, "Return coordinates of all geometry points. Each row is the " "coordinate of a point.") - .def_prop_ro("cmap", &dolfinx::mesh::Geometry::cmap, - "The coordinate map") + .def_prop_ro( + "cmap", [](dolfinx::mesh::Geometry& self) { return self.cmap(); }, + "The coordinate map") .def_prop_ro("input_global_indices", &dolfinx::mesh::Geometry::input_global_indices);