Nodal coordinates and values with MPI

dolfinx
, ,
1

Hello,

I recently installed fenicsx with dolfinx 0.9.0.0. I’m working on a 3D problem, where I read a gmsh file with the gmshio.read_from_msh command. The definition of the variational problem is working and the visualization in paraview works as tought. To speed up the simulation I do some parallelization, which again works out of the box. Now I want to add a step and evaluate the nodal values and coordinates of a surface, which i defined via PhysicalGroup in gmsh.

The following code snipped runs perfectly on a single core:

facets = facet_tags.find(surf_tag)

indices = mesh.entities_to_geometry(domain, domain.topology.dim-1, facets, False)
indices = np.unique(indices)

coordinates = domain.geometry.x[indices]
vertex_realvalues = E.x.array[indices].real
vertex_imagvalues = E.x.array[indices].imag

outpfile = open('results/ScatterE_1CPU.dat','w')
outpfile.write('# %11s %22s %22s %22s %22s\n'%('x','y','z','E_real','E_imag'))
for coord, Erealval, Eimagval in zip(coordinates, vertex_realvalues, vertex_imagvalues):
    outpfile.write('%11.15f %22.15f %22.15f %22.15f %22.15f\n'%(coord[0], coord[1], coord[2], Erealval, Eimagval))
outpfile.close()

However, when parallelizing this, I’m facing errors with the message, division by zero. I’ve read, that its related to a bug in mesh.entities_to_geometry which is fixed in dolfinx 0.10.0.0, which is not yet in the apt repositories to install. Is there a workaround for now? I figured, that due to the parallelization, the facets variable may be empty, hence I followed with a allCPUfacets = domain.comm.gather(facets,root=0) command, which is later reduced to a list again via

facetslistoflist = [list(nodeid) for nodeid in allCPUfacets if list(nodeid)]
facetslist = [nodeid for nodeidlist in facetslistoflist for nodeid in nodeidlist]
facets = np.asarray(facetslist)

on MPI.COMM_WORLD.Get_rank() == 0 . This procedure is repeated with the indices variable and so forth, but did not reach a satisfactory result. With the mesh.compute_incident_entities command the division by zero error is not exiting the program but the indices then point to different nodes, since the evaluated coordinates are not on a plane anymore.

Can anyone point me towards the correct evaluation of coordinates and values of nodes for parallelized objects? Maybe a .gather or .receive combined with the evaluation on a single core is all I need since it’s just a post-processing step. But I’m fairly new to fenics and mpi4py and I seem to struggle to adapt from old dolfin, when I find some clues in the forum.

Thanks in advance.

2

A workaround would be to adapt this to:

if len(facets) == 0:
    indices = np.empty([], dtype=np.int32)
else: 
    indices = mesh.entities_to_geometry(domain, domain.topology.dim-1, facets, False)
indices = np.unique(indices)

However, the code:

is not really safe in its current status, as you have some heavy implicit assumptions on the mesh geometry and the function space.

I am assuming you are using a first order Lagrange space for your computations?

If so, one should rather use a vertex_to_dofmap construction, as for instance shown in Extract the coordinates of the boundary points and the value of a function on them - #4 by dokken

3

Thank you for the workaround and the hint for the evaluation. I’m working on latter and get another division by zero error but I’m still figuring out its origin. I will post, when I fixed it or with another question.

I’m working with Nedelec elements “N1curl” but the results are projected onto first order Lagrange space for visualization in paraview as in this PML demo. Hence, I planned to evaluate the quantity here as well. What would differ with other kinds of elements? By evaluating the dofs as you suggested, I would get the corresponding (non-interpolated) values and coordinates from my understanding.

4

Moving N1curl into a first order continuous Lagrange space is not a great idea. The space should be discontinuous, but then you end up with multiple degrees of freedom per vertex, and an averaging strategy (or similar) should be consider.

N1curl spaces additionally have no degrees of freedom at the vertices, but rather defined through integral moments of the faces: DefElement

5

