Dear developers,
I am currently using ngsxfem to study unfitted methods and to experiment with writing my own codes. However, I have some questions regarding the usage of dFacetPatch.
In the demos, the ghost penalty term is usually added in the following way:
dw = dFacetPatch(definedonelements=ba_gp_facets, deformation=deformation)
a += gamma_stab / h**2 * (u - u.Other()) * (v - v.Other()) * dw
This works fine for the standard case. The problem arises when I try to implement a higher-order stabilization term, for example the penalty with respect to the jump of the normal gradient:
Ah += gamma_stab * h * ((Grad(u) - Grad(u.Other()))*ne) * ((Grad(v) - Grad(v.Other()))*ne) * dw
I find that this term evaluates to zero. This makes me confused about what exactly is integrated when using dFacetPatch.
Later, I suspected that I might need to use instead:
dw = dCut(lsetp1, NEG, skeleton=True, definedonelements=ba_gp_facets, deformation=deformation)
But I am not sure about the correct domain specification here. In particular, the domain is set to NEG rather than HASNEG, and from the literature I understand that the ghost penalty term should be added on the active mesh rather than on the physical domain.
Could you please clarify the intended usage of dFacetPatch versus dCut(..., skeleton=True, ...) for implementing ghost penalty terms, especially for higher-order cases involving jumps of normal derivatives?
Hi,
Great to see the software being used; welcome to the community! The dFacetPatch symbol is used to implement our patchwise or “direct” version of the Ghost Penalty. See e.g. Eq. (4.43) - (4.46) in my thesis [1], or (4.3), (4.4) in our recent paper [2]. The idea being that also in the high order case the direct volumetric difference between u on “my” side and u on “the facet F neighbour’s” side is assembled. See e.g. Lemma 4.4. in my thesis for the application. Hence, we are often interested in integrating over let’s say \omega_F = T_1 \cup T_2, where F = T_1 \cap T_2, the facet patch, and this is represented in dFacetPatch. It should suffice also for higher order, and we applied it to some problems, but it is hard to generalise about all possible discretisations/ PDEs.
There is the alternative of integrating on the facet F itself and summing over derivatives in accordance with the discretisation order, as e.g. in [3]. (See also Section 4.3 in [4].) To implement that in ngsxfem, you might want to use a differential symbol such as dxtref(mesh, time_order=time_order,skeleton=True, vb=BND, definedonelements=ba_facets) for space-time integrals or just dx(skeleton=True, vb=BND, definedonelements=ba_facets) in a time stepping/ stationary problem application. Might be inconvenient to calculate higher order discrete derivatives though, which is why we would not necessarily recommend this approach for orders like e.g. 6,7,8.
If you say more about the PDE/ discretisation of interest, we might be able to say more in particular.
Best,
Fabian
[1] Higher Order Unfitted Space-Time Finite Element Methods for Moving Domain Problems
[2] [2504.08608] Discretization Error Analysis of a High Order Unfitted Space-Time Method for moving domain problems
[3] https://doi.org/10.1093/imanum/drz021
[4] Cut finite element methods | Acta Numerica | Cambridge Core
1 Like
Dear Fabian,
Thanks a lot for your reply, it really helps me.
I had misunderstood the definition of dFacetPatch, as I initially thought it was meant for integration over faces. Now I understand why my code did not work.
While studying ghost penalty methods, I often find these terms quite challenging. At the moment I am working on a linear elasticity problem in 3D and trying to apply ghost penalty stabilization. However, I observed that the error does not decrease when refining the mesh. My intention was to implement the higher-order discrete derivative approach, which I often see in the literature.
Thank you also for pointing me to your thesis and papers. I will study them carefully to learn more about the patchwise version, and I plan to try this approach in my experiments.
Best regards,
Shuming
Hi Shuming,
Great, very interesting. With regards to linear elasticity, I can also refer to our script [1] from some years ago; might be interesting for comparison. More generally, I think the direct Ghost penalty would be one good candidate stabilisation for such applications.
Best regards,
Fabian
[1] esgi-soft-tissue/elast_unf2_3D.py at main · fabianheimann/esgi-soft-tissue · GitHub