Numerical values from ufl.SpatialCoordinate

dolfinx
2

ufl.SpatialCoordinate does not imply the “mesh-nodes” (meaning the vertices for linear meshes). The x=ufl.SpatialCoordinate(mesh) represents any quadrature point used during integration in its physical space.

What I would suggest to do, is to interpolate your function (non-analytical) into a suitable function space. Say, if you want to actually use the mesh nodes (vertices), you can create a first order Lagrange function and interpolate into that, i.e.

V  = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(V)

def your_expression(x):
     # map your values here
     return value
u.interpolate(your_expression)

If you instead you want to use the quadrature points (as described by x=ufl.SpatialCoordinate(mesh), I suggest something like the following:

from IPython import embed
import dolfinx
from mpi4py import MPI
import ufl
import basix
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)


def height(x):
    # Add whatever function you want to call here, for now choose identity
    return x[0]


quadrature_degree = 4
dx = ufl.Measure("dx", domain=mesh, metadata={
                 "quadrature_degree": quadrature_degree})
quadrature_points, wts = basix.make_quadrature(
    basix.cell.string_to_type(mesh.ufl_cell().cellname()), quadrature_degree)
Q_element = ufl.FiniteElement(
    "Quadrature", mesh.ufl_cell(), quadrature_degree, quad_scheme="default")
Q = dolfinx.fem.FunctionSpace(mesh, Q_element)
height_function = dolfinx.fem.Function(Q)

x = ufl.SpatialCoordinate(mesh)
x_expr = dolfinx.fem.Expression(x, quadrature_points)
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
integration_points = x_expr.eval(np.arange(num_cells, dtype=np.int32))

num_points_per_cell = quadrature_points.shape[0]
gdim = mesh.geometry.dim
assert(integration_points.shape[1] == num_points_per_cell*gdim)
assert(integration_points.shape[0] == num_cells)
for i in range(num_cells):
    for j in range(num_points_per_cell):
        coord = integration_points[i, gdim*j:gdim*(j+1)]
        value = height(coord)
        height_function.x.array[i*num_points_per_cell +
                                j:i*num_points_per_cell + (j+1)] = value

for i in range(mesh.geometry.dim):
    # verify that height_function maps to x[0]
    print(dolfinx.fem.assemble_scalar(
        dolfinx.fem.form((height_function-x[0])**2*dx)))


# Use heigh function in general()
V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = dolfinx.fem.form(height_function*ufl.inner(ufl.grad(u), ufl.grad(v))*dx)
A = dolfinx.fem.petsc.assemble_matrix(a)
1 Like