This code builds triangular unstructured meshes for a 2D rectangular domain in parallel with MPI, designed for use with PETSc. It depends on PETSc configured with triangle (--download-triangle) and returns a parallel DMPlex object which stores the unstructured mesh. To generate large meshes ensure PETSc is configured with 64-bit integers (--with-64-bit-indices).

Domain Size

The domain size can be specified by passing in the -domain_width and -domain_height command line arguments. The default domain size is [0,1] x [0,1].

Motivation

This code was built to enable easy testing of numerical methods (see PFLARE) on large, fully unstructured meshes in parallel. Most freely available meshing tools only feature OpenMP parallelism, require file I/O or require load-balancing with a mesh partitioning tool (such as ParMETIS) which can be very expensive.

At scale, unstructured meshes are often generated by first partitioning a coarse, fully unstructured mesh and then performing structured refinement (e.g., by using -dm_refine in PETSc). This generates an increasingly structured mesh and there are some applications where this is not desirable, such as advection problems, where alignment between advection velocity and mesh faces can affect the results.

This code instead:

• Generates a fully unstructured mesh.
• Uses MPI for parallelism.
• Allows for an in-memory representation of a mesh without file I/O.
• Produces load-balanced meshes without explicitly calling a mesh partitioner such as ParMETIS.

To enable this several compromises were made, namely:

• Can only produce meshes on a simple rectangular domain.
• Uniform resolution throughout the domain.
• Produces good elements (e.g., with reasonable angles and volume ratios) but not necessarily optimal.
• The mesh has reasonably low communication volume, but is not necessarily communication minimising. A mesh partitioner like ParMETIS can be explicitly called by using DMPlexDistribute on the returned DM to further minimise the communication volume.
• Not fully optimised for speed/memory.
• Not robust when the number of elements per MPI rank is small (say <100k).
• If differing numbers of MPI ranks are used, the mesh produced is not identical.

Building the code as an executable

To build an executable which can be called from the command line for small scale testing, ensure PETSC_DIR and PETSC_ARCH environmental variables are set and then call make clean && make.

There are five input variables that can be changed from the command line:

Command line argument	Default value	Details
-target_edge_length	0.0025	Target edge length for elements. Resulting mesh will have edges close to this value.
-final_smooth_its	4	How many iterations of smoothing (LLoyds + springs) to do. More iterations will increase mesh quality and runtime.
-write_mesh	false	Output the PETSc DM to disk in HDF5 format.
-integrity_check	true	Run mesh integrity checks and return NULL if not valid. This takes extra memory and time. Recommend disabling this for production runs.
-print_stats	true	Print global mesh statistics on rank 0. This takes extra memory and time. Recommend disabling this for production runs.

For example, after building the executable we can generate a mesh using 2 MPI ranks on the command line by calling:

 mpiexec -n 2 ./BoxMeshDM

which will generate a mesh with the default parameters. To decrease the edge length for example,

 mpiexec -n 2 ./BoxMeshDM -target_edge_length 0.001

If you wish to write out the meshes generated this way, ensure PETSc has been configured with HDF5 (--download-hdf5) and run the code with -write_mesh true. The resulting .h5 file can be read into a PETSc DMPlex with -dm_plex_filename box_mesh.h5.

To visualise the mesh, from the command line run ${PETSC_DIR}/lib/petsc/bin/petsc_gen_xdmf.py box_mesh.h5. The resulting .xmf file can be visualised in Paraview with the XDMF reader.

Building the code as a library

For large scale use, the code can be compiled as a library. Hence instead of writing out the mesh, the routine GenerateBoxMeshDM can be called directly from existing code. This returns a parallel, load balanced PETSc DMPlex object that can be used without I/O.

Ensure PETSC_DIR and PETSC_ARCH environmental variables are set and then call make clean && make lib. You then need to include the .h file in your code and link to the output library libboxmeshdm.

In your code, to generate a PETSc DM that can then be used as normal, you can call in C/C++ (see test_lib.c):

 #include "BoxMeshDM.h"

 // ...
 
 DM dm;
 PetscErrorCode ierr;

 // Set target mesh edge length
 double target_edge_length = 0.0025;
 // Set the number of smoothing iterations
 PetscInt final_smooth_its = 4;
 // Check the integrity of the mesh and error if not valid
 PetscBool integrity_check = PETSC_TRUE;
 // Print global mesh statistics from MPI rank 0
 PetscBool print_stats = PETSC_TRUE;

 // Generate the unit-box mesh stored in a parallel PETSc DM on the MPI_Comm PETSC_COMM_WORLD
 dm = GenerateBoxMeshDM(PETSC_COMM_WORLD, target_edge_length, 1.0, 1.0, final_smooth_its, integrity_check, print_stats);
 ierr = DMDestroy(&dm);

 // Or specify a rectangular domain explicitly
 dm = GenerateBoxMeshDM(PETSC_COMM_WORLD, target_edge_length, 2.0, 1.5, final_smooth_its, integrity_check, print_stats);

 // Enable the use of command line options for this DM
 ierr = DMSetFromOptions(dm);

Weak scaling

Weak scaling results on ARCHER2 show the code is reasonably performant when generating large meshes in parallel:

Compute nodes	MPI ranks	Target edge length	Elements	Min angle (degrees)	Min/max volume ratio	Load imbalance	Time (s)
256	32768	7.8125e-6	32B	15.7	19.3	1.00131	114
64	8192	1.5625e-5	8.2B	19.4	13.9	1.00133	102
16	2048	3.125e-5	2B	19.4	14.2	1.00126	108
4	512	6.25e-5	512M	21.5	13.6	1.00125	92
1	128	1.25e-4	128M	16.7	11.4	1.00111	82

To improve the mesh quality increase the number of final smoothing iterations. To improve the runtime, disable both the integrity check and statistics printing.