From e8c5c2fe69e65c03da4ce26c949d6c572079edb2 Mon Sep 17 00:00:00 2001
From: Florian <florian.gacon@inria.fr>
Date: Tue, 15 Sep 2020 16:10:16 +0200
Subject: [PATCH] Finalize boundary_normal
---
src/bvpy/domains/geometry.py | 21 ++++++++++++++++++---
1 file changed, 18 insertions(+), 3 deletions(-)
diff --git a/src/bvpy/domains/geometry.py b/src/bvpy/domains/geometry.py
index cb5aed3..8c4d7a2 100644
--- a/src/bvpy/domains/geometry.py
+++ b/src/bvpy/domains/geometry.py
@@ -263,7 +263,9 @@ def boundary_normal_bkp(mesh):
def boundary_normal(mesh):
mesh.init()
bnd_facets = get_boundary_facets(mesh)
- normals = list()
+ facets_index = [f.index() for f in bnd_facets]
+ facet_normals = np.zeros((mesh.num_entities(1), 3))
+ normals = np.zeros((mesh.num_entities(0), 3))
for fct in bnd_facets:
cell = [c for c in fe.cells(fct)][0]
@@ -274,9 +276,22 @@ def boundary_normal(mesh):
vec = midpoint_vertex-midpoint_facet
vec = vec/np.linalg.norm(vec)
n = np.cross(normal_cell, vec)
- normals.append(n)
+ facet_normals[fct.index()] = n
+
+ for fct in bnd_facets:
+ for vtx in fe.vertices(fct):
+ norm = np.array([facet_normals[f.index()]
+ for f in fe.facets(vtx)
+ if f.index() in facets_index])
+ norm = np.sum(norm, axis=0)
+ norm = norm / np.linalg.norm(norm)
+ normals[vtx.index()] = norm
+
+ V = fe.VectorFunctionSpace(mesh, 'P', 1, dim=3)
+ u = fe.Function(V)
+ d2v = fe.dof_to_vertex_map(V)
+ u.vector().set_local(normals.flatten()[d2v])
- normals = np.asarray(normals)
return normals
# geometry.py ends here
--
GitLab