Integrate a flux over an internal curve

raboynics · May 25, 2025, 8:09am

I have a 2D geometry with an annulus shape and 2 internal curves, one at \theta=0, the other at \theta=\pi/2 (polar coordinates). I solve 2 coupled PDEs for the variables T and V (approximated as continuous Lagrange polynomials of degree 4). Then, after solving the FEM problem, I form \vec J in terms of T and V as \vec J = -\sigma \nabla V -\sigma S \nabla T where \sigma is a scalar that depends on the region:

R = functionspace(the_mesh.mesh, ("DG", 0))
sigma = Function(R)
sigma.x.array[three_quarters_cells] = np.full_like(three_quarters_cells, 1/(1e0), dtype=default_scalar_type)
sigma.x.array[quarter_cells] = np.full_like(quarter_cells, 1/1, dtype=default_scalar_type)

And S is a tensor (a diagonal 2x2 matrix) that also depends on the region:


        s_xx, s_yy = self.input_parameters.seebeck_components

        # Dokken's solution below
        tensor_el = element(family="DG", cell=str(the_mesh.mesh.ufl_cell()), degree=0, shape=(2, 2))

        T = functionspace(the_mesh.mesh, tensor_el)
        Seebeck_tensor = Function(T)

        def seebeck_three_quarters(x):
            tensor = np.array([[0, 0], [0, 0]])
            values = np.repeat(tensor, x.shape[1])
            #print(values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1]))
            return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])


        def seebeck_quarter(x):
            tensor = np.array([[s_xx, 0], [0, s_yy]])
            values = np.repeat(tensor, x.shape[1])
            #print(values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1]))
            return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])

        Seebeck_tensor.interpolate(seebeck_three_quarters, cells0=self.the_mesh.cell_markers.find(20))
        Seebeck_tensor.interpolate(seebeck_quarter, cells0=self.the_mesh.cell_markers.find(21))
        self.Seebeck_tensor = Seebeck_tensor

The usual


I_theta0 = assemble_scalar(form(dot(J_vector, n) * dS(22)))

print(f'I_theta0: {I_theta0}')

yields an error about restrictions whose traceback finishes by


  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/apply_restrictions.py", line 61, in _require_restriction

    raise ValueError(f"Discontinuous type {o._ufl_class_.__name__} must be restricted.")

ValueError: Discontinuous type Coefficient must be restricted.

I tried several things, like using (‘+’), but nothing worked.

For example

    J_integrand_on_facet = (-sigma('+') * grad(volt)('+') 
                                - sigma('+') * self.thermoelec_sim.Seebeck_tensor('+') * grad(temp)('+'))
        
        # Or, relying on UFL to correctly handle traces of gradients of CG functions:
        # J_integrand_on_facet = (-sigma('+') * grad(volt) 
        #                         - sigma('+') * self.thermoelec_sim.Seebeck_tensor('+') * grad(temp))

        I_theta0 = assemble_scalar(form(dot(J_integrand_on_facet, n) * dS(22)))
        print(f'I_theta0: {I_theta0}')

yields

Traceback (most recent call last):

  File "/home/isaac/Documents/fenicsx/Bridgman/2d_bridgman_textbook_example/debug_quarter_annulus_renormalize_PF/full_annulus_Lukosz_v0dot9.py", line 338, in <module>

    post_proc.post_proc()

    ~~~~~~~~~~~~~~~~~~~^^

  File "/home/isaac/Documents/fenicsx/Bridgman/2d_bridgman_textbook_example/debug_quarter_annulus_renormalize_PF/full_annulus_Lukosz_v0dot9.py", line 322, in post_proc

    I_theta0 = assemble_scalar(form(dot(J_integrand_on_facet, n) * dS(22)))

                               ~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py", line 337, in form

    return _create_form(form)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py", line 331, in _create_form

    return _form(form)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py", line 254, in _form

    ufcx_form, module, code = jit.ffcx_jit(

                              ~~~~~~~~~~~~^

        mesh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options

        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    )

    ^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/jit.py", line 62, in mpi_jit

    return local_jit(*args, **kwargs)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/jit.py", line 212, in ffcx_jit

    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py", line 225, in compile_forms

    raise e

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py", line 205, in compile_forms

    impl = _compile_objects(

        decl,

    ...<9 lines>...

        visualise=visualise,

    )

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py", line 330, in _compile_objects

    _, code_body = ffcx.compiler.compile_ufl_objects(

                   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^

        ufl_objects, prefix=module_name, options=options, visualise=visualise

        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    )

    ^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/compiler.py", line 108, in compile_ufl_objects

    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects

    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/analysis.py", line 94, in <genexpr>

    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)

                      ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/analysis.py", line 180, in _analyze_form

    form_data = ufl.algorithms.compute_form_data(

        form,

    ...<6 lines>...

        complex_mode=complex_mode,

    )

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/compute_form_data.py", line 337, in compute_form_data

    form = apply_restrictions(form)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/apply_restrictions.py", line 183, in apply_restrictions

    return map_integrand_dags(rules, expression, only_integral_type=integral_types)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/map_integrands.py", line 86, in map_integrand_dags

    return map_integrands(

        lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type

    )

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/map_integrands.py", line 30, in map_integrands

    map_integrands(function, itg, only_integral_type) for itg in form.integrals()

    ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands

    return itg.reconstruct(function(itg.integrand()))

                           ~~~~~~~~^^^^^^^^^^^^^^^^^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/map_integrands.py", line 87, in <lambda>

    lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type

                 ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag

    (result,) = map_expr_dags(

                ~~~~~~~~~~~~~^

        function, [expression], compress=compress, vcache=vcache, rcache=rcache

        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    )

    ^

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/corealg/map_dag.py", line 101, in map_expr_dags

    r = handlers[v._ufl_typecode_](v)

  File "/home/isaac/.conda/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algorithms/apply_restrictions.py", line 61, in _require_restriction

    raise ValueError(f"Discontinuous type {o._ufl_class_.__name__} must be restricted.")

ValueError: Discontinuous type Jacobian must be restricted.

Some more information, about the mesh, it comes from generating it from this file.geo:

// Gmsh project created on Sat Mar  8 22:37:56 2025
SetFactory("OpenCASCADE");
//+
Mesh.CharacteristicLengthMax = 0.0009;
Mesh.CharacteristicLengthMin = 0.0001;

Disk(1) = {0, 0, 0, 0.05, 0.05};
//+
Disk(2) = {0, 0, 0, 0.035, 0.035};


BooleanDifference{ Surface{1}; Delete;}{ Surface{2}; Delete;}
// Surface 3 is the full annulus. I need to split it into 2.

Point(4) = {0, 0.05, 0};
//+
Point(5) = {0, 0.035, 0};
//+
// Origin.
Point(6) = {0, 0, 0};


Line(4) = {2,3};
Line(5) = {4,5};

Circle(6) = {5, 6, 2};

Circle(7) = {4, 6, 3};
Curve Loop(8) = {5,6,4,7};




Plane Surface(9) = {8};

// 3-quarter surface
BooleanDifference{ Surface{1}; Delete;}{ Surface{9};}

BooleanFragments{ Surface{1}; Surface{9}; Delete; }{ }

Physical Surface("three_quarter_material", 20) = {1};
Physical Surface("quarter_material", 21) = {9};
Physical Curve("theta_0", 22) = {4};
Physical Curve("theta_pi_over_2", 23) = {5};
Physical Curve("r_inner_quarter", 24) = {6};
Physical Curve("r_outer_quarter", 25) = {7};
Physical Curve("r_inner_three_quarter", 26) = {8};
Physical Curve("r_outer_three_quarter", 27) = {9};

Maybe it is faulty but I don’t see it?
I read it with


class TheMesh:
    def __init__(self, mesh_file="meshes/quarter_annulus.msh") -> None:
        self.mesh, self.cell_markers, self.facet_markers = gmshio.read_from_msh(mesh_file, MPI.COMM_WORLD, gdim=2)
        self.x = SpatialCoordinate(self.mesh)
        self.dx = Measure("dx", domain=self.mesh, subdomain_data=self.cell_markers)
        self.ds = Measure("ds", domain=self.mesh, subdomain_data=self.facet_markers)
        self.dS = Measure("dS", domain=self.mesh, subdomain_data=self.facet_markers)
        self.n = FacetNormal(self.mesh)

My whole code is kind of gigantic, the MWE would be big too, since I have to recreate the tensor part, the mixed element part, etc.
I realize it’s very hard to offer guidance, but any idea is still welcome.

dokken · May 25, 2025, 9:15am

It doesn’t seem like you have restricted the FacetNormal (n) in either case

