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

Mesh preprocessing operations #260

Merged
merged 14 commits into from
Jul 1, 2024
Merged

Mesh preprocessing operations #260

merged 14 commits into from
Jul 1, 2024

Conversation

sebastiangrimberg
Copy link
Contributor

@sebastiangrimberg sebastiangrimberg commented Jun 5, 2024

  • Convert all elements to simplices (enable with "MakeSimplex")
  • Convert all elements to hexahedra (enable with "MakeHexahedral", only for tet meshes currently, compatible with "MakeSimplex")
  • Serial refinement levels ("SerialUniformLevels")

All mesh conversions work with high-order curved meshes by interpolating the nodal grid function between meshes using GSLIB. Also exposes via configuration file a number of mesh preprocessing options which used to be hard-coded.

@sebastiangrimberg sebastiangrimberg added enhancement New feature or request mesh Related to meshes and mesh generation labels Jun 5, 2024
@sebastiangrimberg sebastiangrimberg force-pushed the sjg/mesh-simplex-dev branch 2 times, most recently from 24fc0d8 to 77cbc97 Compare June 21, 2024 18:44
Copy link
Collaborator

@hughcars hughcars left a comment

Choose a reason for hiding this comment

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

A few minor questions and comments, but LGTM.

int npts = 0, offset = 0;
for (int i = 0; i < dest_mesh.GetNE(); i++)
{
npts += V.FESpace()->GetFE(i)->GetNodes().GetNPoints();
Copy link
Collaborator

Choose a reason for hiding this comment

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

All points are interior to elements or if there is duplication at intersections we don't mind? I would guess it makes indexing etc way more complicated so probably not worth it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah the latter, we just find the duplicate shared points nodes between elements and they are disregarded because we SetSubVector later on. Certainly a future optimization could be to remove these duplicates from the GSLIB search and just maintain the indexing separately.

@@ -212,7 +285,7 @@ void RefineMesh(const IoData &iodata, std::vector<std::unique_ptr<mfem::ParMesh>
// reorient only the coarse mesh so that the refinements are still true refinements of
// the original mesh (required for geometric multigrid). Otherwise, it happens after
// refinement.
if (iodata.model.reorient_tet && mesh.capacity() > 1)
if (iodata.model.reorient_tet_mesh && mesh.capacity() > 1)
Copy link
Collaborator

Choose a reason for hiding this comment

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

If mesh.empty() == true this could still be triggered as capacity is the largest value seen. I think not an actual problem as we're always in the case where size() == capacity(), but worth fixing.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The check at https://github.com/awslabs/palace/pull/260/files#diff-0642b48eb210cdd7e98a30345d4857eb87c192985610bc00696294583cd09a43R259 should prevent that from ever being the case. The idea here is that if we are storing a multigrid hierarchy of meshes, we want to reorient the coarse mesh, otherwise we will wait until the end to do the reorientation.

@@ -360,7 +433,7 @@ void RefineMesh(const IoData &iodata, std::vector<std::unique_ptr<mfem::ParMesh>
// not required after MFEM's PR #1046, but in case you want to be absolutely sure, we
// reorient only the mesh after refinement if there is a single mesh (doesn't work with
// h-refinement geometric multigrid).
if (iodata.model.reorient_tet && mesh.size() == 1)
if (iodata.model.reorient_tet_mesh && mesh.capacity() == 1)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Changing from size to capacity, implies a vector where pop() etc. had been called shouldn't on this branch, is that intentional?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

See my comment above, but in particular the change here is just to be consistent in use of capacity for both checks.

Comment on lines +87 to +91
if (!use_mesh_partitioner)
{
MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, Mpi::Rank(comm), MPI_INFO_NULL,
&node_comm);
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

To check I am getting this correctly, the idea is to split the communicator into N communicators, one for each node. Then within each communicator, only the root rank loads the serial, thereby avoiding the memory overloading issue of every rank on a node loading the mesh for preprocessing purposes.

Copy link
Contributor Author

@sebastiangrimberg sebastiangrimberg Jun 26, 2024

Choose a reason for hiding this comment

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

Correct. Just one process on each node does the serial mesh preprocessing operations to save memory overhead, and then broadcasts the final serial mesh to all other processes on the node.

palace/utils/geodata.cpp Show resolved Hide resolved
palace/utils/geodata.cpp Outdated Show resolved Hide resolved
palace/utils/geodata.cpp Show resolved Hide resolved
Comment on lines +1862 to +1978
if (end == start + 1)
{
new_fespace.GetElementVDofs(start, vdofs);
new_nodes.SetSubVector(vdofs, vals);
start = end;
continue;
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

In the MakeSimplex and MakeHex methods, neither allow mixed presently, so is this an edge case that can ever be hit? I.e. all tets get converted, and all hexes get converted. Is this just anticipating mixed meshes?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually MakeSimplex allows for a mixed mesh, so this code gets hit for a mixed tet-prism mesh, for example.

Comment on lines 1882 to 1999
for (int j = 0; j < pointmat.Width(); j++)
{
for (int d = 0; d < sdim; d++)
{
// Use default ordering byNODES.
xyz(d * npts + offset + j) = pointmat(d, j);
}
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

For this loop, rather than jumping back and forwards by npts in xyz for each write, it might be more cache coherent to have the d loop on the outer. This loops the slow way through pointmat (assuming column major), but pointmat is smaller than xyz, so the npt jumps about in xyz might be causing more reads etc.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah good call about pointmat being small so the jumps not being a big issue. Fixed.

@sebastiangrimberg sebastiangrimberg merged commit 385804b into main Jul 1, 2024
17 checks passed
@sebastiangrimberg sebastiangrimberg deleted the sjg/mesh-simplex-dev branch July 1, 2024 16:16
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request mesh Related to meshes and mesh generation
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants