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.

What you sketch out here should work, could you comment on what is not working, and post an acompanying reproducible example

Hi there , I am facing the same problem as you indicated here before. Could you please share you final code for the solution?

My specific usage required me to reorder facets at only one boundary and it was not necessary to add more than one tag to dS. Poisson submesh DG coupling · GitHub has some utility functions that you may find helpful.