Bug? report: specialcf.normal orientation

Main Forum Old-Forum - do not post here
1

Hello everybody,

I’m curious whether it is a bug and/or known, that specialcf.normal(3) doesn’t necessarily point outwards for non-trivial geometries.

Attached you find a minimal (non) working example.
Best, Carl

2

https://ngsolve.org/media/kunena/attachments/972/normal_bug.py [attachment=undefined]normal_bug.jpg[/attachment] [attachment=undefined]normal_bug.jpg[/attachment]

3

Attachment doesn’t work…
Here the example:

import netgen.gui
from ngsolve import Draw, Mesh
from netgen.csg import OrthoBrick, Pnt, CSGeometry

big = OrthoBrick(Pnt(0,0,0), Pnt(4,4,4))
small = OrthoBrick(Pnt(1, 1, 2), Pnt(3, 3, 4))

geo = CSGeometry()
geo.Add(big-small)
ngmesh = geo.GenerateMesh(maxh=0.5)
mesh = Mesh(ngmesh)

from ngsolve import specialcf
Draw(specialcf.normal(3), mesh, “nu”)
input(“drawn; enable Visual >> Draw Surface Vectors”)

4

Hi Carl,

that’s not a bug. The normal vector is defined outward for you ‘small’ brick.
Later you subtract it to get the subdomain, this does not change the normal vector.

You can use this code snipped to obtain the sub-domain numbers in front / behind the individual faces (taking care of 1-based enumeration in Netgen).

for i in range(ngmesh.GetNFaceDescriptors()):
    fd = ngmesh.FaceDescriptor(i+1)
    print (fd)
    print (fd.domin, fd.domout)

Joachim

5

Hey Joachim,

thanks for the quick reply.

Naturally, I suppose one wants to consider (orientation of) normal vectors with respect to the meshed domain. Of course, there opinions might differ, BUT:

Normal vectors not pointing outwards everywhere, (I suppose) is a consequence of boundary faces not being oriented consistently. And therefore, this simple example of big_brick minus small_brick above, breaks your tutorial on Working with meshes:

Since boundary faces are not oriented consistently, creating a volume mesh from the boundary mesh fails with the correct error message

ERROR: Edge 307 - 308 multiple times in surface mesh ... same error for some other edges ... netgen.libngpy._meshing.NgException: Stop meshing since surface mesh not consistent
Since boundary face orientation is incosistent, alle edges belonging to the boundary of the top face of the small brick appear twice in the mesh – instead of once, plus once in reverse direction!

Using the mesh above, and a snippet from the tutorial, reproduce error via:

[code]from ngsolve import Mesh as ngsMesh
from netgen.csg import OrthoBrick, Pnt, CSGeometry

big = OrthoBrick(Pnt(0,0,0), Pnt(4,4,4))
small = OrthoBrick(Pnt(1, 1, 2), Pnt(3, 3, 4))

geo = CSGeometry()
geo.Add(big-small)
ngmesh = geo.GenerateMesh(maxh=0.5)
mesh = ngsMesh(ngmesh)

new from here

from netgen.meshing import FaceDescriptor, Element2D, Mesh
new_ngmesh = Mesh()
fd_outside = new_ngmesh.Add(FaceDescriptor(bc=1,domin=1,surfnr=1))

pmap = {}
for e in mesh.ngmesh.Elements2D():
for v in e.vertices:
if (v not in pmap):
pmap[v] = new_ngmesh.Add(mesh.ngmesh[v])

for e in mesh.ngmesh.Elements2D():
new_ngmesh.Add(Element2D(fd_outside, [pmap[v] for v in e.vertices]))

new_ngmesh.GenerateVolumeMesh()
[/code]

If it’s not a bug, there’s definitely an inconsistency.
(Nothing I couldn’t work around with the badhack of changing the list-order of

[pmap[v] for v in e.vertices]

for some boundary faces. Anyhow, I still wanted to report this, since I could imagine other problems implied by the incosistent orientation.)

Best, Carl

6

yes, that’s a good example where you can use the code snipped from my previous answer:

The surface element has an index pointing into the array of FaceDescriptors.
Whenever the FaceDescriptor has domout = 0, you use the orientation as is, if fd.domin = 0 you flip the orientation.

If both are non-zero, you found an element belonging to an internal interface, which we don’t want to copy to the new mesh.

7

I have a similar question to this. Does one really have to do the hack proposed above to get a consistent outward pointing normal vector on the boundary? Shouldn’t there be a method to make the normal vector outward pointing everywhere?

In some other threads, it was claimed that specialcf.normal gave the outward pointing normal on an element. But when I interpolate specialcf.normal so I can visualize it, it doesn’t do that on all parts of the boundary.

8

There is a quite recent feature that you can do:

specialcf.normal(mesh.Materials("mat"))

to get an outward pointing normal for mat. For boundaries inside mat or not adjacent to mat this normal is also randomly facing (as the occ face is oriented).

best

9

Thanks a lot! For interior interfaces, that makes sense that the normal is random facing (the normal direction is more arbitrary).