Brief question on using dolfinx.fem.Expression to evaluate UFL expressions

dolfinx
1

I am trying to use Dolfinx Expression to evaluate a UFL expression at an arbitrary point inside of a mesh. I have had a look at various examples, but I cannot seem to get it to work. The present code worked for evaluating the solution directly, but the UFL expression does not work. What am I doing incorrectly?

import sys
import numpy as np
from scipy import optimize as opt
from scipy import special as sp
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import gmsh
import meshio
import ufl
import basix.ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx.fem.petsc import LinearProblem
from dolfinx import io
from dolfinx import geometry

eps = 1.0e-4 # Geometric tolerance

class Hinc:
    def __init__(self, kc, H):
        self.kc = kc 
        self.H = H
    def eval(self, x):
        hy = self.H * np.sin(self.kc * x[1])
        hx = np.full_like(hy, 0.0+0.0j, dtype=np.complex128)
        hz = np.full_like(hy, 0.0+0j, dtype=np.complex128)
        return(hx, hy, hz)

def Model(x):
    comm = nMPI.COMM_WORLD
    mpiRank = comm.rank
    modelRank = 0
    a = 2.286 # WR90 width (all dims in cm)
    b = 1.016 # WR90 height
    r = 1.115 # Circ WG radius
    l1 = 2.0 # WR90 length
    l2 = 2.0 # Circ WG length

    k = 2.0 * np.pi * x[0] / 3.0e10
    eta0 = 377.0
    kc = 1.84 / r
    kr = np.pi / a
    beta_c = np.sqrt(k * k - kc * kc) # circ WG imp
    beta_r = np.sqrt(k * k - kr * kr) # rect WG imp

    lm = 0.3
    ls = 0.3

    Plot = 1

    if mpiRank == modelRank:
        gmsh.initialize()
        gmsh.option.setNumber('General.Terminal', 1)
        gmsh.model.add('WG Transition')
        gmsh.model.occ.addBox(-b/2, 0, -r-l1, b, a, r+l1, 1)
        gmsh.model.occ.addCylinder(0, 0, 0, 0, l2+a, 0, r, 2, 2*np.pi)
        gmsh.model.occ.fuse([(3,1)],[(3,2)], 3, removeObject=True, removeTool=True)
        gmsh.model.occ.addBox(0, 0, -l1, r, l2+a, l1+2*r, 4)
        gmsh.model.occ.intersect([(3,3)], [(3,4)], 1, removeObject=True, removeTool=True)

        gmsh.model.occ.synchronize()

        pt = gmsh.model.getEntities(dim=0)
        gmsh.model.mesh.setSize(pt, lm)
        pt = gmsh.model.getEntitiesInBoundingBox(b/2-eps, -eps, -r-eps, b/2+eps, a+eps, eps)
        gmsh.model.mesh.setSize(pt, ls)
        pt
        gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
        gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
        gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
        gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
        gmsh.option.setNumber('Mesh.Format', 1)
        gmsh.option.setNumber('Mesh.MinimumCirclePoints', 30)
        gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
        gmsh.option.setNumber('Mesh.Tetrahedra', 0)

        bb = gmsh.model.getEntities(dim=3)
        pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
        WR90Port = []
        CircWGPort = []
        PEC = []
        PMC = []
        NoPMC = 1
        for bnd in pt:
            CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
            if np.abs(CoM[0]) < eps and NoPMC == 0:
                PMC.append(bnd[1])
            elif np.abs(CoM[2] + l1) < eps:
                WR90Port.append(bnd[1])
            elif np.abs(CoM[1] - l2 - a) < eps:
                CircWGPort.append(bnd[1])
            else:
                PEC.append(bnd[1])

        gmsh.model.addPhysicalGroup(3, [1], 1) # Wg interior
        gmsh.model.setPhysicalName(3, 1, "WG Interior")
        gmsh.model.addPhysicalGroup(2, PEC, 1)
        gmsh.model.setPhysicalName(2, 1, "PEC")
        gmsh.model.addPhysicalGroup(2, PMC, 2)
        gmsh.model.setPhysicalName(2, 2, "PMC")
        gmsh.model.addPhysicalGroup(2, WR90Port, 3)
        gmsh.model.setPhysicalName(2, 3, "InPort")
        gmsh.model.addPhysicalGroup(2, CircWGPort, 4)
        gmsh.model.setPhysicalName(2, 4, "Circ WG Port")

        gmsh.model.mesh.generate(3)
#        gmsh.fltk.run()


    mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
    if mpiRank == modelRank:
        gmsh.finalize()

    mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
    elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
    V = fem.functionspace(mesh, elem)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    n = ufl.FacetNormal(mesh)
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)

