Tabulating the bilinear form a

General
1 
from mpi4py import MPI
from dolfinx import mesh
import numpy
from dolfinx.fem import FunctionSpace 
from dolfinx import fem 
import ufl 
from petsc4py.PETSc import ScalarType 

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("CG", 1))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V) 
f = fem.Constant(domain, ScalarType(-6)) 
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx

From what I gather a is the bilinear form and v is a test function which at this point is unknown what it is.

So I would guess that a can be expressed with a symbolic expression. Is there a way to tabulate and print out what a is for each element u(x,y) what that a amounts to?

2

Yes, you can tabulate a bilinear form at any point x (in the reference element) for any cell with index i by using dolfinx.fem.Expression

https://docs.fenicsproject.org/dolfinx/v0.6.0/python/generated/dolfinx.fem.html#dolfinx.fem.Expression.eval

1 Like
3 
from mpi4py import MPI
from dolfinx import mesh
import numpy as np
from dolfinx.fem import FunctionSpace 
from dolfinx import fem 
import ufl 
from petsc4py.PETSc import ScalarType 


domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("CG", 1))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V) 
f = fem.Constant(domain, ScalarType(-6)) 
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx 
v = fem.Constant(domain.ufl_cell(), np.float32)

# from test file...
xtype = np.float32.real.dtype
domain2 = mesh.create_unit_square(MPI.COMM_WORLD, 4, 4, dtype=xtype)
v = fem.Constant(domain2.ufl_cell(), dtype(1))
u = fem.Function(FunctionSpace(domain2, ("Lagrange", 1)), dtype=dtype)
fem.Expression(v, u.function_space.element.interpolation_points())
fem.Expression(v, u.function_space.element.interpolation_points(), comm=MPI.COMM_WORLD)
fem.Expression(v, u.function_space.element.interpolation_points(), comm=MPI.COMM_SELF)

I brought this over from the test file. So far the difficulties I will have to work through are:

- The data typing (dtype) and how to free that from a pytest format.
- The ufl that I would like an expression to is ‘a’ not ‘v’, how to swap v, for a and get that to work.
- How to swap domain as opposed to domain2 which came from the pytest example.

Anyone can shed light on how to get those things working?

…
…
…

Just an aside but at the current time:

root@1bb2dcd10d7b:~/shared/dolfinx/python#  /usr/bin/env /bin/python3 /root/.vscode-server/extensions/ms-python.python-2023.14.0/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher 57547 -- /root/shared/dolfinx/python/test/unit/fem/test_expression.py 
Traceback (most recent call last):
  File "/root/shared/dolfinx/python/test/unit/fem/test_expression.py", line 12, in <module>
    from dolfinx.fem import (Constant, Expression, Function, FunctionSpace,
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/__init__.py", line 29, in <module>
    from dolfinx import common
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/common.py", line 11, in <module>
    from dolfinx import cpp as _cpp
ImportError: /usr/local/dolfinx-real/lib/libdolfinx.so.0.6: undefined symbol: _ZNK5basix13FiniteElement25interpolation_is_identityEv
root@1bb2dcd10d7b:~/shared/dolfinx/python#

when running the test file in question but at the current time I don’t know if I should try to compile cpp into the container or what might be wrong with symbols that are missing… I am using the r1 cpp container.

4

So far of test_expression.py it currently has a trouble importing basix.ufl with my dolfinx v0.6.0 installation from source on Ubuntu jammy. It appears that the test_expression_comm function is new to main brach of dolfinx. I don’t really know yet fully what:

@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128])

actually means yet. So far I would have to get the suggested code working on v0.6.0 since that I tried to get dolfinx main branch going but couldn’t yet get basix, ffcx and ufl to compile from source for a main branch source installation…

I guess to pose a question is the highlighted function feasible to implement in dolfinx v0.6.0 for python or do I have to recompile for main branch?

So far in v0.6.0 the fourth parameter to create_unit_square is for CellType. I feel a bit confused so far about what the data type is supposed to be defining in terms of the square and later on the call to Function…

I would so far expect type(0) to be a numpy float32 is this correct?

5

The only change that you should care about from main to v0.6.0 is that one can create a mesh with reduced precision. Removing the dtype argument will use float64.
You can simply look at the same file when released at v0.6.0, i.e.

which has similar examples (here tabulating the matrix at interpolation points, and generating a global interpolation matrix)

6

