Prescribing discontinuity using hybridizable discontinuous Galerkin (HDG)

General
1

Using symmetric interior penalty DG, I have solved a Poisson problem \nabla\cdot(-\kappa\nabla u)=0 where I prescribe a discontinuity across an interface between two subdomains of the form
-\kappa\nabla u\cdot \pmb{n}|_+=\alpha(u_+ - u_-) where “+” refers to the right side and “-” to the left side. The domain can be thought of as a a unit square vertically divided into two halves.

As my goal is to solve the problem in a much larger 3D domain with at least two subdomains, I’m exploring hybridizable discontinuous Galerkin because of it’s superconvergence and sparsity over usual DG. I have compiled main branch of FEniCSx and can run the demo dolfinx/python/demo/demo_hdg.py at main · FEniCS/dolfinx · GitHub.

I have separated out the interface integration facets to have:

ds_c = ufl.Measure("ds", subdomain_data=[(2, non_interface_facets),(3, left_interface_facets), (4, right_interface_facets) ], domain=msh)

where non_interface_facets is an array of cell_id1, local_facet_id1,... etc for facets that are not at the interface and left_interface_facets and right_interface_facets are similarly ordered for facets that are at the interface.

When I prescribe the discontinuity that involves ds_c(3) and ds_c(4) at the interface, I get the error at the point of computing the Jacobian. Basically, I cannot integrate on a subset of the cell boundaries.

Traceback (most recent call last):
  File "/home/lmolel/work/ssb/tertiary_current_hdg.py", line 318, in <module>
    J10 = fem.form(jac10, entity_maps=entity_maps)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/lmolel/python3-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 280, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/home/lmolel/python3-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 275, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/home/lmolel/python3-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 183, in _form
    assert all([d is data[0] for d in data])
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError

Any idea of how I can solve the issue?

My Minimum Working Example is rather long because of the code in extracting the cell integration facets and that I use a custom Newton solver (from Poisson submesh DG coupling · GitHub), but I can provide my code.