As you are working on an interior facet (dS) measure, the normal can go in either direction.

For one-sided integrals, see for instance

github.com/jorgensd/dolfinx-tutorial

Add one-sided integration example as example of custom integration

opened 08:04PM - 28 Nov 23 UTC

jorgensd

```python # Computing a onesided integral over an interior facet with consisten…t orientation # Enabled by: https://github.com/FEniCS/dolfinx/pull/2269 # Copyright 2023 Jørgen S. Dokken # SPDX: MIT import dolfinx import numpy as np import ufl from mpi4py import MPI mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet) # Create connectivties required for defining integration entities tdim = mesh.topology.dim fdim = tdim - 1 mesh.topology.create_entities(fdim) mesh.topology.create_connectivity(fdim, tdim) mesh.topology.create_connectivity(tdim, fdim) # Get number of cells on process cell_map = mesh.topology.index_map(tdim) num_cells = cell_map.size_local + cell_map.num_ghosts # Create markers for each size of the interface cell_values = np.ones(num_cells, dtype=np.int32) cell_values[dolfinx.mesh.locate_entities( mesh, tdim, lambda x: x[0] <= 0.5+1e-13)] = 2 ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange( num_cells, dtype=np.int32), cell_values) facet_map = mesh.topology.index_map(fdim) num_facets = facet_map.size_local + facet_map.num_ghosts # Create facet markers facet_values = np.ones(num_facets, dtype=np.int32) facet_values[dolfinx.mesh.locate_entities( mesh, fdim, lambda x: np.isclose(x[0], 0.5))] = 2 ft = dolfinx.mesh.meshtags(mesh, fdim, np.arange( num_facets, dtype=np.int32), facet_values) # Give a set of facets marked with a value (in this case 2), get a consistent orientation for an interior integral facets_to_integrate = ft.find(2) f_to_c = mesh.topology.connectivity(fdim, tdim) c_to_f = mesh.topology.connectivity(tdim, fdim) # Compute integration entities for a single facet of a cell. # Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process # local_facet_index is the local indexing of a facet for a given cell integration_entities = [] for i, facet in enumerate(facets_to_integrate): # Only loop over facets owned by the process to avoid duplicate integration if facet >= facet_map.size_local: continue # Find cells connected to facet cells = f_to_c.links(facet) # Get value of cells marked_cells = ct.values[cells] # Get the cell marked with 2 correct_cell = np.flatnonzero(marked_cells == 2) assert len(correct_cell) == 1 # Get local index of facet local_facets = c_to_f.links(cells[correct_cell[0]]) local_index = np.flatnonzero(local_facets == facet) assert len(local_index) == 1 # Append integration entities integration_entities.append(cells[correct_cell[0]]) integration_entities.append(local_index[0]) # Create custom integration measure for one-sided integrals breakpoint() ds = ufl.Measure("ds", domain=mesh, subdomain_data=[ (8, np.asarray(integration_entities, dtype=np.int32))]) n = ufl.FacetNormal(mesh) x = ufl.SpatialCoordinate(mesh) # Exact integral is [y/2**2]_0^1= 1/2 L = ufl.dot(ufl.as_vector((x[1], 0)), n)*ds(8) L_compiled = dolfinx.fem.form(L) print( f"Correct integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L_compiled), op=MPI.SUM)}") # Create reference implementation where we use a restricted two-sided integral with no notion of orientation dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft, subdomain_id=2) n = ufl.FacetNormal(mesh) x = ufl.SpatialCoordinate(mesh) L2 = ufl.dot(ufl.as_vector((x[1], 0)), n("+"))*dS L2_compiled = dolfinx.fem.form(L2) print( f"Wrong integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_compiled), op=MPI.SUM)}") ``` Source: https://fenicsproject.discourse.group/t/wrong-facetnormal-vector-on-internal-boundaries/12887/2?u=dokken

raboynics · May 25, 2025, 1:34pm

Thanks a lot dokken, that worked like a charm!

Topic		Replies	Views
Calculate flux over boundary variational formulation	6	481	November 20, 2023
Integrate over edges variational formulation	10	1516	July 30, 2019
Compute a boundary flux after solving the PDE dolfinx	3	231	June 12, 2024
Unable to compute the curl of a vector field dolfinx	21	253	June 26, 2024
DG variational form	0	667	May 6, 2019

Integrate a flux over an internal curve

Related topics