dolfinx
dolfinx copied to clipboard
Orient mesh entities consistently and generalise `entities_to_geometry`
In main
, the orientation of mesh entities is defined by whichever cell appears last in the entity-to-cell connectivity. This assumption is relied on in a few places in the code (see e.g. entities_to_geometry
), despite the assumption not being well documented.
This PR ensures that mesh entities are created such that their local (to the process) orientation agrees with their global orientation, which is unique, can be computed locally, and is a more natural definition. It will also greatly simplify mixed-dimensional coupling.
This PR also generalises entities_to_geometry
for high-order elements. In main
, this function computes the geometric DOFs associates with the vertices of a particular entity, whereas in this PR, it computes the geometric DOFs associated with the closure of the entities. I've removed the option orient
, which only worked for the facets of tets, but let me know if this is problematic. I've added an option permute
; this permutes the DOFs so that they agree with the orientation of mesh entities, but it requires create_entity_permutations
to be called. In many cases, you don't need to do this (e.g. when computing the diameter of an entity). Alternatively, this could be implemented as a separate function.
Thanks @jorgensd for the idea to deal with high-order geometry via closure DOFs and @mscroggs for help with permuting the DOFs.
Example from: http://jsdokken.com/FEniCS23-tutorial/src/mesh_generation.html now slightly modified:
from mpi4py import MPI
import dolfinx.mesh
import dolfinx.plot
import basix.ufl
import ufl
import numpy as np
import pyvista
pyvista.start_xvfb(wait=0.1)
def plot_mesh(mesh, dim):
mesh.topology.create_connectivity(dim, mesh.topology.dim)
mesh.topology.create_entity_permutations()
vtk_mesh = dolfinx.plot.vtk_mesh(mesh, dim=dim)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
linear_grid = dolfinx.plot.vtk_mesh(V)
lgrid = pyvista.UnstructuredGrid(*linear_grid)
ugrid = pyvista.UnstructuredGrid(*vtk_mesh)
tgrid = ugrid.tessellate()
plotter = pyvista.Plotter()
plotter.add_mesh(ugrid,style="points", color="b", point_size=10, render_points_as_spheres=True)
plotter.add_mesh(tgrid, show_edges=False, opacity=0.5)
plotter.add_mesh(lgrid, style="wireframe", color="b")
plotter.show()
plotter.screenshot(f"grid_{dim}.png")
nodes = np.array([[1., 0.], [2., 0.], [3., 2.],
[2.9, 1.3], [1.5, 1.5], [1.5, -0.2]], dtype=np.float64)
connectivity = np.array([[0, 1, 2, 3, 4, 5]], dtype=np.int64)
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 2, shape=(2, )))
domain = dolfinx.mesh.create_mesh(MPI.COMM_SELF, connectivity, nodes, c_el)
plot_mesh(domain, 1)
plot_mesh(domain, 2)
yields:
Facet plot:
Cell plot: