Find adjacent DOF

dolfinx
1

Is it possible to find the adjacent DOFs connected by edges?

In my case, I have a marked domain with a Function of 0 and 1 (Fig a). I try to find DOFs on the exterior layer of a marked region (Fig b). I pick all those points by looping all DOFs with 0 and at least 1 adjacent DOF with 1.

To do so, I made a specific case where P1 element is used. Dim 0 vertex and Dim 1 edge connectivity are created, and I examine all connected vertices based on the current vertex. This would be insufficient in other cases, e.g. high order element and mixed element. Any advice for improvement?

adjacent_DOFs
adjacent_DOFs1434×622 364 KB

MWE:

from mpi4py import MPI
from dolfinx import fem, mesh
import numpy as np

def mark(x):
    dist = np.sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2)
    values = np.ones_like(dist, dtype=np.int32)
    values[dist <= 0.25] = 0
    return values

comm = MPI.COMM_WORLD
_mesh = mesh.create_unit_square(comm, 64, 64, mesh.CellType.quadrilateral)
V = fem.functionspace(_mesh, ("CG", 1))
f = fem.Function(V)
f.interpolate(mark)

geometry = _mesh.geometry
topology = _mesh.topology
num_vertices = topology.index_map(0).size_local + topology.index_map(0).num_ghosts

vertices = []

topology.create_entities(1)
topology.create_connectivity(0, 1)
vertex2edge = topology.connectivity(0, 1)
edge2vertex = topology.connectivity(1, 0)
for vertex in range(num_vertices):
    adjacent_vertices = []
    edges = vertex2edge.links(vertex)
    for edge in edges:
        found_vertices = edge2vertex.links(edge)
        adjacent_vertices.append(list(found_vertices))
    
    adjacent_vertices = np.unique(adjacent_vertices)
    adjacent_vertices_excluded = np.setdiff1d(adjacent_vertices, vertex)
    
    if (f.x.array[vertex] == 0) and np.any(f.x.array[adjacent_vertices_excluded] == 1):
        vertices.append(vertex)

# Validation
import matplotlib.pyplot as plt
plt.figure()
for vertex in vertices:
    x, y, z = geometry.x[vertex, :]
    plt.plot(x, y, 'bo')
plt.gca().set_aspect(1)
1 Like
2

Ti me it seems like you want to use locate_dofs_topological on a subset of facets that you have found given some criteria.

However, is is such that your use-case has to go from a function (with 0,1 values) to the dofs? Or do you have some other application/use-case in mind?

3

Indeed, the above approach is my unmature way in mind to deal with the real problem. I am also looking for a more sophisticated method.

In my FSI application, it involves two meshes: the foreground solid mesh and the background fluid mesh (Fig a). The shape of foreground mesh could be not that regular.

To increase the accuracy in the system between two sets of meshes, I want to pick either the vertices or DOFs on the most exterior layer of background mesh in the overlapped region to either
(1) apply boundary condition for solving an equation to move mesh to fit the foreground mesh shape,
or
(2) change the interpolation process (e.g. L2 projection) only for the most exterior layer of mesh, for instance, only use the overlapped region data to interpolate onto foreground mesh.

So, I design a function on two meshes. And then interpolate foregound one onto background one to get a 0-1 mask function (Fig b). By looking up the situation of current DOF and its adjacent DOFs to pick up the most exterior layer of DOFs.

marked_mask_mesh
marked_mask_mesh1254×669 93.3 KB

4

As you are dealing with FSI, are we only talking Lagrange/DG function spaces, where a dof is associated with a coordinate (that being a vertex or a point on a facet, edge or internally in the cell)? If so you can use compute_closest_entity on the dof coordinates of the foreground mesh: dolfinx/python/dolfinx/geometry.py at 58fd47824ab920fd6576e4d321782c683b168532 · FEniCS/dolfinx · GitHub

5

Yes, Lagrange/DG functionspace is all I need.

I give it a go as you adviced: I set up two overlapped square meshes as fluid and solid. I build bb_tree and midpoint_tree based on background fluid mesh, then find the DOF coordinates of exterior boundary of solid mesh via locate_dofs_topological.

