I am currently implementing a cut DG method on a surface embedded in either a 2D og a 3D mesh. The surface is discretized using the elements in the background mesh, and I would like to access the normals that are coplanar to the intersection between an element T and the discrete surface Gamma_h (see image). I can compute the surface normal by

n = 1.0 / sqrt(InnerProduct(grad(lsetip), grad(lsetip))) * grad(lsetip)

but is there any convenient way to compute the co-normals?

I don’t know what “convenient” way you are searching for, but it should be quiet obvious.

You have the facet normals from ngsolve in the embedded space. You can use the cross product to obtain a co-normal vector. Here is an example snippet that we use for our CutDG calculations:

This is (almost) what I have done in my code as well, but I have some questions.

How is the facet normal defined inside the element? It does not really matter, since the normalised projection should give the same output regardless, but I am just curious.

What happens if the facet normal and the surface normal are parallell? Then the projection gives 0, right?

Does this work when the background mesh is 3D? The projection does not necessarily give a vector which is orthogonal to E when E is a line.

well, the facet normal is the unit outer normal vector. If you are on one element and use an integration rule that acts on the boundary of the element, then the normal points away from that element. If you are using skeleton-integrators there is one unique direction.

Assuming you use a P1 level set function and straight elements, the facet normal and the surface normal can only be parallel if the level set function is zero on all three vertices. In ngsxfem (which you are using?!) we avoid handling these situation by perturbing the P1 levelset function whenever it hits zero at a vertex. This is done in InterpolateP1. The perturbation shifts every zero at a vertex to a small (around machine precision, but you can also set this yourself) non-zero value. The geometrical error is typically negligible, but this step simplifies dealing with different cut topologies.

This works in 3D, yes. And the result is orthogonal to E (if E ist the intersection of the zero level set and the facet). Note that the conormal is a linear combination of the surface normal nphi and the facet normal nF which are both orthogonal to E. Hence, the conormal is also orthogonal to E.