I have a UnitCubeMesh(2,2,2) and I wish to average the dof spatial coordinates for each facet on the outer surface of the mesh in order to get the coordinates of the center of each facet. From there I also wish to take the three vertices of the triangular facets and determine the normal vector via a cross produce (ie for three vertices (a, b and c, the normal n could be calculated by n = (a-b)x(a-c) where the subtraction is element-wise.)
What I would like to know is, is there a way to be able to average the dof coordinates for each facet on the outer surface of the cubic mesh. Then be able to find the desired normal vectors on each facet via the vertices as described earlier?
Note that this cannot be done via FacetNormal as this produces normal vectors that point away from the center of the mesh (see image).
My current code can produce the dof coordinates and produces the normal vectors via FacetNormal, but this obviously should be changed.
from dolfin import *
import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
mesh = UnitCubeMesh(2,2,2)
Space = FunctionSpace(mesh,'P', 1)
Temperature = project(Expression("x[0]+0.1*x[1] + 0.1*x[2]", degree=3), Space) # Temperature solution across mesh
x_Space = Space.tabulate_dof_coordinates().reshape((-1,mesh.geometry().dim()))
T_values = abs(Temperature.vector().get_local().reshape((-1,1))) # obtain the local temperature values at the correponding dof coordinates
# OLD METHOD OF OBTAINING NORMAL VECTORS
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*ds
l = inner(n, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
num_dofs_per_component = int(V.dim()/V.num_sub_spaces())
num_sub_spaces = V.num_sub_spaces()
vector = np.zeros((num_sub_spaces, num_dofs_per_component))
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
x = V.sub(0).collapse().tabulate_dof_coordinates()
for i in range(num_sub_spaces):
vector[i] = nh.sub(i, deepcopy=True).vector().get_local()
vector = vector.T
for coord, vec, temp in zip(x_Space, vector, T_values):
print(coord, vec, temp)
File("nh.pvd") << nh
filename = "output.txt"
data = zip(x_Space, vector, T_values)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter3D(x_Space[:,0],x_Space[:,1],x_Space[:,2])
ax.quiver(x_Space[:,0],x_Space[:,1],x_Space[:,2], vector[:,0],vector[:,1],vector[:,2], length=0.1, normalize=True)
ax.set_xlim(-0.1,1.1)
ax.set_ylim(-0.1,1.1)
ax.set_zlim(-0.1,1.1)
plt.savefig('sourceplot.pdf', bbox_inches='tight', format='pdf')
plt.savefig('sourceplot.png', dpi=300, bbox_inches='tight', format='png')
The plot below shows the current location of the dof coordinates that I would like to average over for each facet and the incorrect normal vectors output by FacetNormal