The problem is solved in Nedelec space and I interpolate it by a corresponding DG space via

V_dg = fem.functionspace(domain, ("DG", degree, (3,)))
E_dg = fem.Function(V_dg, name="E")
E_dg.interpolate(E)

, is this the correct approach? My plan is to evaluate the corresponding nodal coordinates and values for further processing. Is it easier or the better approach to evaluate the integral moments? I thought identifying the faces and the corresponding surface and moment for deriving the quantity is more complicated.

6

Is there a specific reason you want the values at the vertices? Would midpoints be sufficient for your case, as those values are uniquely defined on exterior facets (as opposed to values at vertices).

7

midpoint is sufficient as well, as long as i can have the coordinates.

8

Then you can use a code like:

from mpi4py import MPI
import dolfinx
import ufl
import basix
import numpy as np
import numpy.typing as npt


def evaluate_expression_at_facet_midpoints(
    expr: ufl.core.expr.Expr, mesh: dolfinx.mesh.Mesh, facets: npt.NDArray[np.int32]
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.inexact]]:
    # Get coordinate element for facets of cell
    facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
    unique_facet_types = np.unique(facet_types)
    assert len(unique_facet_types) == 1, (
        "All facets must have the same type for interpolation."
    )
    facet_type = facet_types[0]
    facet_geom = basix.cell.geometry(facet_type)
    facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]

    # Compile code for evaluation
    bndry_expr = dolfinx.fem.Expression(expr, facet_midpoint)

    # Map facet indices to (cell, local_facet) pairs
    try:
        boundary_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets
        )
    except TypeError:
        boundary_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet,
            mesh.topology,
            facets,
            mesh.topology.dim - 1,
        )

    x = ufl.SpatialCoordinate(mesh)
    coord_expr = dolfinx.fem.Expression(x, facet_midpoint)

    return coord_expr.eval(mesh, boundary_entities), bndry_expr.eval(
        mesh, boundary_entities
    )


mesh = dolfinx.mesh.create_unit_cube(
    MPI.COMM_WORLD, 10, 10, 10, cell_type=dolfinx.mesh.CellType.hexahedron
)
V = dolfinx.fem.functionspace(mesh, ("N1curl", 4))
u = dolfinx.fem.Function(V)
u.interpolate(
    lambda x: (x[1] * x[0] ** 2 + x[1] ** 2, np.zeros(x[0].shape), np.zeros(x[0].shape))
)

n = ufl.FacetNormal(mesh)
fdim = mesh.topology.dim - 1
facet_for_evaluation = dolfinx.mesh.locate_entities_boundary(
    mesh, fdim, lambda x: np.isclose(x[0], 1)
)
expr = ufl.dot(ufl.grad(u), n)

coordinates, values = evaluate_expression_at_facet_midpoints(
    expr, mesh, facet_for_evaluation
)
print(coordinates)
print(values)

ref_values = (2 * coordinates.T[0] * coordinates.T[1]).T
np.testing.assert_allclose(values.T[0], ref_values)
np.testing.assert_allclose(values.T[1], 0, atol=1e-12, rtol=1e-12)
np.testing.assert_allclose(values.T[2], 0, atol=1e-12, rtol=1e-12)
9

This works great, thank you very much. Now I gotta adapt it to my problem. Is there an issue with overlap due to the parallelization? Interestingly, the coordinates are complex valued as well. however, the imaginary part is zero.

10

You can always do a filter on “owned” facets after calling locate_entities_boundary, which would remove any potential overlap from ghost facets.

That is because I used dolfinx.fem.Expression to get the facet midpoints, instead of calling

dolfinx.fem.compute_midpoints (dolfinx.mesh — DOLFINx 0.10.0.0 documentation) which you could use if you would like to.

11

Thank you very much for the clarification! It’s a great help and boosted my understanding of fenics.

12

