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


211
views
0
6 months ago by
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
6 months ago by
micro  
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.

Similar posts:
Search »