Bug? report: specialcf.normal orientation

More
4 years 1 month ago - 4 years 1 month ago #3219 by cpfeiler
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
Last edit: 4 years 1 month ago by cpfeiler. Reason: somehow attachment got lost first try
More
4 years 1 month ago #3220 by cpfeiler

Attachment not found

[attachment=undefined]normal_bug.jpg[/attachment] [attachment=undefined]normal_bug.jpg[/attachment]
More
4 years 1 month ago #3221 by cpfeiler
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")
More
4 years 1 month ago #3223 by joachim
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).
Code:
for i in range(ngmesh.GetNFaceDescriptors()): fd = ngmesh.FaceDescriptor(i+1) print (fd) print (fd.domin, fd.domout)

Joachim
More
4 years 1 month ago #3224 by cpfeiler
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
Code:
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()

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
Code:
[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
More
4 years 1 month ago #3225 by joachim
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.
Time to create page: 0.102 seconds