Plus and negative side of an internal facet

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(``-")?

I am basically enforcing the flux condition according to Imposing a discontinuity at interface using DG method - #19 by smesc

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 “-”.

You may want to have a look at dolfinx/python/test/unit/fem/test_assemble_domains.py at b16ef8abc20d70c28d8f1850fe42a0ee0bcf5106 · FEniCS/dolfinx · GitHub . Rather than say + and -, it allows to choose yourself which faces to integrate on.

2 Likes

Another example can be found here: Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub

1 Like

This is what I needed.

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.

what’s wrong with defining your own jump function that reads n[+] - n[-], or the other way around?

1 Like

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 have checked ufl/ufl/algorithms/apply_restrictions.py at main · FEniCS/ufl · GitHub but cannot figure out how the “+” and “-” sides get selected and whether at the point of that selection, there is access to information about which subdomains lie on each side of the facet.

Basically, what I require is the old fix Fix consistent restriction of one-sided interior facet integrals by jorgensd · Pull Request #2127 · FEniCS/dolfinx · GitHub. I feel that the PR introducing that change was prematurely closed. Situations where one needs to integrate a function whose input values is the jump of a value in a specified direction across the interface cannot be solved by the one-sided integral fix.

Though the fix by Dokken was deemed implicit, I believe it would be tremendously useful for folks whose problem required such a feature and would not have negatively impacted anyone who did not need it. Therefore, both Dokken’s fix and the alternate one-sided integration can both exist. By the description of the test at dolfinx/python/test/unit/fem/test_assemble_domains.py at b3cb5c8127274aaeec2ddea88340651b6068df7e · FEniCS/dolfinx · GitHub, I believe this would solve my problem.

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

Thanks, I think I’ve gotten it now. I noticed that I can mark only one boundary with such a method. For example, I can do:

int_facet_domains = [(3, int_facet_domain)]
dS_manual = ufl.Measure("dS", subdomain_data=int_facet_domains, domain=msh)

but cannot do

int_facet_domains = [(3, int_facet_domain), (5, other_int_facet_domain)]
dS_manual = ufl.Measure("dS", subdomain_data=int_facet_domains, domain=msh)

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.