# Tangent vectors at facet centers

**URL:** <https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816>\
**Category:** General\
**Created:** [November 21, 2023, 4:27pm UTC](https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816 "2023-11-21T16:27:48Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![deshik](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@deshik](https://fenicsproject.discourse.group/u/deshik)\
**Post date:** [November 21, 2023, 4:27pm UTC](https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816/1 "2023-11-21T16:27:48Z")

</div>

Hi,  
I am trying to calculate tangent vectors at cell facet centers using dolfin.  
Tangent vectors,

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/2d57729351c0577ac836fbef77b347bfb074d622.png)

I computed this using FacetJacobian (J) that is available in UFL package.

Now, I would like to compute and print them at facet centers that are on a particular boundary.  
Here is my minimal code,

```auto
from fenics import *
from ufl import Jacobian, replace
from ufl.classes import FacetJacobian, ReferenceGrad, ReferenceValue

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

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]
G2 = J[:,1]
N = cross(G1,G2) # normal vector
N = N/sqrt(dot(N,N)) # normalize N

# print G1, G2, N at each facet center
for f in facets(mesh):
    if boundaries[f] == 1:
    ..............

```

Can someone help me how to print G1, G2 on facet centers? Or am I understanding the G1 and G2 incorrectly??

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 21, 2023, 4:45pm UTC](https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816/2 "2023-11-21T16:45:21Z")

</div>

`G1` and `G2` are symbolic expressions, that can be used in variational forms. They are not easily printed, as it requires code generation to get expressions for them. What you can do is a projection of each of them to a function space (say DG-1), and visualize them in Paraview, similar to: [How to plot normal unit vector of faces in a 2D mesh? - #4 by dokken](https://fenicsproject.discourse.group/t/how-to-plot-normal-unit-vector-of-faces-in-a-2d-mesh/3912/4)

---

<div class="post-metadata">

**Author:** ![deshik](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@deshik](https://fenicsproject.discourse.group/u/deshik)\
**Post date:** [November 21, 2023, 9:19pm UTC](https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816/3 "2023-11-21T21:19:17Z")

</div>

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.

```auto
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.

```auto
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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 21, 2023, 9:44pm UTC](https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816/4 "2023-11-21T21:44:54Z")

</div>

> [@deshik](#):
>
> `param_value = N_out(mp)`

What is N\_out? I cant find its definition in the code you supplied

---

<div class="post-metadata">

**Author:** ![deshik](https://avatars.discourse-cdn.com/v4/letter/d/3e96dc/32.png) [@deshik](https://fenicsproject.discourse.group/u/deshik)\
**Post date:** [November 21, 2023, 9:53pm UTC](https://fenicsproject.discourse.group/t/tangent-vectors-at-facet-centers/12816/5 "2023-11-21T21:53:54Z")

</div>

Sorry!!  
It is a typo. Here is the correct code.

```auto
from fenics import *
from ufl import Jacobian, replace
from ufl.classes import FacetJacobian, ReferenceGrad, ReferenceValue

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)
W = FunctionSpace(mesh, 'CG', 1)

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]
G2 = J[:,1]
N = cross(G1,G2) # normal vector
N = N/sqrt(dot(N,N)) # normalize N
n = FacetNormal(mesh) 

param = N # 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
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)

        mp = f.midpoint()
        print("Facet midpoint:", mp.array())
        param_value = param_out(mp)
        print("param at midpoint:", param_value)    

N_file = File("param.pvd")
N_file << param_out

```

I did check by making param = n, where n is FacetNormal function. Even, this even gives me wrong results.