I’m on Dolfinx 0.10.0 using the nightly build, and I’m currently tinkering around with previous Dolfinx tutorials to implement the usage of cross products to see how they work and how to use them for another program I’m working on. Currently I’m playing around with this older example for a mixed formulation of Poisson’s Equation which uses Dolfinx 0.7.1, and the code seems to work 1:1 between the two different Dolfinx versions. Upon retrieving the vector solution sigma_h by running
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions,
div, exp, inner, as_vector, cross, CellNormal, FacetNormal)
import ufl
import pyvista as pv
domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)
k = 1
Q_el = element("BDMCF", domain.basix_cell(), k)
P_el = element("DG", domain.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(domain, V_el)
Ve = fem.functionspace(domain, ("Lagrange", k + 1))
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
x = SpatialCoordinate(domain)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)
dx = Measure("dx", domain)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx
fdim = domain.topology.dim - 1
facets_top = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V.sub(0).collapse()
dofs_top = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_top)
def f1(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = np.sin(5 * x[0])
return values
f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V.sub(0))
facets_bottom = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V.sub(0), Q), fdim, facets_bottom)
def f2(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = -np.sin(5 * x[0])
return values
f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V.sub(0))
bcs = [bc_top, bc_bottom]
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
'petsc_options_prefix': 'mixed_poisson',
"pc_factor_mat_solver_type": "mumps"},
petsc_options_prefix="mixed_poisson_",
kind="mpi")
w_h = problem.solve()
sigma_h, u_h = w_h.split()
I wanted to try creating and plotting \boldsymbol \sigma^\perp = \hat{\textbf{k}} \times \boldsymbol \sigma. The result should be the same as \boldsymbol \sigma^\perp = (-\boldsymbol \sigma_y,\boldsymbol\sigma_x) which I ran below for reference
# Cross product
sigma_v = as_vector(sigma_h)
k = CellNormal(domain) # normal vector
sigma_cross = cross(k,sigma_v) # k x sigma
sigma_cross_true = as_vector((-sigma_h[1],sigma_h[0])) # True result of cross product
# Plotting
W = fem.functionspace(domain, ("DG", 0, (domain.geometry.dim, )))
sigma_v = ufl.as_vector(sigma_h)
sigma_d = fem.Function(W)
expr_d = fem.Expression(sigma_cross_true,W.element.interpolation_points)
sigma_d.interpolate(expr_d)
u_d = fem.Function(W)
expr_d = fem.Expression(sigma_h,W.element.interpolation_points)
u_d.interpolate(expr_d)
num_dofs = W.dofmap.index_map.size_local
num_verts = W.mesh.topology.index_map(0).size_local
values_d = np.zeros((num_dofs, 3), dtype=np.float64)
values_d[:, :domain.geometry.dim] = sigma_d.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs)
values = np.zeros((num_dofs, 3), dtype=np.float64)
values[:, :domain.geometry.dim] = u_d.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs)
grid = pv.UnstructuredGrid(*plot.vtk_mesh(Ve))
grid["sigma"] = values_d
glyphs = grid.glyph(orient = "sigma", scale = "sigma", factor = 0.1)# tolerance=0.1)
grid.set_active_vectors("sigma")
grid["u"] = values
grid.set_active_scalars("u")
u_plotter = pv.Plotter()
u_plotter.add_mesh(grid,scalars="u",show_edges= False)
u_plotter.add_mesh(glyphs, cmap = 'bwr')#, clim = [np.min(grid["sigma"]),np.max(grid["sigma"])])
u_plotter.view_xy()
u_plotter.show()
which yields the following plot
However, this code fails to compile when attempting to create a Dolfinx expression of the cross product, hence running the following
expr_d = fem.Expression(sigma_cross,W.element.interpolation_points)
returns an error
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[9], line 5
3 sigma_v = ufl.as_vector(sigma_h)
4 sigma_d = fem.Function(W)
----> 5 expr_d = fem.Expression(sigma_cross,W.element.interpolation_points)
6 sigma_d.interpolate(expr_d)
8 u_d = fem.Function(W)
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/function.py:151, in Expression.__init__(self, e, X, comm, form_compiler_options, jit_options, dtype)
149 form_compiler_options = dict()
150 form_compiler_options["scalar_type"] = dtype
--> 151 self._ufcx_expression, module, self._code = jit.ffcx_jit(
152 comm, (e, _X), form_compiler_options=form_compiler_options, jit_options=jit_options
153 )
154 self._ufl_expression = e
156 # Prepare coefficients data. For every coefficient in expression
157 # take its C++ object.
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:61, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
57 @functools.wraps(local_jit)
58 def mpi_jit(comm, *args, **kwargs):
59 # Just call JIT compiler when running in serial
60 if comm.size == 1:
---> 61 return local_jit(*args, **kwargs)
63 # Default status (0 == ok, 1 == fail)
64 status = 0
File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:218, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
216 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
217 elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
--> 218 r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
219 else:
220 raise TypeError(type(ufl_object))
File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:306, in compile_expressions(expressions, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
304 except Exception:
305 pass
--> 306 raise e
308 obj, module = _load_objects(cache_dir, module_name, expr_names)
309 return obj, module, (decl, impl)
File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:286, in compile_expressions(expressions, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
283 for name in expr_names:
284 decl += expression_template.format(name=name)
--> 286 impl = _compile_objects(
287 decl,
288 expressions,
289 expr_names,
290 module_name,
291 p,
292 cache_dir,
293 cffi_extra_compile_args,
294 cffi_verbose,
295 cffi_debug,
296 cffi_libraries,
297 visualise=visualise,
298 )
299 except Exception as e:
300 try:
301 # remove c file so that it will not timeout next time
File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:331, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
327 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
329 # JIT uses module_name as prefix, which is needed to make names of all struct/function
330 # unique across modules
--> 331 _, code_body = ffcx.compiler.compile_ufl_objects(
332 ufl_objects, prefix=module_name, options=options, visualise=visualise
333 )
335 # Raise error immediately prior to compilation if no support for C99
336 # _Complex. Doing this here allows FFCx to be used for complex codegen on
337 # Windows.
338 if sys.platform.startswith("win32"):
File /dolfinx-env/lib/python3.12/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
106 # Stage 1: analysis
107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"]) # type: ignore
109 _print_timing(1, time() - cpu_time)
111 # Stage 2: intermediate representation
File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:103, in analyze_ufl_objects(ufl_objects, scalar_type)
101 for original_expression, points in expressions:
102 elements += ufl.algorithms.extract_elements(original_expression)
--> 103 processed_expression = _analyze_expression(original_expression, scalar_type)
104 processed_expressions += [(processed_expression, points, original_expression)]
106 elements += ufl.algorithms.analysis.extract_sub_elements(elements)
File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:132, in _analyze_expression(expression, scalar_type)
130 """Analyzes and preprocesses expressions."""
131 preserve_geometry_types = (ufl.classes.Jacobian,)
--> 132 expression = ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(expression)
133 expression = ufl.algorithms.apply_derivatives.apply_derivatives(expression)
134 expression = ufl.algorithms.apply_function_pullbacks.apply_function_pullbacks(expression)
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_algebra_lowering.py:150, in apply_algebra_lowering(expr)
148 def apply_algebra_lowering(expr):
149 """Expands high level compound operators to equivalent representations using basic operators."""
--> 150 return map_integrand_dags(LowerCompoundAlgebra(), expr)
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py:86, in map_integrand_dags(function, form, only_integral_type, compress)
84 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
85 """Map integrand dags."""
---> 86 return map_integrands(
87 lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
88 )
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py:79, in map_integrands(function, form, only_integral_type)
77 elif isinstance(form, (Expr, BaseForm)):
78 integrand = form
---> 79 return function(integrand)
80 else:
81 raise ValueError("Expecting Form, Integral or Expr.")
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py:87, in map_integrand_dags.<locals>.<lambda>(expr)
84 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
85 """Map integrand dags."""
86 return map_integrands(
---> 87 lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
88 )
File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:35, in map_expr_dag(function, expression, compress, vcache, rcache)
15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
16 """Apply a function to each subexpression node in an expression DAG.
17
18 If the same function is called multiple times in a transformation
(...) 33 The result of the final function call
34 """
---> 35 (result,) = map_expr_dags(
36 function, [expression], compress=compress, vcache=vcache, rcache=rcache
37 )
38 return result
File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:103, in map_expr_dags(function, expressions, compress, vcache, rcache)
101 r = handlers[v._ufl_typecode_](v)
102 else:
--> 103 r = handlers[v._ufl_typecode_](v, *(vcache[u] for u in v.ufl_operands))
105 # Optionally check if r is in rcache, a memory optimization
106 # to be able to keep representation of result compact
107 if compress:
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_algebra_lowering.py:60, in LowerCompoundAlgebra.cross(self, o, a, b)
57 def c(i, j):
58 return Product(a[i], b[j]) - Product(a[j], b[i])
---> 60 return as_vector((c(1, 2), c(2, 0), c(0, 1)))
File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_algebra_lowering.py:58, in LowerCompoundAlgebra.cross.<locals>.c(i, j)
57 def c(i, j):
---> 58 return Product(a[i], b[j]) - Product(a[j], b[i])
File /dolfinx-env/lib/python3.12/site-packages/ufl/exproperators.py:383, in _getitem(self, component)
380 shape = self.ufl_shape
382 # Analyse slices (:) and Ellipsis (...)
--> 383 all_indices, slice_indices, repeated_indices = create_slice_indices(
384 component, shape, self.ufl_free_indices
385 )
387 # Check that we have the right number of indices for a tensor with
388 # this shape
389 if len(shape) != len(all_indices):
File /dolfinx-env/lib/python3.12/site-packages/ufl/index_combination_utils.py:150, in create_slice_indices(component, shape, fi)
148 elif isinstance(ind, int):
149 if int(ind) >= shape[len(all_indices)]:
--> 150 raise ValueError("Index out of bounds.")
151 all_indices.append(FixedIndex(ind))
152 elif isinstance(ind, slice):
ValueError: Index out of bounds.
Not sure what Index out of bounds is supposed to mean in this context. Also, I receive this error whether I use CellNormal or FacetNormal to create my unit vector.