dolfinx icon indicating copy to clipboard operation
dolfinx copied to clipboard

Orient mesh entities consistently and generalise `entities_to_geometry`

Open jpdean opened this issue 9 months ago • 1 comments

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.

jpdean avatar May 13 '24 17:05 jpdean

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: grid_1 Cell plot: grid_2

jorgensd avatar May 16 '24 13:05 jorgensd