By modifying demo_poisson.py of the main brach dolfinx project at github"

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx + inner(g, v) * ds
# -
u = fem.Function(FunctionSpace(msh, ("Lagrange", 1)))
> fem.Expression(a, u.function_space.element.interpolation_points())
fem.Expression(a, u.function_space.element.interpolation_points(), comm=MPI.COMM_WORLD)
fem.Expression(a, u.function_space.element.interpolation_points(), comm=MPI.COMM_SELF)

So far the problem is:
Could not extract MPI communicator for Expression. Maybe you need to pass a communicator?

prusso@yoga:~/dolfinx$  cd /home/prusso/dolfinx ; /usr/bin/env /bin/python /home/prusso/.vscode/extensions/ms-python.python-2023.14.0/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher 48479 -- /home/prusso/dolfinx/python/demo/demo_poisson.py 
Could not extract MPI communicator for Expression. Maybe you need to pass a communicator?
Traceback (most recent call last):
  File "/home/prusso/dolfinx/python/demo/demo_poisson.py", line 134, in <module>
    fem.Expression(a, u.function_space.element.interpolation_points())
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/function.py", line 126, in __init__
    mesh = extract_unique_domain(e).ufl_cargo()
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/domain.py", line 345, in extract_unique_domain
    domains = extract_domains(expr)
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/domain.py", line 338, in extract_domains
    for t in traverse_unique_terminals(expr):
  File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/ufl/corealg/traversal.py", line 135, in traverse_unique_terminals
    if op._ufl_is_terminal_:
AttributeError: 'Form' object has no attribute '_ufl_is_terminal_'
prusso@yoga:~/dolfinx$

What does that line about extracting MPI communicator mean?

…
…
…

I copied petsc4py over from spack to site-distributions:

u = fem.Function(FunctionSpace(msh, ("Lagrange", 1)))
#fem.Expression(a, u.function_space.element.interpolation_points())
#fem.Expression(a, u.function_space.element.interpolation_points(), comm=MPI.COMM_WORLD)
#v = fem.Constant(msh.ufl_cell(),np.float64)
exp = fem.Expression(v, u.function_space.element.interpolation_points(), comm=MPI.COMM_SELF)

does give an expression exp. How is it possible to print out exp in tabulated form? Is this possible for a’s tuple?

7

By modifying your code to

import dolfinx
from mpi4py import MPI
import ufl
msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(msh, ("Lagrange", 1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds
expr = dolfinx.fem.Expression(a, V.element.interpolation_points(), comm=MPI.COMM_WORLD)

I get to

Traceback (most recent call last):
  File "/root/shared/mwe123.py", line 14, in <module>
    expr = dolfinx.fem.Expression(a, V.element.interpolation_points(), comm=MPI.COMM_WORLD)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 153, in __init__
    self._ufcx_expression, module, self._code = jit.ffcx_jit(comm, (e, _X),
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 55, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 211, in ffcx_jit
    raise TypeError(type(ufl_object))
TypeError: <class 'tuple'>

which is because expression actually cannot work with integral measures.

I would just like to ask for a sanity check. When you say tabulate the bilinear form, what do you imply?
Do you want the local element tensor for a single cell?
Or do you actually want to get the local element tensor evaluated at a single point in a single element (if so, what point, which cell, and what is the use-case?)

8

What my aim is to start tracing back out what is happening mathematically starting with the actual co-ordinates of a mesh. So I would like to start to achieve a spreadsheet like output of all the steps possible to look at. It’s so that I have a sheet to look at to better help me understand each step of the mathematical process.

Granted it seems like the bi-linear form is maybe midway through the process and maybe wasn’t the best place to start. I suppose that it would maybe be more appropriate to open a separate process so far I wouldn’t really know how best to title it…

9

I would consider: How to enumerate the cells of a mesh with dolfinx v0.6.0? - #2 by dokken that makes the local element matrix.

You could Also consider
https://simula-sscp.github.io/SSCP_2023_lectures/L12%20(FEniCS%20Intro)/L01_FEM_intro.html
That works through a 1D example

10

I found this useful which is what I was grasping for:

V = VectorFunctionSpace(msh, ("Lagrange", 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = form(inner(σ(u), grad(v)) * dx)  

vv = Function(FunctionSpace(msh, ("Lagrange", 1)))
exp = Expression(u, vv.function_space.element.interpolation_points(), comm=MPI.COMM_SELF)

L = form(inner(f, v) * dx)

#u.geometry.x

There it is possible to get a tabulation in matrix form of u, and also of v in some form to look at. Then it can start to be determined it how that u, and v start to go through the equations…

1 Like