When calculating the boundary mesh normal vectors the sign of the direction is incorrect

1

25 minutes ago by

Terrence

Hello all,
I need to calculate the normal vectors of the boundary mesh for my model. When I do this the orientation of the vector is correct but the sign of the vectors seem to be randomly positive or negative. My model involves a somewhat complicated geometry so testing each cell would be difficult.

Below is a short working example of the problem with a unit cube. The code prints the value of the z component of the normal vector of the two cells on the top face.

from dolfin import *
import time
import numpy as np
import neumann
x_min = y_min = z_min = 0.0
x_max = y_max = z_max= 1.0
nx = ny =1
nz =1
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)
bmesh = BoundaryMesh(mesh,'exterior')
for i, cell in enumerate(cells(bmesh)):
    if cell.midpoint().z()==z_max:
       print('z-direction top surface normal',cell.cell_normal().z())

with the following result:
z-direction top surface normal -1.0
z-direction top surface normal 1.0
Both should be 1.0.

Any help would be greatly appreciated!

Cheers,
Terrence

2

Hi, consider the following

from dolfin import *

mesh = BoxMesh(Point(-1, -1, -1), Point(1, 1, 1), 4, 4, 4)
bmesh = BoundaryMesh(mesh, 'exterior')
# Init w.r.t outward normal
bmesh.init_cell_orientations(Expression(('x[0]', 'x[1]', 'x[2]'), degree=1))

cell_orientations = bmesh.cell_orientations()
for o, cell in zip(cell_orientations, cells(bmesh)):
    x, y, z = cell.midpoint()

    if near(abs(x), 1):
        want = Point(-1 if near(x, -1) else 1, 0, 0)
    elif near(abs(y), 1):
        want = Point(0, -1 if near(y, -1) else 1, 0)
    else:
        want = Point(0, 0, -1 if near(z, -1) else 1)

    cn = cell.cell_normal()
    # Reorient
    if o: cn *= -1

    assert (cn - want).norm() < 1E-14
3

Thanks for the reply. This seems to work well for simple geometry where the boundary coordinates are known (in this case 1/-1/0). I need a method that can do this for a complex geometry. I also don’t fully understand the use of the cell_orientations in your code.

Thanks,
Terrence

4

This example shows that i) the mesh due to BoundaryMesh is not oriented and ii) how to get the orientation using cell.normal and an orientation flag which is computed once you orient the mesh with respect to some normal. The example is simple but the idea should be generic.

1 Like
5

Running this on the unit cube doesn’t produce the correct results as seen below. I understand the local normal definition is not the same as the global normal. Maybe I am not fully understanding the clockwise counterclockwise orientation given some vector defined in Expression() or how to implement this.
from dolfin import *
import time
import numpy as np
import neumann
x_min = y_min = z_min =0.0
x_max = y_max = z_max= 1.0
nx = ny =1
nz =1
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)
bmesh = BoundaryMesh(mesh,‘exterior’)
bmesh.init_cell_orientations(Expression((‘x[0]’, ‘x[1]’, ‘x[2]’), degree=1))
cell_orientations = bmesh.cell_orientations()
for i in range(len(cell_orientations)):
print(‘i=’,i,‘cell orient’,cell_orientations[i])
j = 0
for o, cell in zip(cell_orientations, cells(bmesh)):
x, y, z = cell.midpoint()
cn = cell.cell_normal()
# Reorient
print(o)
print(x,y,z)
print(cn.x(),cn.y(),cn.z())
if o: cn *= -1
print(cn.x(),cn.y(),cn.z())
j+=1
with results:

0
0.6666666666666666 0.3333333333333333 0.0
0.0 0.0 1.0
0.0 0.0 1.0
0
0.6666666666666666 0.0 0.3333333333333333
0.0 -1.0 0.0
0.0 -1.0 0.0
0
0.3333333333333333 0.6666666666666666 0.0
0.0 0.0 1.0
0.0 0.0 1.0
0
0.0 0.6666666666666666 0.3333333333333333
1.0 0.0 0.0
1.0 0.0 0.0
…
1
0.6666666666666666 0.3333333333333333 1.0
0.0 0.0 -1.0
-0.0 -0.0 1.0

