Problem with field extraction on grid


I want to extract field values on a grid following the approach shown at the end of this tutorial:

For some reason, a lot of points are unable to find any colliding cell (and get marked np.nan) which leads to the following result using imshow():


The program is/should be running on only one processor as I did not use mpirun to launch it. I am also quite sure the function was interpolated right as pyvista has no problem visualizing it. What is puzzling is that this behavior seems to be mesh/geometry dependent. On a simple geometry (say square), I get expected results. Only when I use this somewhat skewed geometry which is required for my problem do I run into this issue.

Below is an MWE evaluated on dolfinx 0.5.2 installed through conda. Apologies in advance for the complicated mesh generation function but as I said, a simple geometry is not affected by this anomaly.

import dolfinx, ufl, mpi4py, gmsh, petsc4py, pyvista
from import gmshio
import numpy as np
import matplotlib.pyplot as plt

def create_geometry():
    '''define the solar cell stack. metal blocks are centered at x = fov_x/2. Boundaries are addressed for the mesh element size'''
    gdim = 2
    fov_x = 2e-2
    # stack length in m
    L = np.array([100e-6, # cg, index 0 
                  50e-6, # si, index 1
                  1e-6, # AR, 1um, index 2
                  10e-6, # metal blocks layer, index 3
                  15e-6, # TJC, index 4
                  150e-6, # Ge, index 5
                  10e-6, # metal, index 6
                  50e-6, # Si, index 7
                  0.484e-3, # CFRP, index 8
                  25e-3, # HC core, index 9
                  0.484e-3, # CFRP, index 10 
    met_width = L[3]
    gmsh.clear() # close existing session, if any
    occ = gmsh.model.occ
    L_flat = L.copy()
    L_flat[1] = np.sum(L_flat[1:3])
    L_flat = np.delete(L_flat, (2, 3))
    fov_y = np.sum(L_flat)
    fov = occ.add_rectangle(0, 0, 0, fov_x/2, fov_y)
    layers = []
    for idx, thick in enumerate(L_flat[::-1]):
        layers.append((gdim, occ.add_rectangle(0, sum(L_flat[::-1][:idx]), 0, fov_x/2, thick)))
    # metal contact and ARC
    layers.append((gdim, occ.add_rectangle(0, sum(L_flat[2:]), 0, fov_x/2, L[2]))) # ARC
    layers.append((gdim, occ.add_rectangle(fov_x/2 - met_width/2, sum(L_flat[2:]), 0, met_width/2, L[3]))) # metal contact
    layers.append((gdim, occ.add_rectangle(fov_x/2 - (met_width/2+L[2]), sum(L_flat[2:]), 0, met_width/2+L[2], L[2]+L[3]))) # metal contact + AR
    geom = occ.fragment([(gdim, fov)], layers)

    # tag all domains for later
    all_domains = gmsh.model.getEntities(gdim)
    for j, domain in enumerate(all_domains):
        gmsh.model.addPhysicalGroup(gdim, [domain[1]], j+1)

    boundaries = gmsh.model.getEntities(gdim-1)
    for idx, bnd in enumerate(boundaries):
        gmsh.model.addPhysicalGroup(gdim-1, [bnd[1]], tag=idx+1)

    # generate fine mesh around curved boundaries and course away from them
    # taken from
    surf = [(gdim, 8), (gdim, 6), (gdim, 6), (gdim, 6), ] # by inspection, original gmsh indexing and not the group tag
    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(surf, oriented=False)

    arc = (1, )
    tjc = (39, )
    ge = (36, )
    boundaries = np.hstack((tjc, arc, ge))

    gmsh.model.mesh.field.setNumbers(1, "EdgesList", boundaries)
    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 1e-5)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 1)  # 5 * R_outer
    gmsh.model.mesh.field.setNumber(2, "DistMin", 3e-3)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 25)
#     # # Generate mesh
    gmsh.option.setNumber("Mesh.Algorithm", 7)
    # convert gmsh geometry to fenics geometry
    mpi_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi4py.MPI.COMM_WORLD.rank, gdim)    
    return mesh, ct, ft

mesh, ct, ft = create_geometry()

V = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))
Qs = dolfinx.fem.Function(V)

qs_mu, qs_sigma = (0, 0.0261855), 1e-03 # central location and sigma of the guassian source
qs_func = lambda x : np.exp(-((x[0]-qs_mu[0])**2 + 1*(x[1]-qs_mu[1])**2)/(2*qs_sigma**2))

# extract field in a small region around the gaussian function's peak
xaxis = np.linspace(0, qs_sigma, 200)
yaxis = np.linspace(qs_mu[1]-1e-4, qs_mu[1]+1e-4, 200)
x2, y2 = np.meshgrid(xaxis, yaxis)
points = np.zeros((3, x2.size))
points[0, :], points[1, :] = x2.ravel(), y2.ravel()  # the third axis remains zero

bb_tree = dolfinx.geometry.BoundingBoxTree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = dolfinx.geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates, points.T)

Qs_interp = np.zeros(points.shape[1])
Qs_interp[:] = np.nan
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i)) > 0:
        # points_on_proc.append(point)
        # cells.append(colliding_cells.links(i)[0])
        cell = colliding_cells.links(i)[0]
        Qs_interp[i] = Qs.eval(point, cell)

plt.imshow(np.flipud(Qs_interp.reshape((yaxis.size, xaxis.size))))

The problem here is that you cells are very small, (cell diameter of 1e-7). It is very tricky to have exact collision detection on such elements.
You should preferrably scale your mesh, such that you avoid these issues.

Thanks a lot for your answer. That means if I solve the problem on this mesh then the computations will not suffer due to it, only the postprocessing has limitations. I can manage a life with this situation.