Hi,
I want to extract field values on a grid following the approach shown at the end of this tutorial:
https://jorgensd.github.io/dolfinx-tutorial/chapter1/membrane_code.html
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 dolfinx.io 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.initialize()
gmsh.option.setNumber("General.Terminal",0)
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
occ.synchronize()
geom = occ.fragment([(gdim, fov)], layers)
occ.synchronize()
# 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 https://github.com/jorgensd/dolfinx-tutorial/blob/v0.4.0/chapter3/em.ipynb
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)
gmsh.model.mesh.field.setAsBackgroundMesh(2)
# # # Generate mesh
gmsh.option.setNumber("Mesh.Algorithm", 7)
gmsh.model.mesh.generate(gdim)
# 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)
gmsh.finalize()
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))
Qs.interpolate(qs_func)
# 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.figure()
plt.imshow(np.flipud(Qs_interp.reshape((yaxis.size, xaxis.size))))
plt.colorbar()