So far, I only found all the colliding cells (cell ID) in a serial run, and it is like a geometry.compute_colliding_cells. For a colliding cell, how can I locate only the interior point A but not exterior point B as pictured?

interior_exterior
interior_exterior899×557 14.5 KB

MWE with DOLFINx v0.8.0:

from mpi4py import MPI
from dolfinx import fem, mesh, geometry
import numpy as np

comm = MPI.COMM_WORLD

# Fluid
fluid_mesh = mesh.create_unit_square(comm, 64, 64)
V = fem.functionspace(fluid_mesh, ("CG", 1))
f = fem.Function(V)

_geometry = fluid_mesh.geometry
_topology = fluid_mesh.topology
tdim = _topology.dim
fdim = tdim - 1
num_vertices = _topology.index_map(0).size_local + _topology.index_map(0).num_ghosts

# Solid
solid_mesh = mesh.create_unit_square(comm, 16, 16)
solid_mesh.geometry.x[:, :tdim] -= 0.5
solid_mesh.geometry.x[:, :tdim] *= 0.49
solid_mesh.geometry.x[:, :tdim] += 0.5

Vs = fem.functionspace(solid_mesh, ("CG", 1))
fs = fem.Function(Vs)

num_entities = _topology.index_map(tdim).size_local + _topology.index_map(tdim).num_ghosts
entities = np.arange(num_entities, dtype=np.int32)
tree = geometry.bb_tree(fluid_mesh, tdim)
midpoint_tree = geometry.create_midpoint_tree(fluid_mesh, tdim, entities)

# DOF on the boundary of foreground solid mesh
solid_mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(solid_mesh.topology)
boundary_dofs = fem.locate_dofs_topological(Vs, fdim, boundary_facets)
boundary_dof_coordinates = Vs.tabulate_dof_coordinates()[boundary_dofs]
closest_entities = geometry.compute_closest_entity(tree, midpoint_tree, 
                                                   fluid_mesh,
                                                   boundary_dof_coordinates)

# Validation
import matplotlib.pyplot as plt
plt.figure()
plt.gca().set_aspect(1)
plt.title(f"Process {comm.rank}")
for coordinate in boundary_dof_coordinates:
    x, y = list(coordinate)[:2]
    plt.plot(x, y, 'b.')

for cell in _geometry.dofmap[closest_entities]:
    xs, ys = [], []
    for vertex in _geometry.x[cell]:
        xs.append(vertex[0])
        ys.append(vertex[1])
    
    xs.append(xs[0])
    ys.append(ys[0])
    plt.plot(xs, ys, 'r')

plt.show()
6

If you want your code to work in parallel, I would simply get all nodes of your mesh, i.e. mesh.geometry.x and pass it to: dolfinx/cpp/dolfinx/geometry/utils.h at 58fd47824ab920fd6576e4d321782c683b168532 · FEniCS/dolfinx · GitHub
(it has a Python wrapper, see for instance: Can you help me apply compute_cell_contributions method to 1D cantilever beam example? - #6 by dokken for an example.

It will give you the cell that collides with every point (vertex) of your mesh (and tell you what process it came from.

7

Indeed, I used cpp.geometry.determine_point_ownership for my Lagrange-only parallel interpolation matrix generation instead of create_nonmatching_meshes_interpolation_data.

With it in a serial run, I can pick the interior DOFs in the colliding cells as

from dolfinx import cpp

dofs = []
for cell in closest_entities:
    dofs.append(list(V.dofmap.cell_dofs(cell)))
dofs = np.unique(dofs).astype(np.int32)
dof_coords = _geometry.x[dofs]
overlaps = cpp.geometry.determine_point_ownership(solid_mesh._cpp_object, dof_coords, 1e-6)[0]

for dof, overlap in zip(dofs, overlaps):
    if overlap >= 0:
        coords = _geometry.x[dof]
        x, y = list(coords)[:2]
        plt.plot(x, y, "k*")

yielding black * in the figure:

image
image792×729 26.1 KB