Hi there,
I skimmed over the discourse but could not pick the most proper and easiest way of getting normal vector field projected on a CG-1 space. All I want is;
I checked this, this and this.
Especially this seems the way to go but not sure if it is. Any help would be much appreciated.
Hello,
I managed to get it working; Here is the function that takes function space and returns the facet normal field.
# Copyright (C) 2023 Jørgen S. Dokken and Halvor Herlyng
#
# SPDX-License-Identifier: MIT
import dolfinx as dfx
from petsc4py import PETSc
import ufl
from mpi4py import MPI
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, set_bc
from dolfinx.cpp.fem.petsc import insert_diagonal
def facet_vector_approximation(V: dfx.fem.FunctionSpace,
mt: dfx.mesh.MeshTags = None,
jit_options: dict = None,
form_compiler_options: dict = None) -> dfx.fem.Function:
"""
Approximate the facet normal or tangent on a surface by projecting it into the function space for a set of facets.
Args:
V: The function space to project into. Needs to be defined on discontinuous elements.
mt: The `dolfinx.mesh.MeshTags` containing facet markers.
mt_id: The id for the facets in `mt` we want to represent the normal or tangent on.
tangent: To approximate the tangent of the facets set this flag to `True`.
jit_options: Parameters used in CFFI JIT compilation of C code generated by FFCx.
See https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/jit.py#L22-L37
for all available parameters. Takes priority over all other parameter values.
form_compiler_options: Parameters used in FFCx compilation of this form. Run `ffcx - -help` at
the commandline to see all available options. Takes priority over all
other parameter values, except for `scalar_type` which is determined by
DOLFINx.
"""
# if not V.element.basix_element.discontinuous:
# raise RuntimeError("Discontinuous elements are required for a proper representation of the facet vectors.")
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
nh = dfx.fem.Function(V) # Function for the facet vector approximation
u, v = ufl.TrialFunction(V), ufl.TestFunction(V) # Trial and test functions
#
# Define integral measure to create exterior facet
ds = ufl.ds(domain=V.mesh , subdomain_data=mt)
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
local_val = dfx.fem.form(ufl.inner(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)
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)
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 = assemble_vector(linear_form)
# Apply lifting to the right-hand side vector and set boundary conditions.
apply_lifting(b, [bilinear_form], [[bc_deac]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
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-12
solver.setOperators(A)
# Solve the linear system and perform ghost update.
solver.solve(b, nh.vector)
nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# print(np.unique(nh.vector[:]))
# Normalize the vectors to get the unit facet normal/tangent vector.
nh_norm = ufl.sqrt(ufl.inner(nh, nh)) # Norm of facet vector
if V.mesh.geometry.dim == 1:
nh_norm_vec = ufl.as_vector((nh[0]/nh_norm))
elif V.mesh.geometry.dim == 2:
nh_norm_vec = ufl.as_vector((nh[0]/nh_norm, nh[1]/nh_norm))
elif V.mesh.geometry.dim == 3:
nh_norm_vec = ufl.as_vector((nh[0]/nh_norm, nh[1]/nh_norm, nh[2]/nh_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
Please note that this is a simplified version of @hherlyng’s post: How to compute tangent vectors on facets (2D surface parametrization) - #3 by hherlyng
i.e. the gist:
Facet normal and facet tangent vector approximation in dolfinx · GitHub
Please note that this is distributed under an MIT license, and the license with its copyright holders should be included in your snippet (if you used it)
Thanks, post edited.