Thank you very much @dokken !!
I have looked at the link you provided. Since i am not dealing with interior facets, I can simply solve over ‘ds’ instead of ‘dS’ and without restricting the normal directions using ‘res’.
Here is the code I used, and this prints the tangent vectors and normal vectors.
In the code, I am printing the facet normals that lie in the bottom (base) of cube. Hence, the normals are supposed to have only ‘z’ component.
from fenics import *
from ufl import Jacobian, replace
from ufl.classes import FacetJacobian, ReferenceGrad, ReferenceValue
# generate mesh
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)
# mark bottom boundary
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
class BottomFace(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.0)
bottom_face = BottomFace()
bottom_face.mark(boundaries, 1) # Mark with '1' for bottom face
# Calculate facet normals, tangent vectors
J = FacetJacobian(mesh)
G1 = J[:,0] # \partial{x_i}/\partial{Xf_1}
G2 = J[:,1] # \partial{x_i}/\partial{Xf_2}
N = cross(G1,G2) # normal vector
N = N/sqrt(dot(N,N)) # normalize N
# n = FacetNormal(mesh)
# simple problem to project them and visualize
param = N # (normal) this is the parameter we want to compute
V = VectorFunctionSpace(mesh, 'DG', 1)
dv = TestFunction(V)
v = TrialFunction(V)
a_V = inner(dv, v)*ds + Constant(0)*inner(dv,v)*dx
L_v = inner(dv, param)*ds
param_out = Function(V)
a_V = assemble(a_V, keep_diagonal=True)
L_v = assemble(L_v)
a_V.ident_zeros()
solve(a_V, param_out.vector(), L_v)
# Print out the values of normal vector at each facet midpoint
# (fenicsproject.disclosure.group)
W = FunctionSpace(mesh, 'CG', 1)
for f in facets(mesh):
if boundaries[f] == 1:
print("Facet index:", f.index())
vertex_coords = []
vertex_indices = []
for vertex in vertices(f):
vertex_coords.append(list(vertex.point().array()))
vertex_indices.append(vertex.index())
facet_dofs = W.dofmap().entity_dofs(mesh, mesh.topology().dim()-1, [f.index()])
vertex_dofs = W.dofmap().entity_dofs(mesh, 0, vertex_indices)
print("Vertex Coords and DOFs:", vertex_coords, facet_dofs, vertex_dofs)
# print parameter values at facet midpoint
mp = f.midpoint()
param_value = N_out(mp)
print("N at midpoint:", param_value)
N_file = File("param.pvd")
N_file << param_out
But the results are weird.
Facet index: 0
Vertex Coords and DOFs: [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0]] [] [1, 5, 6]
N at midpoint: [0.16666667 0. 0.83333333]
Facet index: 3
Vertex Coords and DOFs: [[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]] [] [1, 7, 6]
N at midpoint: [ 0. -0.16666667 -0.83333333]
I think this is coming from the simple problem i am solving for projecting. Is this correct?
Anyways, I am going to use these normals in my weak formulation. I just wanted to check whether the normal and tangent vectors I am calculating are correct.
