I am a solving a problem using symmetric interior penalty discontinuous Galerkin and I am running into convergence issues that I am attributing to labeling of internal facets.

I need to impose a flux expression on the internal boundary between the pink and the turquoise colored subdomains.

The expression for flux that I need to impose is: -\kappa \nabla u \cdot n_{\mathrm{turquoise}} = \gamma (u_{\mathrm{turquoise}} - u_{\mathrm{pink}} - \mathrm{constant})

I obtain convergence when the boundary between the pink and turquoise domains is a straight line.

My solution involves replacing the jump of u
[\![ u ]\!] in the internal facet integral with (\mathrm{constant} - \frac{\kappa}{\gamma}\nabla u \cdot n_{\mathrm{turquoise}})n_{\mathrm{turqoise}} .

In such a boundary, is there a way to know which is the subdomain that provides u(``+") and n(``+"), and which subdomain provides u(``-") and n(``-")?

Further investigation shows that it is related to Facets orientation in inner surface with fenicsX - #3 by RaphaelC, but I cannot figure out what the solution is. At the interface between the subdomains, the orientation needs to be definite. For example, I may want the turquoise domain to be “+” and the pink domain to be “-”.

Given now that I have control over direction of ordering, how do I express a function such as ufl.sqrt([\![ u ]\!]) where [\![ u ]\!] is a jump that is non-zero? I am unclear how to combine the labeling that I have made according to Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub to generate a function that takes the square root of the local jump and integrates once over the shared interface.

Nothing wrong with that. I am trying to figure that out from the ufl library, but it’s a slow process and any tips you provide is greatly appreciated. What I need is for “+” side to be on the side of domain with higher cell tag and “-” to be on the side of domain with lower cell tag if the facet is an internal boundary shared by two subdomains. For facets not at the internal boundary between two subdomains, I can use default behavior of “+” and “-” sides.

I wouldn’t necessarily rely on us resurrecting a fix that was dismissed years ago. OK, so since indeed you can’t rely on + and -, can’t you expand the definition of jump, use its linearity, and write the two integrals separately?

In the integral tuple we create with integration entities, (cell_0, local_ facet_index_ cell_0, cell_1, local_facet_index_cell_1) the plus restriction is the first cell (0), while the negative restriction is cell 1. thus by manually packing them in that order you can get a directional jump

Marking only 1 internal boundary where I want deterministic “+” and “-” sides is sufficient for my current problem, but it’d be great to learn if it possible (if so how) to label multiple internal boundaries for deterministic “+” and “-” sides.

I appreciate both of you for your help and patience.