Questions about Jump on Interior Edge and Interior Face in ngsxfem

NGSolve
1

Hi Everyone,

I hope you are doing well.

I have some questions about ngsxfem using in 3D.

  1. How to do integration on the interior edge which is the intersection between the discrete surface and the active background mesh? And how to find the co-normal to the interior edge? Do we need the normal to the discrete surface and the tangential to the interior edge, then take the cross product?
  2. How to do integration on the interior face which is the intersection between the active tetrahedrons? And how to find the normal to the interior face?
  3. How to construct the jump for the tangential gradient over the interior edge, the jump for the full gradient over the interior face, and the jump for the Hessian matrix over the interior face, respectively?

Thank you very much.

Best regards,

Hongzhi

2

Hi Hongzhi,

1a. You do the integration on the surface part of an interior facet through
dCut(lsetp1, IF, skeleton=True, definedonelements=...)
1b. in 3D you can take compute everything based on the normal of the facet of the background mesh (nF=specialcf.normal(..)) and the normal of the level set function nl=Normalize(grad(lsetp1)). If you take t=Normalize(Cross(nF,nl)) you obtain the tangential vector of the cut edge. Afterwards you can compute the co-normal through nc=Normalize(Cross(nl,t)).
2. If you integrate on the whole interior face (= facet here), you can work with the standard ngsolve skeleton-integral and combine it with the right marking of facets in a definedonelements=…-flag (which you should also use before for the edge integrals). To obtain the right marking take a look at our standard CutFEM examples. The normal of the facet is simply the standard ngsolve normal (nF=specialcf.normal(..))
3. The jumps are also created as usual in NGSolve. simply work with the evaulation of the differential operators and combine with the ....Other()-operation. To obtain the tangential gradient multiply the gradient with the tangential direction, etc…

Best,
Christoph

3

Hi Christoph,

Thank you very much for the detailed reply!

I just have some follow-up questions.

  1. What is the right way to define the Laplace-Beltrami operator on the discrete surface? I have tried the following, but it does not work. (I am sorry that I do not know how to type code in the reply.)

Ph = Id(3) - OuterProduct(nl, nl)
def LB(u):
lb = Ph * u.Operator(“hesse”) - InnerProduct(nl, grad(u)) * grad(nl)
return lb[0, 0] + lb[1, 1] + lb[2, 2]

  1. For the jump, I know we can use “grad(u) - grad(u.Other())” for the full gradient over the interior face. Is it correct to use “u.Operator(“hesse”) - u.Other().Operator(“hesse”)” for the Hessian matrix over the interior face?
  2. Do we have the jump for the normal, tangential, or co-normal vector, e.g., nl.other(), nc.other()? If not, is there any way to define the average or jump for the Laplace-Beltrami operator?

Best regards,

Hongzhi

4

Hi,

nl is not a Gridfunction and you can hence not take the gradient. You can try interpolating nl into a GridFunction to take the gradient then. Be careful, if you treat mesh deformation (are you using it?) correctly in your context.

Yes, you can use (u.Operator("hesse") - u.Operator("hesse").Other())…

For a CoefficientFunction cf, you can also use cf.Other().

Best,
Christoph

5

Hi Christoph,

Thank you very much for the clarification!

In your answer for the first question, you mentioned I need to be careful when I use mesh deformation. I feel a little confused about this. Could you give some details or an example?

Best regards,

Hongzhi

6

If you are interested in first order polynomial FEM anyway, ignore that comment. Otherwise you want to have higher order accuracy in the geometry handling. And for that we use a mesh deformation strategy similar to isoparametric FEM.
Take a look at this tutorial for more details: 8.2 Integration on level set domains — NGS-Py 6.2.2501 documentation

Best,
Christoph

7

Hi Christoph,

Thank you for the reference. If we use the deformation, how can we get the normal vector to the higher order discrete surface (or geometric approximation)? Can we use the following, e.g., for the 2nd order, (I am not sure if the syntax is correct)?

lsetp2 = GridFunction(H1(mesh, order = 2))
InterpolateToP2(levelset, lsetp2)
nl = Normalize(grad(lsetp2))

Or, should we directly get the information from the deformation?

Best regards,

Hongzhi

8

You can simply use the lsetp1 and normalize. One you apply the mesh deformation this will do the correct thing. The gradient of the p1-function ist mapped according to the deformation which will yield the higher order accurate normal direction. Normalize then acts on this vector.

Best,
Christoph

9

Hi Christoph,

I see. Thank you very much for your help and I really appreciate it!

Best regards,

Hongzhi

11

Hi Christoph,

I hope you are doing well.

I attach two versions of code in which I am trying to implement the method (2.25)-(2.29) in the paper https://doi.org/10.1093/imanum/drv068 whose formulation involves the stuff I asked before, but the results from the two versions are very different. The main difference happens when I define the measures on the surface, on the surface edge, and the interior face around the surface (i.e., ds, de, and df in the code). I am not sure which one is correct. Could you give me some advice?

Thank you very much.

Best regards,

Hongzhi
exampleV2.py (3.9 KB)
exampleV1.py (3.9 KB)