Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow promotion of geometry dimension for built-in mesh creators. #3471

Open
wants to merge 17 commits into
base: main
Choose a base branch
from
Open
2 changes: 1 addition & 1 deletion cpp/dolfinx/io/XDMFFile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ XDMFFile::read_mesh(const fem::CoordinateElement<double>& element,
// Create mesh
const std::vector<double>& _x = std::get<std::vector<double>>(x);
mesh::Mesh<double> mesh
= mesh::create_mesh(_comm.comm(), cells, element, _x, xshape, mode);
= mesh::create_mesh(_comm.comm(), cells, element, _x, xshape, xshape[1], mode);
mesh.name = name;
return mesh;
}
Expand Down
37 changes: 19 additions & 18 deletions cpp/dolfinx/mesh/Geometry.h
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,8 @@ class Geometry
_input_global_indices(std::forward<W>(input_global_indices))
{
assert(_x.size() % 3 == 0);
std::cout << _x.size() << std::endl;
std::cout << _input_global_indices.size() << std::endl;
if (_x.size() / 3 != _input_global_indices.size())
throw std::runtime_error("Geometry size mis-match");
}
Expand Down Expand Up @@ -94,6 +96,8 @@ class Geometry
_input_global_indices(std::forward<W>(input_global_indices))
{
assert(_x.size() % 3 == 0);
std::cout << _x.size() << std::endl;
std::cout << _input_global_indices.size() << std::endl;
if (_x.size() / 3 != _input_global_indices.size())
throw std::runtime_error("Geometry size mis-match");
}
Expand Down Expand Up @@ -254,7 +258,8 @@ Geometry(std::shared_ptr<const common::IndexMap>, U,
/// 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] xshape The shape of `x`.
/// @param[in] gdim Euclidean dimension of ambient space.
/// @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
Expand All @@ -266,7 +271,7 @@ create_geometry(
const std::vector<fem::CoordinateElement<
std::remove_reference_t<typename U::value_type>>>& elements,
std::span<const std::int64_t> nodes, std::span<const std::int64_t> xdofs,
const U& x, int dim,
const U& x, std::array<std::size_t, 2> xshape, int gdim,
std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
reorder_fn
= nullptr)
Expand Down Expand Up @@ -338,20 +343,17 @@ create_geometry(
[&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<T> xg(3 * shape0, 0);
for (std::size_t i = 0; i < shape0; ++i)
std::vector<T> xg(3 * xshape[0], 0);
for (std::size_t i = 0; i < xshape[0]; ++i)
{
std::copy_n(std::next(x.begin(), shape1 * l2l[i]), shape1,
std::copy_n(std::next(x.begin(), xshape[1] * l2l[i]), xshape[1],
std::next(xg.begin(), 3 * i));
}

spdlog::info("Creating geometry with {} dofmaps", dof_layouts.size());

return Geometry(dof_index_map, std::move(dofmaps), elements, std::move(xg),
dim, std::move(igi));
gdim, std::move(igi));
}

/// @brief Build Geometry from input data.
Expand All @@ -373,7 +375,8 @@ create_geometry(
/// 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] xshape Shape of the `x` data.
/// @param[in] gdim Geometric dimension of the ambient space.
/// @param[in] reorder_fn Function for re-ordering the degree-of-freedom
/// map associated with the geometry data.
/// @return A mesh geometry.
Expand All @@ -384,7 +387,7 @@ create_geometry(
const fem::CoordinateElement<
std::remove_reference_t<typename U::value_type>>& element,
std::span<const std::int64_t> nodes, std::span<const std::int64_t> xdofs,
const U& x, int dim,
const U& x, std::array<std::size_t, 2> xshape, int gdim,
std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)>
reorder_fn
= nullptr)
Expand Down Expand Up @@ -430,18 +433,16 @@ create_geometry(
[&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<T> xg(3 * shape0, 0);
for (std::size_t i = 0; i < shape0; ++i)
assert(x.size() % xshape[1] == 0);
std::vector<T> xg(3 * xshape[0], 0);
for (std::size_t i = 0; i < xshape[0]; ++i)
{
std::copy_n(std::next(x.cbegin(), shape1 * l2l[i]), shape1,
std::copy_n(std::next(x.cbegin(), xshape[1] * l2l[i]), xshape[1],
std::next(xg.begin(), 3 * i));
}

return Geometry(dof_index_map, std::move(dofmaps.front()), {element},
std::move(xg), dim, std::move(igi));
std::move(xg), gdim, std::move(igi));
}

} // namespace dolfinx::mesh
39 changes: 22 additions & 17 deletions cpp/dolfinx/mesh/generation.h
Original file line number Diff line number Diff line change
Expand Up @@ -42,13 +42,13 @@ create_interval_cells(std::array<T, 2> p, std::int64_t n);

template <std::floating_point T>
Mesh<T> build_tri(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
std::array<std::int64_t, 2> n,
std::array<std::int64_t, 2> n, int gdim,
const CellPartitionFunction& partitioner,
DiagonalType diagonal);

template <std::floating_point T>
Mesh<T> build_quad(MPI_Comm comm, const std::array<std::array<T, 2>, 2> p,
std::array<std::int64_t, 2> n,
std::array<std::int64_t, 2> n, int gdim,
const CellPartitionFunction& partitioner);

template <std::floating_point T>
Expand Down Expand Up @@ -161,6 +161,7 @@ Mesh<T> create_box(MPI_Comm comm, std::array<std::array<T, 3>, 2> p,
/// @param[in] p Bottom-left and top-right corners of the rectangle.
/// @param[in] n Number of cells in each direction.
/// @param[in] celltype Cell shape.
/// @param[in] gdim Geometric dimension of ambient space.
/// @param[in] partitioner Partitioning function for distributing cells
/// across MPI ranks.
/// @param[in] diagonal Direction of diagonals
Expand All @@ -169,6 +170,7 @@ template <std::floating_point T = double>
Mesh<T> create_rectangle(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
std::array<std::int64_t, 2> n, CellType celltype,
CellPartitionFunction partitioner,
int gdim = 2,
DiagonalType diagonal = DiagonalType::right)
{
if (std::ranges::any_of(n, [](auto e) { return e < 1; }))
Expand All @@ -186,9 +188,9 @@ Mesh<T> create_rectangle(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
switch (celltype)
{
case CellType::triangle:
return impl::build_tri<T>(comm, p, n, partitioner, diagonal);
return impl::build_tri<T>(comm, p, n, gdim, partitioner, diagonal);
case CellType::quadrilateral:
return impl::build_quad<T>(comm, p, n, partitioner);
return impl::build_quad<T>(comm, p, n, gdim, partitioner);
default:
throw std::runtime_error("Generate rectangle mesh. Wrong cell type");
}
Expand All @@ -206,14 +208,16 @@ Mesh<T> create_rectangle(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
/// @param[in] p Two corner points
/// @param[in] n Number of cells in each direction
/// @param[in] celltype Cell shape
/// @param[in] gdim Geometric dimension of ambient space
/// @param[in] diagonal Direction of diagonals
/// @return Mesh
template <std::floating_point T = double>
Mesh<T> create_rectangle(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
std::array<std::int64_t, 2> n, CellType celltype,
int gdim = 2,
DiagonalType diagonal = DiagonalType::right)
{
return create_rectangle<T>(comm, p, n, celltype, nullptr, diagonal);
return create_rectangle<T>(comm, p, n, celltype, nullptr, gdim, diagonal);
}

/// @brief Interval mesh of the 1D line `[a, b]`.
Expand All @@ -225,12 +229,13 @@ Mesh<T> create_rectangle(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
/// @param[in] comm MPI communicator to build the mesh on.
/// @param[in] n Number of cells.
/// @param[in] p End points of the interval.
/// @param[in] gdim Geometric dimension of ambient space
/// @param[in] ghost_mode ghost mode of the created mesh, defaults to none
/// @param[in] partitioner Partitioning function for distributing cells
/// across MPI ranks.
/// @return A mesh.
template <std::floating_point T = double>
Mesh<T> create_interval(MPI_Comm comm, std::int64_t n, std::array<T, 2> p,
Mesh<T> create_interval(MPI_Comm comm, std::int64_t n, std::array<T, 2> p, int gdim = 1,
mesh::GhostMode ghost_mode = mesh::GhostMode::none,
CellPartitionFunction partitioner = nullptr)
{
Expand All @@ -254,12 +259,12 @@ Mesh<T> create_interval(MPI_Comm comm, std::int64_t n, std::array<T, 2> p,
{
auto [x, cells] = impl::create_interval_cells<T>(p, n);
return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x,
{x.size(), 1}, partitioner);
{x.size(), 1}, gdim, partitioner);
}
else
{
return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL,
std::vector<T>{}, {0, 1}, partitioner);
std::vector<T>{}, {0, 1}, gdim, partitioner);
}
}

Expand Down Expand Up @@ -383,7 +388,7 @@ Mesh<T> build_tet(MPI_Comm comm, MPI_Comm subcomm,
}

return create_mesh(comm, subcomm, cells, element, subcomm, x,
{x.size() / 3, 3}, partitioner);
{x.size() / 3, 3}, 3, partitioner);
}

template <std::floating_point T>
Expand Down Expand Up @@ -426,7 +431,7 @@ mesh::Mesh<T> build_hex(MPI_Comm comm, MPI_Comm subcomm,
}

return create_mesh(comm, subcomm, cells, element, subcomm, x,
{x.size() / 3, 3}, partitioner);
{x.size() / 3, 3}, 3, partitioner);
}

template <std::floating_point T>
Expand Down Expand Up @@ -473,12 +478,12 @@ Mesh<T> build_prism(MPI_Comm comm, MPI_Comm subcomm,
}

return create_mesh(comm, subcomm, cells, element, subcomm, x,
{x.size() / 3, 3}, partitioner);
{x.size() / 3, 3}, 3, partitioner);
}

