### If statement to calculate the face normal vectors only on the geometric boundary

211

views

0

Hey everyone I figured out a way to explicitly calculate the facet normal vectors but I only want to calculate them on a prescribed boundary. My simple code below:

```
from dolfin import *
import time
import numpy as np
h = [1,1,0]
#define domain
x_min = y_min = z_min = 0.0
x_max = y_max = z_max = 1.0
#define number of elements
nx = ny = 1
nz = 1
#create mesh
mesh = BoxMesh(Point(x_min, y_min, z_min), Point(x_max, y_max, z_max), nx, ny, nz)
# define function space
V = FunctionSpace(mesh, 'P', 1)
#define boundary
def boundary_sides(x, on_boundary):
return near(x[0],x_min) or near(x[0],x_max) or near(x[1],y_max) or near(x[1],y_min)
bmesh = BoundaryMesh(mesh, "exterior", True)
print(bmesh.coordinates())
b = len(bmesh.coordinates())
face_numbers=mesh.num_faces()
print('number of faces', face_numbers)
i = 0
while i < face_numbers:
b_f = Facet(mesh,i)
n = [b_f.normal().x(), b_f.normal().y(), b_f.normal().z()]
#print('type n', type(n))
print('node # =', i,' ','n = ', n)
#print('size n', len(n))
i=i+1
```

This calculates every facet normal which is not what I need.

when I try to only include the boundary mesh:

```
bmesh = BoundaryMesh(mesh, "exterior", True)
print(bmesh.coordinates())
b = len(bmesh.coordinates())
face_numbers=mesh.num_faces()
print('number of faces', face_numbers)
i = 0
while i < face_numbers:
b_f = Facet(bmesh,i)
n = [b_f.normal().x(), b_f.normal().y(), b_f.normal().z()]
```

I get an error:

"

```
Error: Unable to find normal.
*** Reason: Normal vector is not defined in dimension 3 (only defined when the triangle is in R^2.
*** Where: This error was encountered inside TriangleCell.cpp.
```

is there a better way to only calculate the unit normals only on a prescribed boundary facets?

Thanks in advance!

Terrence

Community: FEniCS Project

### 1 Answer

3

Hi, consider

```
from dolfin import *
mesh = UnitCubeMesh(4, 4, 4)
bmesh = BoundaryMesh(mesh, 'exterior')
# bmesh cells to mesh cells
cell_map = bmesh.entity_map(2)
msg = 'Cell %d in bmesh(%d in mesh) has normal %r'
for bci, cell in enumerate(cells(bmesh)):
n = cell.cell_normal().array()
print msg % (bci, cell_map[bci], list(n))
```

Thanks so much for the response this works perfectly!

written
6 months ago by
Terrence

Please login to add an answer/comment or follow this question.