Normal vector in manifold

General
5

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