# Dirichlet BCs on PECs
    facets = fm.find(1)
    ubc = fem.Function(V)
    ubc.interpolate(lambda x:np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
    dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
    bc = fem.dirichletbc(ubc, dofs)
# Build RHS
    uh = fem.Function(V, dtype=np.complex128)
    f = Hinc(kr, 1.0+0j)
    uh.interpolate(f.eval)
    finc = 2.0j * k * eta0 * ufl.inner(ufl.cross(n, uh), v) * ds(3)
    uh.name = "Hinc"


# Weak form
    L = (ufl.inner(ufl.curl(u), ufl.curl(v)) - k * k * ufl.inner(u, v)) * ufl.dx
# Outgoing at inport
    bc1 = 1j * beta_r * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(3)
    L += bc1
# Circ WG  boundary
    bc2 = 1j * beta_c * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(4) 
    L += bc2

    num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
    num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
    print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
    if mpiRank == modelRank:
        print(f"Number of dofs global: {num_dofs_global}")

# Solve system of equations
    Lin_system = LinearProblem(L, finc, bcs=[bc], petsc_options={"ksp_type":"preonly", "pc_type":"lu"})
    E = Lin_system.solve()

# Try out the point evaluation
    bb_tree = geometry.bb_tree(mesh, mesh.topology.dim)
    points = np.array([[eps, l2+a-eps, 0.0], [r/2, l2+a-eps, 0.0]])
    potential_colliding_cells = geometry.compute_collisions_points(bb_tree, points)
    colliding_cells = geometry.compute_colliding_cells(mesh, potential_colliding_cells, points)

    points_on_proc = []
    cells = []
    for i, pp in enumerate(points):
        if len(colliding_cells.links(i)) > 0:
            points_on_proc.append(pp)
            cells.append(colliding_cells.links(i)[0])
    p_proc = np.array(points_on_proc, dtype=np.float64).reshape(-1, 3)
    St = fem.Expression(ufl.cross(ufl.conj(E), ufl.curl(E)), points)
    u_val = St.eval(p_proc, cells) / (1.0j * k * eta0)
    print(p_proc, u_val)
    
    return np.real(Pinc - Pref)

Model([1.0e10])

I get the following error.

Traceback (most recent call last):
  File "/home/bill/Cluster/Fenics2020/FeedHorn/WR90CircTransition.py", line 204, in <module>
Traceback (most recent call last):
  File "/home/bill/Cluster/Fenics2020/FeedHorn/WR90CircTransition.py", line 204, in <module>
    Model([1.0e10])
  File "/home/bill/Cluster/Fenics2020/FeedHorn/WR90CircTransition.py", line 200, in Model
    u_val = St.eval(p_proc, cells) / (1.0j * k * eta0)
  File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 222, in eval
    Model([1.0e10])
  File "/home/bill/Cluster/Fenics2020/FeedHorn/WR90CircTransition.py", line 200, in Model
    if (tdim := mesh.topology.dim) != (expr_dim := self._cpp_object.X().shape[1]):
AttributeError: 'numpy.ndarray' object has no attribute 'topology'
2

Update: I now understand that the problem lies in the first argument in the St.eval statement.

    St = fem.Expression(ufl.cross(ufl.conj(E), ufl.curl(E)), points)
    u_val = St.eval(p_proc, cells) / (1.0j * k * eta0)

When evaluating UFL expressions, replacing the “p_proc” array with the “mesh” designator, viz.

    St = fem.Expression(ufl.cross(ufl.conj(E), ufl.curl(E)), points)
    u_val = St.eval(mesh, cells) / (1.0j * k * eta0)

makes the code work. However, the points where the evaluations occur must be local reference element coordinates. I need to find a way to convert the global x,y,z coordinates to the (cell, reference_coordinate) pair. I already have the cell information, but I am at a loss as to how one goes about retrieving the local element coordinates from the global mesh coordinates.

3

I am getting closer, but I am running into a type mismatch error the mesh.geometry.cmap.pull_back(point, Vertices). The error messages are

(0, 181, array([1.0000e-04, 4.2859e+00, 0.0000e+00]))
[48 53 57 59]
[[ 0.00000000e+00  4.09367097e+00  1.32767104e-01]
 [ 2.40386685e-01  4.28600000e+00 -1.29733706e-01]
 [ 0.00000000e+00  4.02327166e+00 -1.33273093e-01]
 [ 0.00000000e+00  4.28600000e+00 -4.44089210e-16]]
Traceback (most recent call last):
  File "/home/bill/Cluster/Fenics2020/FeedHorn/WR90CircTransition.py", line 218, in <module>
    Model([1.0e10])
  File "/home/bill/Cluster/Fenics2020/FeedHorn/WR90CircTransition.py", line 208, in Model
    x_ref[i] = mesh.geometry.cmap.pull_back(point, Vertices[DoFs])
  File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/element.py", line 101, in pull_back
    return self._cpp_object.pull_back(x, cell_geometry)
TypeError: pull_back(): incompatible function arguments. The following argument types are supported:
    1. pull_back(self, x: ndarray[dtype=float64, shape=(*, *), order='C', writable=False], cell_geometry: ndarray[dtype=float64, shape=(*, *), order='C', writable=False]) -> numpy.ndarray[dtype=float64]

Invoked with types: dolfinx.cpp.fem.CoordinateElement_float64, ndarray, ndarray

The code snippet for the cell, point search is

# Try out the point evaluation
    bb_tree = geometry.bb_tree(mesh, mesh.topology.dim)
    points = np.array([[eps, l2+a-eps, 0], [r/3, l2+a-eps, 0], [2*r/3, l2+a-eps, 0.0]], dtype=mesh.geometry.x.dtype)
    potential_colliding_cells = geometry.compute_collisions_points(bb_tree, points)
    colliding_cells = geometry.compute_colliding_cells(mesh, potential_colliding_cells, points)

    points_on_proc = []
    cells = []
    for i, pp in enumerate(points):
        if len(colliding_cells.links(i)) > 0:
            points_on_proc.append(pp)
            cells.append(colliding_cells.links(i)[0])

    Vertices = mesh.geometry.x
    x_ref = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
    for i, (cell, point) in enumerate(zip(cells, points)):
        print((i, cell, point))
        DoFs = mesh.geometry.dofmap[cell]
        print(DoFs)
        print(Vertices[DoFs])
        x_ref[i] = mesh.geometry.cmap.pull_back(point, Vertices[DoFs])
    print(x_ref)

The error comes in the last line in the final for loop. Does anyone have an idea of where the mistyping error lies? I tried reshaping the arrays to no effect.

4

OK..I finally managed to get things to work. The final error had to do with incorrect shaping of the point array passed to the cmap.pull_back routine. Hopefully this is helpful. The final working code is

import sys
import numpy as np
from scipy import optimize as opt
from scipy import special as sp
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import gmsh
import meshio
import ufl
import basix.ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx.fem.petsc import LinearProblem
from dolfinx import io
from dolfinx import geometry

eps = 1.0e-4 # Geometric tolerance

class Hinc:
    def __init__(self, kc, H):
        self.kc = kc 
        self.H = H
    def eval(self, x):
        hy = self.H * np.sin(self.kc * x[1])
        hx = np.full_like(hy, 0.0+0.0j, dtype=np.complex128)
        hz = np.full_like(hy, 0.0+0j, dtype=np.complex128)
        return(hx, hy, hz)

def Model(pp_arg):
    comm = nMPI.COMM_WORLD
    mpiRank = comm.rank
    modelRank = 0
    a = 2.286 # WR90 width (all dims in cm)
    b = 1.016 # WR90 height
    r = 1.115 # Circ WG radius
    l1 = 2.0 # WR90 length
    l2 = 2.0 # Circ WG length

    k = 2.0 * np.pi * pp_arg[0] / 3.0e10
    eta0 = 377.0
    kc = 1.84 / r
    kr = np.pi / a
    beta_c = np.sqrt(k * k - kc * kc) # circ WG imp
    beta_r = np.sqrt(k * k - kr * kr) # rect WG imp

    lm = 0.3
    ls = 0.15

    Plot = 1

    if mpiRank == modelRank:
        gmsh.initialize()
        gmsh.option.setNumber('General.Terminal', 1)
        gmsh.model.add('WG Transition')
        gmsh.model.occ.addBox(-b/2, 0, -r-l1, b, a, r+l1, 1)
        gmsh.model.occ.addCylinder(0, 0, 0, 0, l2+a, 0, r, 2, 2*np.pi)
        gmsh.model.occ.fuse([(3,1)],[(3,2)], 3, removeObject=True, removeTool=True)
        gmsh.model.occ.addBox(0, 0, -l1, r, l2+a, l1+2*r, 4)
        gmsh.model.occ.intersect([(3,3)], [(3,4)], 1, removeObject=True, removeTool=True)

        gmsh.model.occ.synchronize()

        pt = gmsh.model.getEntities(dim=0)
        gmsh.model.mesh.setSize(pt, lm)
        pt = gmsh.model.getEntitiesInBoundingBox(b/2-eps, -eps, -r-eps, b/2+eps, a+eps, eps)
        gmsh.model.mesh.setSize(pt, ls)
        pt
        gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
        gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
        gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
        gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
        gmsh.option.setNumber('Mesh.Format', 1)
        gmsh.option.setNumber('Mesh.MinimumCirclePoints', 30)
        gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
        gmsh.option.setNumber('Mesh.Tetrahedra', 0)

        bb = gmsh.model.getEntities(dim=3)
        pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
        WR90Port = []
        CircWGPort = []
        PEC = []
        PMC = []
        NoPMC = 1
        for bnd in pt:
            CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
            if np.abs(CoM[0]) < eps and NoPMC == 0:
                PMC.append(bnd[1])
            elif np.abs(CoM[2] + l1) < eps:
                WR90Port.append(bnd[1])
            elif np.abs(CoM[1] - l2 - a) < eps:
                CircWGPort.append(bnd[1])
            else:
                PEC.append(bnd[1])

        gmsh.model.addPhysicalGroup(3, [1], 1) # Wg interior
        gmsh.model.setPhysicalName(3, 1, "WG Interior")
        gmsh.model.addPhysicalGroup(2, PEC, 1)
        gmsh.model.setPhysicalName(2, 1, "PEC")
        gmsh.model.addPhysicalGroup(2, PMC, 2)
        gmsh.model.setPhysicalName(2, 2, "PMC")
        gmsh.model.addPhysicalGroup(2, WR90Port, 3)
        gmsh.model.setPhysicalName(2, 3, "InPort")
        gmsh.model.addPhysicalGroup(2, CircWGPort, 4)
        gmsh.model.setPhysicalName(2, 4, "Circ WG Port")

        gmsh.model.mesh.generate(3)
#        gmsh.fltk.run()


    mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
    if mpiRank == modelRank:
        gmsh.finalize()

    mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
    elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
    V = fem.functionspace(mesh, elem)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    n = ufl.FacetNormal(mesh)
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)