Is there a possibility to evaluate a planar surface somewhere within the domain? For a transition between media, i.e. two cuboids next to each other, one may evaluate the common plane as well. Hence, the Integral in the evaluate_expression_at_facet_midpoints function is not of exterior_facet type anymore if i got this correctly. The simple change to interior_facets seems insufficient. The obvious change to facet_for_evaluation = dolfinx.mesh.locate_entities(…) instead of dolfinx.mesh.locate_entities_boundary(…) provided.

The corresponding mesh:

import gmsh
import sys
gmsh.initialize(sys.argv)
dim=3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(dim, [box], tag=801, name='box')
gmsh.model.addPhysicalGroup(dim, [box2], tag=802, name='box2')

surfaces = gmsh.model.getEntities(dim=dim-1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:# to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name='evalsurf')
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim-1, othersurf, tag=902, name='other')

msize = 0.1
points = gmsh.model.getEntities(dim=dim-3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim-1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write('doublebox.msh')


mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("doublebox.msh", MPI.COMM_WORLD, 0, gdim=3)

example and read of said mesh to include in the above code.

EDIT: Nevermind, it works with exterior_facet and even with facet_for_evaluation = facet_tags.find(901) .

13

If you want this to work on interior facets, then the logic here should be revised (as one needs to pick a given side of the facet to evaluate at).

14

In this answer of you, you made an example where you differentiated the sides. It is all based on a first construction of the entities of the interior_facet integration domain. However, this array is empty for the given gmsh example. Maybe its due to the construction with two boxes, so it is somehow regarded as exterior facet, even though this would still contradict the logic since the side is not neccessarily picked correctly.

Here is the full code, where the above function is copied in the code for debugging:

from mpi4py import MPI
import dolfinx
import ufl
import basix
import numpy as np
import numpy.typing as npt
import gmsh
import sys


gmsh.initialize(sys.argv)
dim=3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(dim, [box], tag=801, name='box')
gmsh.model.addPhysicalGroup(dim, [box2], tag=802, name='box2')

surfaces = gmsh.model.getEntities(dim=dim-1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:# to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name='evalsurf')
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim-1, othersurf, tag=902, name='other')

msize = 0.1
points = gmsh.model.getEntities(dim=dim-3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim-1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write('doublebox.msh')


mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("doublebox.msh", MPI.COMM_WORLD, 0, gdim=3)

V = dolfinx.fem.functionspace(mesh, ("N1curl", 3))
u = dolfinx.fem.Function(V)
u.interpolate(
    lambda x: (x[1] * x[0] ** 2 + x[1] ** 2, np.zeros(x[0].shape), np.zeros(x[0].shape))
)

n = ufl.FacetNormal(mesh)
fdim = mesh.topology.dim - 1

facet_for_evaluation = facet_tags.find(901)

expr = ufl.dot(ufl.grad(u), n)

# for renaming due to copying the function
facets = facet_for_evaluation
meshdom = mesh
# function
facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
unique_facet_types = np.unique(facet_types)
assert len(unique_facet_types) == 1, (
    "All facets must have the same type for interpolation."
)
facet_type = facet_types[0]
facet_geom = basix.cell.geometry(facet_type)
facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]

# Compile code for evaluation
bndry_expr = dolfinx.fem.Expression(expr, facet_midpoint)

# Map facet indices to (cell, local_facet) pairs
try:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets
    )
except TypeError:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.exterior_facet,
        mesh.topology,
        facets,
        mesh.topology.dim - 1,
    )

x = ufl.SpatialCoordinate(mesh)
coord_expr = dolfinx.fem.Expression(x, facet_midpoint)
# return of function
coordinates = coord_expr.eval(mesh, boundary_entities)
values  =bndry_expr.eval(mesh, boundary_entities)

# handling of interior facets
interior_indices = dolfinx.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.interior_facet, mesh.topology, facets,
    mesh.topology.dim - 1)

#print(coordinates)
#print(values)

The interior_indices variable gives an empty array. Hence, the picking of the side is impossible. The handling as exterior facets even though the facet with tag 901 is the interior connection of the two created cubes.