FacetNormal vector components and director vector

[bagla0]

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.

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!

[dokken]

@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

[bagla0]

@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 doesn’t work.

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

Try2: To trying to match FacetNormal(mesh) result in jump and avg Biharmonic Example, I generate facet normal function separately using facet normal generation. 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 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

Thanks a lot!

[dokken]

@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 doesn’t work.

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

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

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 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 doesn’t.

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

[bagla0]

@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) is compared with created normal vector from facet vector approx. . The final results fails to match. Can you provide how fact vector created normal gives same solution as FacetNormal(mesh).

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.

[dokken]

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:

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

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

without having to put it in a function space.

[bagla0]

Dear @dokken,

@dokken

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

Using your way, I tried tangential_projection of FacetNormal(mesh) as

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

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:

 # 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

[bagla0]

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)

[dokken]

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

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)

Taking the tangential projection one gets

(i.e. with code)

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)