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?
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?
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?
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…
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.)
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?
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?
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().
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?
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
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)?
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.
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?
Hi Hongzhi, when skeleton=True is used in an integrator, one will want to supply a BitArray related to the facets in particular, so in this regard I think it should be your V1 script, and not V2. I implemented methods like this a few years ago (Higher order Discontinuous Galerkin methods for the Laplace-Beltrami problem on unfitted smooth surfaces | GRO.publications), and append a script from that time (syntax partially outdated) which might be helpful for reference. When comparing the two, I noted that you used the H1 space. Depending on whether you want to perform a test or actually go ahead with “full DG”, you might want to start a comparison of the two scripts/ further debugging around this point. Best regards, Fabian low_order_dg2d.py (4.9 KB)
I hope you are doing well. Thank you for your reply with detailed explanation and helpful script!
I have a follow-up question about the higher order geometric approximation. I am trying to implement this with respect to my exampleV1, but when I use different mesh sizes I almost get the same error (I attach the script, or for some other example even get larger error for smaller mesh size). Could you help me with this?
Hi Hongzhi, interesting, of course the method with higher order geometry should converge with higher order as well. I attach one of my scripts as well for higher order and 3D ambient geometry, so you can compare. On a brief look for comparison, I noticed apart from the L2/ H1 space difference, we have different construction principles for e.g. the edge normal. Not sure whether these are ultimately equivalent or there are subtle differences, but that would be among the things I would have a look next… Best regards, Fabian tracedg3d.py (5.3 KB)