# FacetNormal vector components and director vector

**URL:** <https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906>\
**Category:** mesh\
**Created:** [November 30, 2023, 7:35pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906 "2023-11-30T19:35:21Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [November 30, 2023, 7:35pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/1 "2023-11-30T19:35:21Z")

</div>

Hi,

I have a quad mesh over 3 dimensions. I want normal vector for each element. I tried FacetNormal and CellNormal but doesn’t find how to output the normal vector.

```auto
from dolfinx import *
mesh = UnitSquareMesh(4,4) # consider 2D mesh for this case
n = CellNormal(mesh)

```

1. get\_local doesn’t work here to find vector components information.
2. Where is normal vector defined? is it at all guass points or is it average of nodal normal vectors.
3. Since, it is a 2-dimension element. How can I find the initial curvature of mesh element. Is here a way a function in ufl like director\_vector() to find normal and curvatures?

Any help is greatly appreciated!

---

<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 30, 2023, 8:59pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/2 "2023-11-30T20:59:13Z")

</div>

> [@bagla0](#):
>
> `n = CellNormal(mesh)`

Is a symbolic expression, which can be evaluated at any point in the cell. It is not in a function space. You could project it into a function space (for linear grids I would use DG-0).

Usually curvature is better to import from a mesh generator if you are using linear meshes.  
You could have a look at [How to compute curvature of a boundary - #2 by dokken](https://fenicsproject.discourse.group/t/how-to-compute-curvature-of-a-boundary/3193/2)

---

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [December 17, 2024, 11:25pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/3 "2024-12-17T23:25:55Z")

</div>

@dokken

Try1: I want to generate the FacetTangent(mesh), similar to symbolic FacetNormal(mesh). Is there a way I could use inbuilt feature? The purpose is to use facettangent(mesh) in jump and avg for C0-penalty method to ensure C1 continuity. Self created tanger vector functions from [here](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f) doesn’t work.

Is FacetJacobian(mesh) same as FacetTangent(mesh)?

Try2: To trying to match FacetNormal(mesh) result in jump and avg [Biharmonic Example](https://docs.fenicsproject.org/dolfinx/v0.9.0/python/demos/demo_biharmonic.html), I generate facet normal function separately using [facet normal generation](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f). The results are much different.

I therefore don’t understand how FacetNormal(mesh) works? As only FacetNormal(mesh) gives correct answer and my generated function from [facet normal generation](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f) doesn’t.

Can you provide the github code where FacetNormal(mesh) symbolic mathematics is shown ?. The max that I could know is [this code Line 677](https://github.com/FEniCS/ufl/blob/main/ufl/geometry.py)

Thanks a lot!

---

<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:** [December 18, 2024, 6:30am UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/4 "2024-12-18T06:30:00Z")

</div>

> [@bagla0](#):
>
> Try1: I want to generate the FacetTangent(mesh), similar to symbolic FacetNormal(mesh). Is there a way I could use inbuilt feature? The purpose is to use facettangent(mesh) in jump and avg for C0-penalty method to ensure C1 continuity. Self created tanger vector functions from [here](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f) doesn’t work.
> 
> Is FacetJacobian(mesh) same as FacetTangent(mesh)?

In 2D, you can compute the tangent based on the FacetNornal, as

```python
n = FacetNornal(mesh)
t = as_vector((-n[1], n[0]))

```

In 3D it is a bit more involved, as the tangent isn’t unique.

> [@bagla0](#):
>
> The purpose is to use facettangent(mesh) in jump and avg for C0-penalty method to ensure C1 continuity. Self created tanger vector functions from [here](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f) doesn’t work.

Why don’t they work? You would have to be more specific, or show a failing example.

> [@bagla0](#):
>
> I therefore don’t understand how FacetNormal(mesh) works? As only FacetNormal(mesh) gives correct answer and my generated function from [facet normal generation](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f) doesn’t.

Again, not giving actual codes that fail makes it hard to help you.

---

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [January 19, 2025, 2:11am UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/5 "2025-01-19T02:11:24Z")

</div>

> [@dokken](#):
>
> Why don’t they work? You would have to be more specific, or show a failing example

@dokken

Thanks for your response. The belwoMWE show how, the [Biharmonic example solution]

([Biharmonic equation — DOLFINx 0.9.0 documentation](https://docs.fenicsproject.org/dolfinx/v0.9.0/python/demos/demo_biharmonic.html)) is compared with created normal vector from [facet vector approx.](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f). The final results fails to match. Can you provide how fact vector created normal gives same solution as `FacetNormal(mesh)`.

```auto
import importlib.util
from mpi4py import MPI
import basix
import dolfinx
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem, assemble_matrix
from dolfinx.mesh import CellType, GhostMode
from ufl import CellDiameter, FacetNormal, avg, div, dS, dx, grad, inner, jump, pi, sin
import basix.ufl
import numpy as np
import dolfinx as dfx
from petsc4py import PETSc

msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=((0.0, 0.0), (1.0, 1.0)),
    n=(2, 2),
    cell_type=CellType.triangle,
    ghost_mode=GhostMode.shared_facet,
)
V = fem.functionspace(msh, ("Lagrange", 2))

tdim = msh.topology.dim
msh.topology.create_connectivity(tdim - 1, tdim)
facets = mesh.exterior_facet_indices(msh.topology)

dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)

bc = fem.dirichletbc(value=0.0, dofs=dofs, V=V)

alpha = 8.0
h = CellDiameter(msh)
n = FacetNormal(msh)
h_avg = (h("+") + h("-")) / 2.0

###################################################################
###################### Copied Code ####################################

def tangential_projection(u: ufl.Coefficient, n: ufl.FacetNormal) -> ufl.Coefficient:
    return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u

def facet_vector_approximation(V: dfx.fem.FunctionSpace,
                               mt: dfx.mesh.MeshTags | None = None,
                               mt_id: int | None = None,
                               tangent: bool = False,
                               interior: bool = False,
                               jit_options: dict | None = None,
                               form_compiler_options: dict | None = None) -> dfx.fem.Function:
    
    jit_options = jit_options if jit_options is not None else {}
    form_compiler_options = form_compiler_options if form_compiler_options is not None else {}

    comm = V.mesh.comm # MPI Communicator
    n = ufl.FacetNormal(V.mesh) # UFL representation of mesh facet normal
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V) # Trial and test functions

    # Define integral measure
    if interior:
        # Create interior facet integral measure
        dS = ufl.dS(domain=V.mesh) if mt is None else ufl.dS(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
    else:
        # Create exterior facet integral measure
        ds = ufl.ds(domain=V.mesh) if mt is None else ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id) 
    
    # If tangent==True, the right-hand side of the problem should be a tangential projection of the facet normal vector.
    if tangent:
        if V.mesh.geometry.dim == 2:
            if interior:
                a = (ufl.inner(u('+'), v('+')) + ufl.inner(u('-'), v('-'))) * dS
                L = ufl.inner(ufl.as_vector([-n('+')[1], n('+')[0]]), v('+')) * dS \
                  + ufl.inner(ufl.as_vector([-n('-')[1], n('-')[0]]), v('-')) * dS
            else:
                a = ufl.inner(u, v) * ds
                L = ufl.inner(ufl.as_vector([-n[1], n[0]]), v) * ds
        else:
            c = dfx.fem.Constant(V.mesh, (1.0, 1.0, 1.0)) # Vector to tangentially project the facet normal vectors on 
#####################################################################
         
            if interior:
                a = (ufl.inner(u('+'), v('+')) + ufl.inner(u('-'), v('-'))) * dS
                L = ufl.inner(tangential_projection(c, n('+')), v('+')) * dS \
                  + ufl.inner(tangential_projection(c, n('-')), v('-')) * dS
            else:
                a = ufl.inner(u, v) * ds
                L = ufl.inner(tangential_projection(c, n), v) * ds

    # If tangent==false the right-hand side is simply the facet normal vector.
    else:
        if interior:
            a = (ufl.inner(u('+'), v('+')) + ufl.inner(u('-'), v('-'))) * dS
            L = (ufl.inner(n('+'), v('+')) + ufl.inner(n('-'), v('-'))) * dS
        else:
            a = ufl.inner(u, v) * ds
            L = ufl.inner(n, v) * ds

    # Find all boundary dofs, which are the dofs where we want to solve for the facet vector approximation.
    # Start by assembling test functions integrated over the boundary integral measure.
    ones = dfx.fem.Constant(V.mesh, dfx.default_scalar_type((1,) * V.mesh.geometry.dim)) # A vector of ones
    if interior:
        local_val = dfx.fem.form((ufl.dot(ones, v('+')) + ufl.dot(ones, v('-')))*dS)
    else:
        local_val = dfx.fem.form(ufl.dot(ones, ufl.TestFunction(V))*ds) 
    local_vec = dfx.fem.assemble_vector(local_val)

    # For the dofs that do not lie on the boundary of the mesh the assembled vector has value zero.
    # Extract these dofs and use them to deactivate the corresponding block in the linear system we will solve.
    bdry_dofs_zero_val = np.flatnonzero(np.isclose(local_vec.array, 0))
    deac_blocks = np.unique(bdry_dofs_zero_val // V.dofmap.bs).astype(np.int32)

    # Create sparsity pattern by manipulating the blocks to be deactivated and set
    # a zero Dirichlet boundary condition for these dofs.
    bilinear_form = dfx.fem.form(a, jit_options=jit_options,
                                 form_compiler_options=form_compiler_options)
    pattern = dfx.fem.create_sparsity_pattern(bilinear_form)
    pattern.insert_diagonal(deac_blocks)
    pattern.finalize()
    u_0 = dfx.fem.Function(V)
    u_0.vector.set(0)
    bc_deac = dfx.fem.dirichletbc(u_0, deac_blocks)

    # Create the matrix
    A = dfx.cpp.la.petsc.create_matrix(comm, pattern)
    A.zeroEntries()

    # Assemble the matrix with all entries
    form_coeffs = dfx.cpp.fem.pack_coefficients(bilinear_form._cpp_object)
    form_consts = dfx.cpp.fem.pack_constants(bilinear_form._cpp_object)
    dfx.fem.petsc.assemble_matrix(A, bilinear_form, constants=form_consts, coeffs=form_coeffs, bcs=[bc_deac])

    # Insert the diagonal with the deactivated blocks.
    if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
        A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
        A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
        dfx.cpp.fem.petsc.insert_diagonal(A=A, V=bilinear_form.function_spaces[0], bcs=[bc_deac._cpp_object], diagonal=1.0)
    A.assemble()

    # Assemble the linear form and the right-hand side vector.
    linear_form = dfx.fem.form(L, jit_options=jit_options,
                               form_compiler_options=form_compiler_options)
    b = dfx.fem.petsc.assemble_vector(linear_form)

    # Apply lifting to the right-hand side vector and set boundary conditions.
    dfx.fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    dfx.fem.petsc.set_bc(b, [bc_deac])

    # Setup a linear solver using the Conjugate Gradient method.
    solver = PETSc.KSP().create(MPI.COMM_WORLD)
    solver.setType("cg")
    solver.rtol = 1e-8
    solver.setOperators(A)

    # Solve the linear system and perform ghost update.
    nh = dfx.fem.Function(V) # Function for the facet vector approximation
    solver.solve(b, nh.vector)
    nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

    # Normalize the vectors to get the unit facet normal/tangent vector.
    nh_norm = ufl.sqrt(ufl.inner(nh, nh)) # Norm of facet vector
    cond_norm = ufl.conditional(ufl.gt(nh_norm, 1e-10), nh_norm, 1.0) # Avoid division by zero
    
    if V.mesh.geometry.dim == 1:
        nh_norm_vec = ufl.as_vector((nh[0]/cond_norm))
    elif V.mesh.geometry.dim == 2:
        nh_norm_vec = ufl.as_vector((nh[0]/cond_norm, nh[1]/cond_norm))
    elif V.mesh.geometry.dim == 3:
        nh_norm_vec = ufl.as_vector((nh[0]/cond_norm, nh[1]/cond_norm, nh[2]/cond_norm))

    nh_normalized = dfx.fem.Expression(nh_norm_vec, V.element.interpolation_points())
    
    n_out = dfx.fem.Function(V)
    n_out.interpolate(nh_normalized)

    return n_out

DEFAULT=2
SUBSET=3
mesh=msh
mesh.topology.create_entities(mesh.topology.dim-1)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

# Create an array facet_marker which contains the facet marker value for all facets in the mesh.
facet_dim = mesh.topology.dim-1 # Topological dimension of facets
num_facets = mesh.topology.index_map(facet_dim).size_local + mesh.topology.index_map(facet_dim).num_ghosts # Total number of facets in mesh
facet_marker = np.full(num_facets, DEFAULT, dtype=np.int32) # Default facet marker value is 2

# Find the subset of facets to be marked.
boundary_facets = dolfinx.mesh.locate_entities_boundary(
        mesh, 1, lambda x: np.full(x.shape[1], True, dtype=bool)) # Boundary mesh afcets

facets_ids = dolfinx.mesh.locate_entities(
        mesh, 1, lambda x: np.full(x.shape[1], True, dtype=bool)) # total mesh facets

facets_int=np.setdiff1d(facets_ids,boundary_facets ) #Interior mesh facets

# subset_facets = dfx.mesh.locate_entities(mesh, facet_dim, locator) # Get subset of facets to be marked
facet_marker[facets_int] = SUBSET # Fill facet marker array with the value SUBSET
facet_tags = dfx.mesh.meshtags(mesh, facet_dim, np.arange(num_facets, dtype=np.int32), facet_marker) # Create facet meshtags
ft_id = SUBSET # Set the meshtags id for which we will approximate the facet vectors

# Create a DG1 space for the facet vectors to be approximated.
DG1 = basix.ufl.element(family="Lagrange", cell=mesh.basix_cell(), degree=1, discontinuous=True, shape=(mesh.geometry.dim,))
space = dfx.fem.functionspace(mesh=mesh, element=DG1)

# Computed facet vector approximation. 
nn= facet_vector_approximation(V=space, mt=facet_tags, mt_id=ft_id, interior=True, tangent=False)

############COPIED CODE ENDS ###########################
############################################################

# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 4.0 * pi**4 * sin(pi * x[0]) * sin(pi * x[1])

a_FacetNormal = (
    inner(div(grad(u)), div(grad(v))) * dx
    - inner(avg(div(grad(u))), jump(grad(v), n)) * dS
    - inner(jump(grad(u), n), avg(div(grad(v)))) * dS
    + alpha / h_avg * inner(jump(grad(u), n), jump(grad(v), n)) * dS
)
a_Created_FacetNormal = (
    inner(div(grad(u)), div(grad(v))) * dx
    - inner(avg(div(grad(u))), jump(grad(v), nn)) * dS
    - inner(jump(grad(u), nn), avg(div(grad(v)))) * dS
    + alpha / h_avg * inner(jump(grad(u), nn), jump(grad(v), nn)) * dS
)

L = inner(f, v) * dx
problem_1 = LinearProblem(a_FacetNormal, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

uh_1 = problem_1.solve()
print('Biharmonic Sol: Using n=FacetNormal(msh)\n')
print(uh_1.vector[:])

problem_2 = LinearProblem(a_Created_FacetNormal, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

uh_2 = problem_2.solve()
print('\n Biharmonic Sol: Using nn=facet normal approximation \n')
print(uh_2.vector[:])

```

I apologize for length of problem, but majority of code is taken directly from [facet vector approx.](https://gist.github.com/hherlyng/cb3ab37dc58205bcecbc9a6e15b1267f)

---

<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:** [January 19, 2025, 3:47pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/6 "2025-01-19T15:47:22Z")

</div>

Dear @bagla0,

Why would you expect the results to match?  
You are using two different normals, that vastly differs (especially on coarse grids, as indicated by inspecting the solution:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/1/113a4f0d43338b4839ef5f89e7043b3fac2981f2.png)  
As you can see there, the normals are completely different form the actual facets.

This will always be a case with problems where you want to project a normal into a space it is not in.

As stated previously, you can use a symbolic version of the tangent, by calling a tangential projection, i.e. use

```python
u_t = tangential_projection(u, ufl.FacetNormal(mesh))

```

without having to put it in a function space.

---

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [January 20, 2025, 8:59pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/7 "2025-01-20T20:59:26Z")

</div>

Dear @dokken,

> [@dokken](#):
>
> `u_t = tangential_projection(u, ufl.FacetNormal(mesh))`

Using your way, I tried tangential\_projection of `FacetNormal(mesh)` as

```auto
n = FacetNormal(mesh)
P=ufl.Identity(3) - ufl.outer(n, n) # Projection matrix of tangent plane for 3D space.

# Project reference vector in tangential plane, 
# The reference vector is one of the two tangents of quad mesh, as

t=ufl.Jacobian(mesh)
r=as_vector([t[0, 0], t[1, 0], t[2, 0]]) # reference vector (different for each element)
t1=P*r # Projected tangent vector in tangential plane (normal to n)

t2=ufl.cross(n,t1) # Projected second tangent vector

```

Q1: How can I plot n, t1 and t2 in given mesh (like you did for previous case), to verify t1 and t2 are actually correct tangents to facets.  
On interpolating to `Function(mesh)`, I am getting

```auto
V_l = dolfinx.fem.functionspace(msh, basix.ufl.element(
    "DG", msh.topology.cell_name(),1, shape=(3, )))

el=dolfinx.fem.Function(V_l)
    
fexpr1=dolfinx.fem.Expression(n,V_l.element.interpolation_points(), comm=MPI.COMM_WORLD)
el.interpolate(fexpr1)

```

`RuntimeError: Expressions for cells do not support <class 'ufl.geometry.ReferenceNormal'>`.

Q2: I wanted one of the facet tangent t1/t2 to lie along facet/edge.

The purpose of computing t1 and t2, is to weakly constraint the derivative continuity along the orthogonal triad directions. This will make all first derivatives continuous. I apologize for asking this again, but, I think with this approach the weak continuity is completely implemented when t1 and t2 directions are included as:

```auto
 # Extra terms used in C0-IPM method for weak C1 continuity, 
nn= - inner(avg(div(grad(dvl))), jump(grad(v_l), n))*dS \
          - inner(jump(grad(dvl), n), avg(div(grad(v_l))))*dS \
          + (alpha/h_avg**2)*inner(jump(grad(dvl),n), jump(grad(v_l),n))*dS
tt1=- inner(avg(div(grad(dvl))), jump(grad(v_l), t1))*dS \
          - inner(jump(grad(dvl), t1), avg(div(grad(v_l))))*dS \
          + (beta/h_avg**2)*inner(jump(grad(dvl),t1), jump(grad(v_l),t1))*dS
tt2=- inner(avg(div(grad(dvl))), jump(grad(v_l), t2))*dS \
          - inner(jump(grad(dvl), t2), avg(div(grad(v_l))))*dS \
          + (beta/h_avg**2)*inner(jump(grad(dvl),t2), jump(grad(v_l),t2))*dS

```

Thanking You

---

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [January 28, 2025, 11:59pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/8 "2025-01-28T23:59:51Z")

</div>

Just claryfying some confusion regarding direction, as cannot extract the numerical data from that variable. Sorry for the dumb questions.  
For a quadrilateral mesh defined in 3D space,

1. Does n=ufl.FacetNormal(mesh) has a) normal to edge lying in plane of quad element or b) normal to quad element plane (similar direction as CellNormal(mesh)) ?
2. Does t=ufl.classes.FacetJacobian(mesh) denotes the tangent vector along the edge,
3. Both n=ufl.FacetNormal(mesh) and t=ufl.classes.FacetJacobian(mesh) lie in element plane and cross(n,t) is equivalent to CellNormal(mesh)

---

<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:** [January 31, 2025, 5:20pm UTC](https://fenicsproject.discourse.group/t/facetnormal-vector-components-and-director-vector/12906/9 "2025-01-31T17:20:36Z")

</div>

Here is an example where one can visualize the `FacetNormal` properly:

```python
from mpi4py import MPI
import dolfinx
import basix.ufl
import scifem
import ufl

def move_to_facet_quadrature(ufl_expr, mesh, sub_facets, scheme="default", degree=6):
    fdim = mesh.topology.dim - 1
    # Create submesh
    bndry_mesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, fdim, sub_facets)
    # Create quadrature space on submesh
    q_el = basix.ufl.quadrature_element(
        bndry_mesh.basix_cell(), ufl_expr.ufl_shape, scheme, degree
    )
    Q = dolfinx.fem.functionspace(bndry_mesh, q_el)

    # Compute where to evaluate expression per submesh cell
    integration_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.exterior_facet, mesh.topology, entity_map, fdim
    )
    compiled_expr = dolfinx.fem.Expression(ufl_expr, Q.element.interpolation_points())

    # Evaluate expression
    q = dolfinx.fem.Function(Q)
    q.x.array[:] = compiled_expr.eval(mesh, integration_entities).reshape(-1)
    return q

# Facet expression evaluations

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
n = ufl.FacetNormal(mesh)
mesh.topology.create_connectivity(1, 2)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
n_h = move_to_facet_quadrature(n, mesh, exterior_facets, degree=2)

with dolfinx.io.XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
with scifem.xdmf.XDMFFile("approximation.xdmf", [n_h]) as xdmf:
    xdmf.write(0.0)

```

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4e9510265869cf6be15183ddaaf26e18c3d013c9.png)  
Taking the tangential projection one gets  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/738c482a80e55e347019ad26a25ca4fbfbb24fbf.png)  
(i.e. with code)

```python
from mpi4py import MPI
import dolfinx
import basix.ufl
import scifem
import ufl

def move_to_facet_quadrature(ufl_expr, mesh, sub_facets, scheme="default", degree=6):
    fdim = mesh.topology.dim - 1
    # Create submesh
    bndry_mesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, fdim, sub_facets)
    # Create quadrature space on submesh
    q_el = basix.ufl.quadrature_element(
        bndry_mesh.basix_cell(), ufl_expr.ufl_shape, scheme, degree
    )
    Q = dolfinx.fem.functionspace(bndry_mesh, q_el)

    # Compute where to evaluate expression per submesh cell
    integration_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.exterior_facet, mesh.topology, entity_map, fdim
    )
    compiled_expr = dolfinx.fem.Expression(ufl_expr, Q.element.interpolation_points())

    # Evaluate expression
    q = dolfinx.fem.Function(Q)
    q.x.array[:] = compiled_expr.eval(mesh, integration_entities).reshape(-1)
    return q

# Facet expression evaluations

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (2,)))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: (x[0], x[1]))
n = ufl.FacetNormal(mesh)

def tangential_projection(u: ufl.Coefficient, n: ufl.FacetNormal) -> ufl.Coefficient:
    return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u

u_t = tangential_projection(u, n)

mesh.topology.create_connectivity(1, 2)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
n_h = move_to_facet_quadrature(u_t, mesh, exterior_facets, degree=2)

with dolfinx.io.XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
with scifem.xdmf.XDMFFile("approximation.xdmf", [n_h]) as xdmf:
    xdmf.write(0.0)

```