# Dirichlet BCs on PECs
    facets = fm.find(1)
    ubc = fem.Function(V)
    ubc.interpolate(lambda x:np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
    dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
    bc = fem.dirichletbc(ubc, dofs)
# Build RHS
    uh = fem.Function(V, dtype=np.complex128)
    f = Hinc(kr, 1.0+0j)
    uh.interpolate(f.eval)
    finc = 2.0j * k * eta0 * ufl.inner(ufl.cross(n, uh), v) * ds(3)
    uh.name = "Hinc"


# Weak form
    L = (ufl.inner(ufl.curl(u), ufl.curl(v)) - k * k * ufl.inner(u, v)) * ufl.dx
# Outgoing at inport
    bc1 = 1j * beta_r * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(3)
    L += bc1
# Circ WG  boundary
    bc2 = 1j * beta_c * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(4) 
    L += bc2

    num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
    num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
    print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
    if mpiRank == modelRank:
        print(f"Number of dofs global: {num_dofs_global}")

# Solve system of equations
    Lin_system = LinearProblem(L, finc, bcs=[bc], petsc_options={"ksp_type":"preonly", "pc_type":"lu"})
    E = Lin_system.solve()

# Try out the point evaluation
    bb_tree = geometry.bb_tree(mesh, mesh.topology.dim)
    points = np.array([[eps, l2+a-eps, 0], [r/3, l2+a-eps, 0], [2*r/3, l2+a-eps, 0.0]], dtype=mesh.geometry.x.dtype)
    potential_colliding_cells = geometry.compute_collisions_points(bb_tree, points)
    colliding_cells = geometry.compute_colliding_cells(mesh, potential_colliding_cells, points)
    points_on_proc = []
    cells = []
    for i, pp in enumerate(points):
        if len(colliding_cells.links(i)) > 0:
            points_on_proc.append(pp)
            cells.append(colliding_cells.links(i)[0])
    
    Vertices = mesh.geometry.x
    x_ref = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
    for i, (cell, point) in enumerate(zip(cells, points)):
        DoFs = mesh.geometry.dofmap[cell]
        print("Rank = {0}".format(mpiRank))
        x_ref[i] = mesh.geometry.cmap.pull_back(point.reshape(-1,3), Vertices[DoFs])
        St = fem.Expression(ufl.cross(ufl.conj(E), ufl.curl(E)), x_ref[i], comm=nMPI.COMM_SELF)
        u_val = St.eval(mesh, np.asarray([cell], dtype=np.int32)) / (1.0j * k * eta0)
        print(u_val)

    
    
    return np.real(Pinc - Pref)
Model([1.05e10])
1 Like