I inherited a python implementation of a domain decomposition method for a 2d Poisson problem for a colleague which has been working fine with an older version (0.7.1) of fenics. Now I want to port the problem to the current version and am pretty far along. There is however one problem that I cannot quite figure out:
The formulation contains several jump terms of the form u('+') * n('+') + u('-') * n('-') which should be integrated over a boundary interface with n = ufl.FacetNormal(mesh) based on a function space of DG elements. To this end, a measure dS is created based on a MeshTags object containing the interface. However, when I try to create the form using dolfinx.fem.form(A) I get a
*** ValueError: Discontinuous type CellFacetJacobian must be restricted.
Could you explain to me what this error means (I am not that familiar with fenics in particular)? As far as I can see the term restriction means that a term is restricted to one side of a facet, which it should be, right? Did I not select the correct DOFs?
Without the actual form it is hard to tell you what you should restrict.
In general, the idea is that the dS measure integrates over an interior facet (a facet connected to two cells). Which cell you would like to extract the values from, depends on the restriction, (“+” and “-”) which is arbitrarily assigned.
For most classical finite element schemes, it doesn’t matter which side is first in a jump, as long as it is consistent through the variational form.
Thank you for the information so far. My problem is that I cannot really share the entire code and unfortunately I don’t have a MWE.
I do understand how jump terms work in principle and I don’t think that there is a problem with the signs at all. I just don’t get why the error occurs and what it means. The message reads like I should add a restriction but that seems backwards to me, after all the restrictions +/- seem to be the cause of the problem.
It means that some quantity in the form is not restricted as it should. To debug this you should take your form A, and remove one and one term until the error message from dolfinx.fem.form(A) disappears. This will point you to what terms that are problematic and make it easier to make a reproducible example.
@dokken: Here is a (hopefully) fairly minimal example. Note that this is not a weak formulation of an PDE, but an example which triggers the behavior I observe (present in fenics-dolfinx 0.1.0 from conda).
import numpy as np
import dolfinx
from mpi4py import MPI
import gmsh
from dolfinx.io import gmsh as gmshio
from basix.ufl import element as FiniteElement
import ufl
from dolfinx import fem
from ufl import inner, dS
gmsh.initialize()
p1 = gmsh.model.occ.addPoint(0, 0, 0)
p2 = gmsh.model.occ.addPoint(1, 0, 0)
p3 = gmsh.model.occ.addPoint(1, 1, 0)
p4 = gmsh.model.occ.addPoint(0, 1, 0)
line1 = gmsh.model.occ.add_line(p1, p2)
line2 = gmsh.model.occ.add_line(p2, p3)
line3 = gmsh.model.occ.add_line(p3, p4)
line4 = gmsh.model.occ.add_line(p4, p1)
big_rect_loop = gmsh.model.occ.addCurveLoop([line1, line2, line3, line4])
big_rect = gmsh.model.occ.addPlaneSurface([big_rect_loop])
tools = [(1, line1), (1, line2), (1, line3), (1, line4)]
tools.append((2, gmsh.model.occ.add_rectangle(
0, 0, 0, 1 - 0, 1 - 0)))
gmsh.model.occ.synchronize()
result = gmsh.model.occ.fragment([(2, big_rect)], tools)
gmsh.model.occ.synchronize()
index = 1
boundary_lines = []
for i in range(index, index + 4):
boundary_lines += np.array(result[1][i])[:, 1].tolist()
gmsh.model.addPhysicalGroup(1, boundary_lines, -1000)
index += 4
gmsh.model.addPhysicalGroup(2, [1], 0)
gmsh.model.occ.synchronize()
h = 0.04
gmsh.model.mesh.setSize(gmsh.model.occ.getEntities(dim=0), h)
gdim = 2
gmsh.model.mesh.generate(gdim)
msh_data = gmshio.model_to_mesh(model=gmsh.model,
comm=MPI.COMM_WORLD,
rank=0,
gdim=gdim)
msh = msh_data.mesh
cell_markers = msh_data.cell_tags
facet_markers = msh_data.facet_tags
msh.topology.create_entities(1)
msh.topology.create_connectivity(0, 1)
gmsh.finalize()
P0 = FiniteElement("DG", str(msh.ufl_cell()), 0)
V = dolfinx.fem.functionspace(msh, P0)
u, = ufl.TrialFunctions(V)
v, = ufl.TestFunctions(V)
n = ufl.FacetNormal(msh)
A = inner(u, v) * dS
A_mat = fem.assemble_matrix(fem.form(A)).to_scipy()
Specifically, I get a ValueError: Discontinuous type Argument must be restricted. This has apparently something to do with the measure being used, it works fine with ds
In the actual example the mesh is decomposed into subdomains and I want to integrate over the boundary between two adjacent subdomains.
The error here is quite clear.
The trial functions u and v are discontinuous functions (piecewise constant).
Thus, when integrating over an interior facet, you have to decide on which value you should access for u and v.
In UFL this is described by the "+" and "-" restriction.
However, most people use ufl.jump(u), ufl.avg(u) to get the correct variational forms.
As you haven’t really specified what should happen at the interface, it is not easy for me to say if you want to use u("+") or u("-").
Thank you for your help so far, I think I understand better what the problem is. I narrowed the problem down to something I thought was insignificant before, namely the facet area. Specifically, the form I am looking at is properly restricted expcept for the facet area, which is causing problems. If in the example you replace the form by A = FacetArea(msh) * dS, the error
ValueError: Discontinuous type CellFacetJacobian must be restricted.
appears. I thought that this was a trivial term, but apparently this area value must be restricted as well. I don’t quite understand why though. Shouldn’t the area be just a constant for each facet?
@memoized_handler
def facet_jacobian(self, o):
"""Apply to facet_jacobian."""
if self._preserve_types[o._ufl_typecode_]:
return o
domain = extract_unique_domain(o)
J = self.jacobian(Jacobian(domain))
RFJ = CellFacetJacobian(domain)
i, j, k = indices(3)
return as_tensor(J[i, k] * RFJ[k, j], (i, j))
where we now see the general jacobian of the mesh (J), which is dependent on which side of the cell we are.
Thus from an algorithmic perspective, one should get the same result from either side of the facet, but it will be computed with different values on either side.
Therefore, I suggest you restrict it to the “+” side.
Thank you for the information, the creation of the form seems to work correctly. However, assembling the matrix seems to cause an error on the C++ side of things (based on conda version 0.10.0):
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ffcx/ir/representation.py", line 402, in _compute_integral_ir
integral_ir = compute_integral_ir(
itg_data.domain.ufl_cell(),
...<5 lines>...
visualise,
)
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ffcx/ir/integral.py", line 144, in compute_integral_ir
expression = balance_modifiers(expression)
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ufl/algorithms/balancing.py", line 97, in balance_modifiers
return map_expr_dag(mf, expr)
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
(result,) = map_expr_dags(
~~~~~~~~~~~~~^
function, [expression], compress=compress, vcache=vcache, rcache=rcache
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ufl/corealg/map_dag.py", line 109, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *(vcache[u] for u in v.ufl_operands))
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ufl/algorithms/balancing.py", line 83, in _modifier
return balance_modified_terminal(expr)
File "/.conda/envs/coupled_quadratic/lib/python3.14/site-packages/ufl/algorithms/balancing.py", line 54, in balance_modified_terminal
assert expr._ufl_is_terminal_modifier_
It looks like the FacetArea must be restricted in the last two (also restricted) parts of the form but not in the first part. At least everything seems to work with