Cache cell-facet topological data on meshes - #5298
Conversation
…neighbours and cell-to-cell facet coordinate transforms
…l facet ID to store facets instead of the facet's position in plex cone of the cell
| plex = self.topology_dm | ||
| cStart, cEnd = plex.getHeightStratum(0) # range of DMPlex point numbers representing cells | ||
|
|
||
| for c_plex_point in range(cStart, cEnd): |
There was a problem hiding this comment.
We want to avoid Python loops over the mesh. This should go in some Cython.
| # NOTE: A fresh Dat has its halo entries marked as invalid -> next with_halos access can trigger a collective exchange. | ||
| # An exchange can overrite a rank's data with the owning rank's local data which is wrong. | ||
| # We avoid that by marking halos as valid. | ||
| cell_facet_neighbours.halo_valid = True |
There was a problem hiding this comment.
Does this mean that the Dat is logically over COMM_SELF? If so then that's what we should do. That may be a pyop3 thing (definitely easy there, not sure with PyOP2)
| V = functionspaceimpl.WithGeometry(coordinates.function_space(), self) | ||
| self._coordinates_function = function.Function(V, val=coordinates) | ||
|
|
||
| self._topology_version = 0 |
There was a problem hiding this comment.
Rogue entry from another PR?
| cid = self._cell_numbering.getOffset(c_plex) | ||
|
|
||
| if cid < 0 or cid >= self.cell_set.total_size: | ||
| continue |
There was a problem hiding this comment.
when does this happen?
| ext_plex_points = frozenset() | ||
|
|
||
| cstart, cend = plex.getHeightStratum(0) | ||
| for c_plex in range(cstart, cend): |
There was a problem hiding this comment.
Same comment about Cython
| topology = FIAT.ufc_cell(self.ufl_cell()).get_topology() | ||
| return sum(len(topology[d]) for d in range(facet_dim)) | ||
|
|
||
| def _cell_facet_point(self, cell, local_facet): |
There was a problem hiding this comment.
You are calling this in a tight loop so probably best to not dispatch to a Python call
|
|
||
| for lf in range(num_facets): | ||
| f_point = self._cell_facet_point(cid, lf) | ||
| mask[cid, lf] = f_point in ext_plex_points |
There was a problem hiding this comment.
I wonder if this could all be numpified, instead of cythonified. Use np.intersect1d between a slice of cell_closure and exterior_facets.getStratumIS(1).
|
|
||
| # inverse coord. map | ||
| # o0 -> canonical orientation | ||
| Q0_inv = np.linalg.inv(Q0) |
There was a problem hiding this comment.
Is this not really slow? I think you're calling this for every facet.
There was a problem hiding this comment.
you can batch these things as done in https://github.com/firedrakeproject/firedrake/pull/5034/changes#diff-e45e704e2c1f2f7f0bd632f90a71ce568a47c880cdf097394d8bdf4ccb0332b6R104
However, the geometric computation done there is avoidable and there is a purely topological way of achieving the goal that we inteded there (permuting an array of coordinates). Perhaps you can achieve your goal (which is not clear from the PR description) via a purely topological approach, in a coordinate-free way.
|
|
||
| return (A, b, A_inv) | ||
|
|
||
| def _get_facet_orientation_coord_maps(self): |
There was a problem hiding this comment.
This only uses self in one location, and only to get the UFL cell. This seems like it doesn't need to be a method of the mesh.
I wonder if we should have a firedrake/facets.py file to collect these routines. It could then be folded into the mesh subpackage when we get around to doing that.
There was a problem hiding this comment.
This should definitely live in FIAT
There was a problem hiding this comment.
Specifically FIAT/reference_element.py
| (A[lf], b[lf]) defines the affine embedding of facet reference coordinates into the cell reference frame, | ||
| and A_inv[lf] defines the pseudoinverse for mapping cell coordinates on the facet back to facet-local coordinates. | ||
| """ | ||
| import FIAT |
| and A_inv[lf] defines the pseudoinverse for mapping cell coordinates on the facet back to facet-local coordinates. | ||
| """ | ||
| import FIAT | ||
| ref_cell = FIAT.ufc_cell(self.ufl_cell()) |
There was a problem hiding this comment.
I think these functions live in FIAT, these are very fundamental operations at the reference elemenet level
| return A, b, A_inv | ||
|
|
||
| for lf in range(num_facets): | ||
| phi = ref_cell.get_entity_transform(facet_dim, lf) |
There was a problem hiding this comment.
Is this reverse-engineering get_entity_transform? I.e. are you assuming an affine mapping for the entity transform and rederiving by sampling it at the vertices?
There was a problem hiding this comment.
Tensor product cells will not give you affine mappings. This has to be carefully thought. What is the goal/motivation for storing the facet data on the mesh?
There was a problem hiding this comment.
I store the facet data on the mesh as my particle trajectory algorithm solves for mesh entity collisions (identifies which facet was crossed by each particle and its crossing position on that facet). I would then like to be able to tell, based on the crossed facet ID: 1) which is the neighbouring cell across that facet and 2) transform the coordinates of the point on that facet to the local (reference) coordinates of the "next cell" across that facet.
There was a problem hiding this comment.
I see. Often it helps to follow the existing patterns established in Firedrake. In this case the main pattern to consider is code generation. When we generate code to assemble a form, we don't precompute/store in the mesh the geometric information (Jacobians/affine mappings). We compute everything from the mesh coordinates withtin the generated C kernels (avoiding python loops over the cells).
It seems to me that the right thing is to generate a facet kernel. A facet kernel iterates over pairs of cells sharing a facet, this solves 1). You can then access the coordinates of both cells and generate code from ufl.FacetJacobian(mesh), this solves 2).
Generation of bespoke code is quite painful, but we are here to help you. I think by sticking to the established patterns one can get a very general particle code that works on high-order meshes of any cell type.
There was a problem hiding this comment.
Yeah, I completely understand your point. As this idea of writing kernels to generate code is pretty new to me, could you point out where I could look for these code patterns for example in existing kernels that achieve some sort of similar operation (e.g., the form assembly you mentioned)? It will help me get started on this task!
3 participants
Compute and cache properties derived from cell-facet adjacency data on meshes:
cell_facet_neighbours: maps each cell to its neighbouring cells across each of its local facetscell_facet_coord_transforms: stores affine reference-coordinate transforms between neighbouring cells across interior facetscell_facet_exterior_mask: marks cell-local facets lying on the exterior domain boundary (distinguishes exterior facets from rank partition boundaries)