Normal vector in manifold

Hello everyone,

I have created a 2D curve mesh (1D topological dimension) like this:

I need to use and visualize the normal unit vector, I looked at some post but in these they use a 3D mesh, is there a way to compute the normal unit vector to use in a variational formulation,

Also, here is my code for the mesh generation

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, create_unit_square, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner

# Save all logging to file
log.set_output_file("log.txt")
# -

# Next, various model parameters are defined:

dt = 5.0e-05  # time step

# A unit square mesh with 96 cells edges in each direction is created,
# and on this mesh a
# {py:class}`FunctionSpaceBase <dolfinx.fem.FunctionSpaceBase>` `ME` is built
# using a pair of linear Lagrange elements.

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    Args:
        model: Gmsh model to add the mesh to.
        name: Name (identifier) of the mesh to add.

    Returns:
        Gmsh model with a circle mesh added.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.2)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.

    Args:
        comm: MPI communicator top create the mesh on.
        model: Gmsh model.
        name: Name (identifier) of the mesh to add.
        filename: XDMF filename.
        mode: XDMF file mode. "w" (write) or "a" (append).

    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

thank you

Consider the following MWE (given your mesh has been read in)

import ufl
n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n
import dolfinx.fem

V = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))
u_expr = dolfinx.fem.Expression(n_oriented, V.element.interpolation_points())
u = dolfinx.fem.Function(V)
u.interpolate(u_expr)
with dolfinx.io.XDMFFile(mesh.comm, "normal.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(u)

Thus you can use sign*n in your variational formulation.

1 Like

Thank you, Dokken. It seems to work. However, I am working with an evolving mesh over time and exporting the functions of the normal vec (normal_vec) tor and a newly defined mean curvature (H) function. It appears these are not updating over time.

Here is a MWE of what I am doing:

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, create_unit_square, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div


# Save all logging to file
log.set_output_file("log.txt")
# -


# Next, various model parameters are defined:

dt = 5.0e-04  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    Args:
        model: Gmsh model to add the mesh to.
        name: Name (identifier) of the mesh to add.

    Returns:
        Gmsh model with a circle mesh added.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.2)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.

    Args:
        comm: MPI communicator top create the mesh on.
        model: Gmsh model.
        name: Name (identifier) of the mesh to add.
        filename: XDMF filename.
        mode: XDMF file mode. "w" (write) or "a" (append).

    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 1, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1],gdim=2))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))
V2 = functionspace(mesh, ("Lagrange", 1))

normal_expr = Expression(n_oriented, V1.element.interpolation_points())
normal_vec = Function(V1)


H_expr = Expression(div(sign*n), V2.element.interpolation_points())
H = Function(V2)

t = 0.0

file1 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu1.xdmf", "w")
file2 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu2.xdmf", "w")
file3 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/normal.xdmf", "w")
file4 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/H.xdmf", "w")


#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
    T = 3 * dt
else:
    T = 10000 * dt

x_initial = mesh.geometry.x[:,0]
y_initial = mesh.geometry.x[:,1]

sigma = 0
rho = 0

sigma = 1e-5
rho = 0.8e-5

z = 0
inc = 0 
impresion = 200

while (t < T):
    t += dt
    inc += 1

    normal_vec.interpolate(normal_expr)
    H.interpolate(H_expr)

    num=inc/impresion-1 

    velx = 1 + sigma*math.sqrt(t)
    vely = 1 + rho*math.sqrt(t)

    mesh.geometry.x[:,0] = x_initial*velx
    mesh.geometry.x[:,1] = y_initial*vely


    if (int(num)==z) or (int(num-z)==0):
        #us, vs = ufl.split(u)

        # Save solution to file (VTK)

        name = "mesh_at_t"+str(t)
        mesh.name = name

        #file1.write_mesh(mesh)
        #file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

        file3.write_mesh(mesh)
        file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

        file4.write_mesh(mesh)
        file4.write_function(H, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
        z += 1
    

file1.close()
file2.close()
file3.close()
file4.close()

print(H.x.array)

I intend to use the normal_vec and H values for a variational formulation but wanted to confirm first if it’s working properly.

Thank you.

Note that when you use linear elements (straight elements), the divergence of the normal is 0 on every element (since the normal is constant on every line segment.

If you run your code for a while, you will see that the normal is updated, but H will always remain 0, as the elements are not curved.
If you use a second order mesh, H will be non-zero.

import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot, mesh
from dolfinx.fem import Function, functionspace, Expression
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, create_unit_square, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div


# Save all logging to file
log.set_output_file("log.txt")
# -


# Next, various model parameters are defined:

dt = 5.0e-04  # time step

try:
    import gmsh  # type: ignore
except ImportError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
    """Create a Gmsh model of a circle.

    Args:
        model: Gmsh model to add the mesh to.
        name: Name (identifier) of the mesh to add.

    Returns:
        Gmsh model with a circle mesh added.

    """
    model.add(name)
    model.setCurrent(name)
    
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.2)

    circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

    # Synchronize OpenCascade representation with gmsh model
    model.occ.synchronize()

    # Add physical marker for cells. It is important to call this
    # function after OpenCascade synchronization
    model.add_physical_group(dim=1, tags=[circle])

    # Generate the mesh
    model.mesh.generate(dim=1)
    model.mesh.setOrder(2)

    return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
    """Create a DOLFINx from a Gmsh model and output to file.

    Args:
        comm: MPI communicator top create the mesh on.
        model: Gmsh model.
        name: Name (identifier) of the mesh to add.
        filename: XDMF filename.
        mode: XDMF file mode. "w" (write) or "a" (append).

    """
    msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
    msh.name = name
    ct.name = f"{msh.name}_cells"
    ft.name = f"{msh.name}_facets"
    with XDMFFile(msh.comm, filename, mode) as file:
        msh.topology.create_connectivity(0, 1)
        file.write_mesh(msh)
        file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
        file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="circle")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 1, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1],gdim=2))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim,)))
V2 = functionspace(mesh, ("Lagrange", 2))

normal_expr = Expression(n_oriented, V1.element.interpolation_points())
normal_vec = Function(V1)


H_expr = Expression(div(sign*n), V2.element.interpolation_points())
H = Function(V2)

t = 0.0

file1 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu1.xdmf", "w")
file2 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu2.xdmf", "w")
file3 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/normal.xdmf", "w")
file4 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/H.xdmf", "w")


#  Reduce run time if on test (CI) server
T = 10000 * dt

x_initial = mesh.geometry.x[:,0]
y_initial = mesh.geometry.x[:,1]

sigma = 0
rho = 0

sigma = 1e-5
rho = 0.8e-5

z = 0
inc = 0 
impresion = 200

while (t < T):
    t += dt
    inc += 1

    normal_vec.interpolate(normal_expr)
    H.interpolate(H_expr)

    num=inc/impresion-1 

    velx = 1 + sigma*math.sqrt(t)
    vely = 1 + rho*math.sqrt(t)

    mesh.geometry.x[:,0] = x_initial
    mesh.geometry.x[:,1] = y_initial*vely


    if (int(num)==z) or (int(num-z)==0):
        #us, vs = ufl.split(u)

        # Save solution to file (VTK)

        name = "mesh_at_t"+str(t)
        mesh.name = name

        #file1.write_mesh(mesh)
        #file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

        file3.write_mesh(mesh)
        file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

        file4.write_mesh(mesh)
        file4.write_function(H, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
        z += 1
    

file1.close()
file2.close()
file3.close()
file4.close()

print(H.x.array)

2 Likes

This worked perfectly, thank you Dokken.