After some discussion about the difficulty of iterating over various entities within the mesh, we've observed the following things:
It's difficult to create a "TDyMesh" instance, either by itself from parameters or from a DMPlex.
Looping over the various entities (vertices of a cell, subcells of a vertex, and so on) requires a lot of knowledge of the TDyMesh data structure, which makes looping code brittle.
It would be nice if we had a more consistent interface to the mesh, both to simplify its creation (for testing and within the dycore itself) and to access its entities. Here's a proposed interface for a refactored mesh. All functions return the usual error code, with 0 for success and anything else indicating failure.
Creating and Destroying a Mesh
// Create a mesh from a DM object. We need to be careful about specifying preconditions on the DM here.
TDyMeshCreateFromDM(DM dm, TDyMesh* mesh);
// Destructor function for a mesh.
TDyMeshDestroy(TDyMesh* mesh);
We might also offer some convenience functions for box meshes, reading from a file, etc.
Looping Over Entities within a Mesh
// Given a mesh and a cell index, retrieve an array of cell vertex indices, and their number.
TDyMeshGetCellVertices(TDyMesh* mesh, PetscInt cell, PetscInt** vertices, PetscInt* num_vertices);
// Given a mesh and a cell index, retrieve an array of indices of faces bounding the cell,
// and their number.
TDyMeshGetCellFaces(TDyMesh* mesh, PetscInt cell, PetscInt** faces, PetscInt* num_faces);
// Given a mesh and a cell index, retrieve an array of indices of neighboring cells, and their number.
TDyMeshGetCellNeighbors(TDyMesh* mesh, PetscInt cell, PetscInt** neighbors, PetscInt* num_neighbors);
// Retrieve the centroid for the given cell in the given mesh.
TDyMeshGetCellCentroid(TDyMesh* mesh, PetscInt cell, TDyCoordinate* centroid);
// Retrieve the volume for the given cell in the given mesh.
TDyMeshGetCellVolume(TDyMesh* mesh, PetscInt cell, PetscReal* volume);
// Given a mesh and a subcell index, retrieve an array of indices of faces bounding the subcell,
// and their number.
TDyMeshGetSubcellFaces(TDyMesh* mesh, PetscInt subcell, PetscInt** faces, PetscInt* num_faces)
// Given a mesh and a vertex index, retrieve an array of subcell indices and their number.
TDyMeshGetVertexSubcells(TDyMesh* mesh, int vertex, int** subcells, int* num_subcells);
// Given a mesh and a face index, retrieve an array of vertex indices and their number.
TDyMeshGetFaceVertices(TDyMesh* mesh, int face, int** vertices, int* num_vertices);
This list is probably incomplete, but hopefully you get the point. In most cases, these functions just set a pair of pointers, so they should be cheap.
I believe all of these functions can be draped over the mesh as is, without breaking existing code. Then we can refactor this code to use the new interface.
I think this is finished now. If we adopt the unit testing approach that we've discussed elsewhere (#138, #139), we can factor out a "constructor" function for the mesh as part of #140.
After some discussion about the difficulty of iterating over various entities within the mesh, we've observed the following things:
DMPlex
.TDyMesh
data structure, which makes looping code brittle.It would be nice if we had a more consistent interface to the mesh, both to simplify its creation (for testing and within the dycore itself) and to access its entities. Here's a proposed interface for a refactored mesh. All functions return the usual error code, with 0 for success and anything else indicating failure.
Creating and Destroying a Mesh
We might also offer some convenience functions for box meshes, reading from a file, etc.
Looping Over Entities within a Mesh
This list is probably incomplete, but hopefully you get the point. In most cases, these functions just set a pair of pointers, so they should be cheap.
I believe all of these functions can be draped over the mesh as is, without breaking existing code. Then we can refactor this code to use the new interface.