6

Global normal is defined as a position w.r.t an interior point.

from dolfin import *

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 4, 4, 4)
bmesh = BoundaryMesh(mesh, 'exterior')
# Init w.r.t outward normal
bmesh.init_cell_orientations(Expression(('x[0]-0.5', 'x[1]-0.5', 'x[2]-0.5'), degree=1))

cell_orientations = bmesh.cell_orientations()
for o, cell in zip(cell_orientations, cells(bmesh)):
    x, y, z = cell.midpoint()

    if near(abs(x*(1-x)), 0):
        want = Point(-1 if near(x, 0) else 1, 0, 0)
    elif near(abs(y*(1-y)), 0):
        want = Point(0, -1 if near(y, 0) else 1, 0)
    else:
        want = Point(0, 0, -1 if near(z, 0) else 1)

    cn = cell.cell_normal()
    # Reorient
    if o: cn *= -1

    assert (cn - want).norm() < 1E-14, ((cn - want).norm(), (x, y, z), cn.array(), want.array())
7

I am facing the same question and your answer is helpful. But I do not think it is correct when the boundary is concave.

8

Can you post the problematic case? Here is my (working) attempt with convex boundary

from dolfin import *

mesh = Mesh('concave_mesh.xml')
bmesh = BoundaryMesh(mesh, 'exterior')
# Inside point
bmesh.init_cell_orientations(Expression(('x[0]+0.5', 'x[1]', 'x[2]-0.5'), degree=1))

cell_orientations = bmesh.cell_orientations()
for o, cell in zip(cell_orientations, cells(bmesh)):
    x, y, z = cell.midpoint()

    if near(z, 0):
        want = Point(0, 0, -1)
    elif near(z, 1):
        want = Point(0, 0, 1)
    elif near(y, -1):
        want = Point(0, -1, 0)
    elif near(y, 1):
        want = Point(0, 1, 0)
    elif near(x, -1):
        want = Point(-1, 0, 0)
    elif near(x, y):
        want = Point(1/sqrt(2), -1/sqrt(2), 0)
    else:
        want = Point(1/sqrt(2), 1/sqrt(2), 0)
        
    cn = cell.cell_normal()
    # Reorient
    if o: cn *= -1

    assert (cn - want).norm() < 1E-14, ((cn - want).norm(), (x, y, z), cn.array(), want.array())

The mesh comes from gmsh

// concave_mesh.geo
// gmsh -2 concave_mesh.geo
// dolfin-convert concave_mesh.msh concave_mesh.xml
size = 1;

Point(1) = {-1, -1, 0, size};
Point(2) = {1, -1, 0, size};
Point(3) = {0, 0, 0, size};
Point(4) = {1, 1, 0, size};
Point(5) = {-1, 1, 0, size};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 1};

Line Loop(1) = {1, 2, 3, 4, 5};
Plane Surface(1) = {1};

Extrude {0, 0, 1} { Surface{1}; }
Physical Volume(1) = {1};
9

My teacher wanted me to simulate a liquid drop in electric field. This liquid drop might deform to a shape like number 8. Considering this situation, it is necessary to devise a generic method to calculate outward normal direction for any concave body.

I think I have solved this problem and my idea is quite simple. For a 3D mesh, its boundary is make up of triangles and every triangle is the surface of only one tetrahedron in the mesh. Therefore, original problem is equivalent to calculate the outward normal of the tetrahedrons who contain boundary triangles.

Class MeshEntity and method entities are used when I implemented it. I only have c++ version currentlyhttps://github.com/mapengfei-nwpu/outnormal/blob/master/Untitled-1.cpp, so I do not post my python code today.

Finally, I am sorry for not responding you until now, because I do not receive any notification.