template <std::floating_point T>
Mesh<T> build_tri(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
std::array<std::int64_t, 2> n,
std::array<std::int64_t, 2> n, int gdim,
const CellPartitionFunction& partitioner,
DiagonalType diagonal)
{
Expand Down Expand Up @@ -624,18 +629,18 @@ Mesh<T> build_tri(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
}

return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x,
{x.size() / 2, 2}, partitioner);
{x.size() / 2, 2}, gdim, partitioner);
}
else
{
return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL,
std::vector<T>{}, {0, 2}, partitioner);
std::vector<T>{}, {0, 2}, gdim, partitioner);
}
}

template <std::floating_point T>
Mesh<T> build_quad(MPI_Comm comm, const std::array<std::array<T, 2>, 2> p,
std::array<std::int64_t, 2> n,
std::array<std::int64_t, 2> n, int gdim,
const CellPartitionFunction& partitioner)
{
fem::CoordinateElement<T> element(CellType::quadrilateral, 1);
Expand Down Expand Up @@ -673,12 +678,12 @@ Mesh<T> build_quad(MPI_Comm comm, const std::array<std::array<T, 2>, 2> p,
}

return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x,
{x.size() / 2, 2}, partitioner);
{x.size() / 2, 2}, gdim, partitioner);
}
else
{
return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL,
std::vector<T>{}, {0, 2}, partitioner);
std::vector<T>{}, {0, 2}, gdim, partitioner);
}
}
} // namespace impl
Expand Down
17 changes: 11 additions & 6 deletions cpp/dolfinx/mesh/utils.h
Original file line number Diff line number Diff line change
Expand Up @@ -776,6 +776,7 @@ compute_incident_entities(const Topology& topology,
/// the process offset on`comm`, The offset is the sum of `x` rows on
/// all processed with a lower rank than the caller.
/// @param[in] xshape Shape of the `x` data.
/// @param[in] gdim Geometric dimension of ambient space.
/// @param[in] partitioner Graph partitioner that computes the owning
/// rank for each cell. If not callable, cells are not redistributed.
/// @return A mesh distributed on the communicator `comm`.
Expand All @@ -785,6 +786,7 @@ Mesh<typename std::remove_reference_t<typename U::value_type>> create_mesh(
const fem::CoordinateElement<
typename std::remove_reference_t<typename U::value_type>>& element,
MPI_Comm commg, const U& x, std::array<std::size_t, 2> xshape,
int gdim,
const CellPartitionFunction& partitioner)
{
CellType celltype = element.cell_shape();
Expand Down Expand Up @@ -909,7 +911,7 @@ Mesh<typename std::remove_reference_t<typename U::value_type>> create_mesh(

// Create geometry object
Geometry geometry
= create_geometry(topology, element, nodes1, cells1, coords, xshape[1]);
= create_geometry(topology, element, nodes1, cells1, coords, xshape, gdim);

return Mesh(comm, std::make_shared<Topology>(std::move(topology)),
std::move(geometry));
Expand Down Expand Up @@ -946,6 +948,7 @@ Mesh<typename std::remove_reference_t<typename U::value_type>> create_mesh(
/// the process offset on`comm`, The offset is the sum of `x` rows on
/// all processed with a lower rank than the caller.
/// @param[in] xshape Shape of the `x` data.
/// @param[in] gdim Geometric dimension of ambient space.
/// @param[in] partitioner Graph partitioner that computes the owning
/// rank for each cell. If not callable, cells are not redistributed.
/// @return A mesh distributed on the communicator `comm`.
Expand All @@ -956,6 +959,7 @@ Mesh<typename std::remove_reference_t<typename U::value_type>> create_mesh(
const std::vector<fem::CoordinateElement<
typename std::remove_reference_t<typename U::value_type>>>& elements,
MPI_Comm commg, const U& x, std::array<std::size_t, 2> xshape,
int gdim,
const CellPartitionFunction& partitioner)
{
assert(cells.size() == elements.size());
Expand Down Expand Up @@ -1138,7 +1142,7 @@ Mesh<typename std::remove_reference_t<typename U::value_type>> create_mesh(

// Create geometry object
Geometry geometry
= create_geometry(topology, elements, nodes1, nodes2, coords, xshape[1]);
= create_geometry(topology, elements, nodes1, nodes2, coords, xshape, gdim);

return Mesh(comm, std::make_shared<Topology>(std::move(topology)),
std::move(geometry));
Expand All @@ -1159,21 +1163,22 @@ Mesh<typename std::remove_reference_t<typename U::value_type>> create_mesh(
/// @param[in] elements Coordinate elements for the cells.
/// @param[in] x Geometry data ('node' coordinates). See ::create_mesh
/// for a detailed description.
/// @param[in] xshape The shape of `x`. It should be `(num_points, gdim)`.
/// @param[in] xshape The shape of `x`.
/// @param[in] gdim Euclidean dimension of ambient space.
/// @param[in] ghost_mode The requested type of cell ghosting/overlap
/// @return A mesh distributed on the communicator `comm`.
template <typename U>
Mesh<typename std::remove_reference_t<typename U::value_type>>
create_mesh(MPI_Comm comm, std::span<const std::int64_t> cells,
const fem::CoordinateElement<
std::remove_reference_t<typename U::value_type>>& elements,
const U& x, std::array<std::size_t, 2> xshape, GhostMode ghost_mode)
const U& x, std::array<std::size_t, 2> xshape, int gdim, GhostMode ghost_mode)
{
if (dolfinx::MPI::size(comm) == 1)
return create_mesh(comm, comm, cells, elements, comm, x, xshape, nullptr);
return create_mesh(comm, comm, cells, elements, comm, x, xshape, gdim, nullptr);
else
{
return create_mesh(comm, comm, cells, elements, comm, x, xshape,
return create_mesh(comm, comm, cells, elements, comm, x, xshape, gdim,
create_cell_partitioner(ghost_mode));
}
}
Expand Down
2 changes: 1 addition & 1 deletion cpp/dolfinx/refinement/refine.h
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ refine(const mesh::Mesh<T>& mesh,

mesh::Mesh<T> mesh1 = mesh::create_mesh(
mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(),
mesh.comm(), new_vertex_coords, xshape, partitioner);
mesh.comm(), new_vertex_coords, xshape, mesh.geometry().dim(), partitioner);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn’t xshape[1] and mesh.geometry().dim() be the same?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes - but in terms of the code being correct in meaning, the former is the shape of the geometry data and the latter is the geometric dimension.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we check if they match? We should probably add some sanity checks.


// Report the number of refined cells
const int D = topology->dim();
Expand Down
4 changes: 2 additions & 2 deletions cpp/test/mesh/distributed_mesh.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ constexpr int N = 8;
// Create mesh using all processes and save xdmf
auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
auto mesh = std::make_shared<mesh::Mesh<double>>(
mesh::create_rectangle(comm, {{{0.0, 0.0}, {1.0, 1.0}}}, {N, N},
mesh::create_rectangle(comm, {{{0.0, 0.0}, {1.0, 1.0}}}, {N, N},
mesh::CellType::triangle, part));

// Save mesh in XDMF format
Expand Down Expand Up @@ -139,7 +139,7 @@ void test_distributed_mesh(mesh::CellPartitionFunction partitioner)

// Build mesh
mesh::Mesh mesh = mesh::create_mesh(comm, subset_comm, cells, cmap, comm, x,
xshape, partitioner);
xshape, xshape[1], partitioner);
auto t = mesh.topology();
int tdim = t->dim();
CHECK(t->index_map(tdim)->size_global() == 2 * N * N);
Expand Down
Loading
Loading