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.

Sign up for GitHub

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

Jump to bottom

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

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

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
17c2b90
Allow promotion of unit interval geometric dimension.
jhale Oct 14, 2024
f399466
Syntax fix.
jhale Oct 14, 2024
de08a02
Better docstring.
jhale Oct 14, 2024
90a8fe9
Update mesh.py
jhale Oct 14, 2024
afb1cf1
Update mesh.py
jhale Oct 14, 2024
7c74a5b
Fix formatting.
jhale Oct 15, 2024
ba428ca
Use argument names for optional arguments.
jhale Oct 15, 2024
a061aee
Push out to unit square, other small fixes/tweaks.
jhale Oct 15, 2024
13eade3
Assertion error as dim and gdim are different.
jhale Oct 15, 2024
a7f587e
Propagate xshape through to Geometry
jhale Oct 15, 2024
3191823
Ruff format.
jhale Oct 15, 2024
f909b00
Fixup mesh tests.
jhale Oct 15, 2024
c2e5b31
Partial fixes for C++ unit tests.
jhale Oct 15, 2024
e155b4c
Update distributed_mesh.cpp
jhale Oct 15, 2024
aa9e35c
Add back assertion, does no harm.
jhale Oct 16, 2024
3459155
Debugging.
jhale Oct 16, 2024
2eea134
Merge branch 'jhale/add-gdim-to-mesh-creators' of github.com:FEniCS/d…
jhale Oct 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion cpp/dolfinx/io/XDMFFile.cpp
View file
Open in desktop
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
32 changes: 14 additions & 18 deletions cpp/dolfinx/mesh/Geometry.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,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 +267,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 +339,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 +371,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 +383,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 +429,15 @@ 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.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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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
8 changes: 4 additions & 4 deletions cpp/test/mesh/generation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ TEMPLATE_TEST_CASE("Interval mesh (parallel)", "[mesh][interval]", float,
};

mesh::Mesh<T> mesh = mesh::create_interval<T>(
MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.}, ghost_mode, part);
MPI_COMM_WORLD, 5 * comm_size - 1, {0., 1.}, 1, ghost_mode, part);

{
int comp_result;
Expand Down Expand Up @@ -304,7 +304,7 @@ TEMPLATE_TEST_CASE("Rectangle triangle mesh (right)",
using T = TestType;

mesh::Mesh<T> mesh = dolfinx::mesh::create_rectangle<T>(
MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle,
MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle, 2,
mesh::DiagonalType::right);

// vertex layout:
Expand Down Expand Up @@ -345,7 +345,7 @@ TEMPLATE_TEST_CASE("Rectangle triangle mesh (left)",
using T = TestType;

mesh::Mesh<T> mesh = dolfinx::mesh::create_rectangle<T>(
MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle,
MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle, 2,
mesh::DiagonalType::left);

// vertex layout:
Expand Down Expand Up @@ -387,7 +387,7 @@ TEMPLATE_TEST_CASE("Rectangle triangle mesh (crossed)",
using T = TestType;

mesh::Mesh<T> mesh = dolfinx::mesh::create_rectangle<T>(
MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle,
MPI_COMM_SELF, {{{0, 0}, {1, 1}}}, {1, 1}, mesh::CellType::triangle, 2,
mesh::DiagonalType::crossed);

// vertex layout:
Expand Down
Loading
Loading