# Line Integral of Vector Field

Hi everyone!

Please consider a calculation of a vector field \mathbf{u}, for example via the solution of the static Stokes equations with zero Dirichlet Boundary Conditions with Taylor-Hood finite elements.

from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py import PETSc
import basix, basix.ufl_wrapper
import matplotlib.pylab as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore")
from ufl import (VectorElement, FiniteElement,
SpatialCoordinate, TrialFunction, TestFunction,
as_vector, cos, sin, inner, div, grad, dx, pi)

def u_ex(x):
sinx = sin(pi * x[0])
siny = sin(pi * x[1])
cosx = cos(pi * x[0])
cosy = cos(pi * x[1])
c_factor = 2 * pi * sinx * siny
return c_factor * as_vector((cosy * sinx, - cosx * siny))

def p_ex(x):
return sin(2 * pi * x[0]) * sin(2 * pi * x[1])

def source(x):
u, p = u_ex(x), p_ex(x)

def create_bilinear_form(V, Q):
u, p = TrialFunction(V), TrialFunction(Q)
v, q = TestFunction(V), TestFunction(Q)
a_up = inner(p, div(v)) * dx
a_pu = inner(div(u), q) * dx
return fem.form([[a_uu, a_up], [a_pu, None]])

def create_linear_form(V, Q):
v, q = TestFunction(V), TestFunction(Q)
domain = V.mesh
x = SpatialCoordinate(domain)
f = source(x)
return fem.form([inner(f, v) * dx,
inner(fem.Constant(domain, 0.), q) * dx])

def create_velocity_bc(V):
domain = V.mesh
g = fem.Constant(domain, [0., 0.])
tdim = domain.topology.dim
domain.topology.create_connectivity(tdim - 1, tdim)
bdry_facets = mesh.exterior_facet_indices(domain.topology)
dofs = fem.locate_dofs_topological(V, tdim - 1, bdry_facets)
return [fem.dirichletbc(g, dofs, V)]

def create_nullspace(rhs_form):
null_vec = fem.petsc.create_vector_nest(rhs_form)
null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0)
null_vecs[1].set(1.0)
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
return nsp

def create_preconditioner(Q, a, bcs):
p, q = TrialFunction(Q), TestFunction(Q)
a_p11 = fem.form(inner(p, q) * dx)
a_p = fem.form([[a[0][0], None],
[None, a_p11]])
P = fem.petsc.assemble_matrix_nest(a_p, bcs)
P.assemble()
return P

def assemble_system(lhs_form, rhs_form, bcs):
A = fem.petsc.assemble_matrix_nest(lhs_form, bcs=bcs)
A.assemble()

b = fem.petsc.assemble_vector_nest(rhs_form)
fem.petsc.apply_lifting_nest(b, lhs_form, bcs=bcs)
for b_sub in b.getNestSubVecs():
mode=PETSc.ScatterMode.REVERSE)
spaces = fem.extract_function_spaces(rhs_form)
bcs0 = fem.bcs_by_block(spaces, bcs)
fem.petsc.set_bc_nest(b, bcs0)
return A, b

def create_block_solver(A, b, P, comm):
ksp = PETSc.KSP().create(comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")

nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(("u", nested_IS[0][0]),
("p", nested_IS[0][1]))

# Set the preconditioners for each block
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

# Monitor the convergence of the KSP
ksp.setFromOptions()
return ksp

def solve_stokes(u_element, p_element, domain):
V = fem.FunctionSpace(domain, u_element)
Q = fem.FunctionSpace(domain, p_element)

lhs_form = create_bilinear_form(V, Q)
rhs_form = create_linear_form(V, Q)

bcs = create_velocity_bc(V)
nsp = create_nullspace(rhs_form)
A, b = assemble_system(lhs_form, rhs_form, bcs)
assert nsp.test(A)
A.setNullSpace(nsp)

P = create_preconditioner(Q, lhs_form, bcs)
ksp = create_block_solver(A, b, P, domain.comm)

u, p = fem.Function(V), fem.Function(Q)
w = PETSc.Vec().createNest([u.vector, p.vector])
ksp.solve(b, w)
assert ksp.getConvergedReason() > 0
u.x.scatter_forward()
p.x.scatter_forward()
return (u, p)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.triangle)
element_u = VectorElement("CG", "triangle", 3)
element_p = FiniteElement("CG", "triangle", 2)

u, p = solve_stokes(element_u, element_p, domain)


Consider now that I need to know the value of a line integral of the vector field I found:

I=\int_c{d\mathbf{x}\cdot \mathbf{u}}

For example consider the path c to be a straight line from point \mathbf{x_0} to \mathbf{x_1} where those points are points of the mesh that may or may not be in the boundary.

What is the easiest way to do this, that may also apply if I have a custom 3d mesh from PyMesh, for example?

If the lines do not align with the mesh, you would have to set up a custom integral, i.e.
Define the quadrature scheme (points and weights). Assuming that the quadrature points have been pushed to the physical space, then evaluate u at these points, multiply by the quadrature weight an sum the contributions.

Thanks for the suggestion. I could parametrise the line integral, define a proper quadrature scheme and evaluate u at the quadrature points as explained, for example, here Implementation — FEniCSx tutorial and procced accordingly

However, considering that we have a line aligned in the mesh, i.e. in a unit square a line from (0.5,0) to (1,0.5) is there a way to define the integral with the help of ufl.dx?

The line integral would then be a ds = ufl.Measure("ds", subdomain_data=facet_tags, subdomain_id=5 ) where facet_tags is a mesh-tag marking the facets who intersect the line (0.5,0), (1, 0.5).
i.e. something along the lines of

def facet_eq(x):
return np.isclose(x[1], x[0]-0.5)
facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, facet_eq)
facet_tags = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, facets, np.full_like(facets, 5))


Thanks for the help so far. I would like to try the line integral with the help of ufl.measure. The idea is that I may need to add this integral as a constrain (with a lagrange multiplier) in a FEniCSx calculation in the next part of the model I am trying to make.

So, consider a simple example. Lets say I want to calculate the line integral of the scalar field u(x,y)=1+x^2+2y^2 in a somewhat trivial path from (x=0.5,y=0) to (x=0.5,y=1). The minimal code is as follows

from dolfinx import fem, mesh
from mpi4py import MPI
import numpy as np
import ufl
import dolfinx

domain = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, mesh.CellType.triangle)
mesh=domain

V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
u_test = dolfinx.fem.Function(V)
u_test.interpolate(lambda x: 1 + x[0]**2 + 2*x[1]**2)

# Points of the domain
points = domain.geometry.x

def facet_eq(x):
return np.isclose(x[0], 0.5)
facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, facet_eq)
facet_tags = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, facets, np.full_like(facets, 5))

# Make sure you are in the right line by printing the points
for i in range(len(points)):
a=facet_eq(points[i])
if a == True : print(points[i])

# Define the measure
ds = ufl.Measure("ds", subdomain_data=facet_tags, subdomain_id=5 )
#Calculate the integral
intergal = dolfinx.fem.form(u_test * ufl.ds)
print(fem.assemble_scalar(intergal))


Running the code outputs approximately 9 that is not the actual result, that is approximately 1.92.

Please , if you could provide more feedback I would appreciate it.

You are not using the ds you defined in your integral, so it is currently integrating over all exterior facets.

intergal = dolfinx.fem.form(u_test * ds)


would be correct

Thank you very much!!
Turns out there was one more mistake, the right _integral_types command would be dS with capital S, I shouldn’t have missed that one. So, I post the complete code in case somebody needs a similar line integral calculation. By the way the integral is calculated very accurately.

Thanks again!

from dolfinx import fem, mesh
from mpi4py import MPI
import numpy as np
import ufl
import dolfinx

domain = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, mesh.CellType.triangle)
mesh=domain

V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
u_test = dolfinx.fem.Function(V)
u_test.interpolate(lambda x: 1 + x[0]**2 + 2*x[1]**2)

# Points of the domain
points = domain.geometry.x

def facet_eq(x):
return np.isclose(x[0], 0.5)# & np.isclose(x[1], 0.5))
facets = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, facet_eq)
facet_tags = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, facets, np.full_like(facets, 5))

# Make sure you are in the right line
for i in range(len(points)):
a=facet_eq(points[i])
if a == True : print(points[i])

# Define the measure
ds = ufl.Measure("dS", subdomain_data=facet_tags, subdomain_id=5 )

intergal = dolfinx.fem.form(u_test * ds)
print(fem.assemble_scalar(intergal))


Hopefully a final followup question. Seems that the 3d version of the previous code is not as easy as I thought. Let me consider the line integral of the scalar function u=1+x^2+2y^2+3z^2 across a trivial straight line starting from (0.5,0.5,0) to (0.5,0.5,1). I changed the dimensions accordingly and also the _integral_types to dx but if I print the facets seems that no facet is found, therfore obviously the integral is zero, which is not correct. Minimal code follows:

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl

domain = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10, dolfinx.mesh.CellType.tetrahedron)

V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
u_test = dolfinx.fem.Function(V)
u_test.interpolate(lambda x: 1 + x[0]**2 + 2*x[1]**2 + 3*x[2]**2)

def facet_eq(x):
return np.logical_and(np.isclose(x[0], 0,5), np.isclose(x[1], 0,5))

facets = dolfinx.mesh.locate_entities(domain, domain.topology.dim, facet_eq)
facet_tags = dolfinx.mesh.meshtags(domain, domain.topology.dim, facets, np.full_like(facets, 5))

# Define the measure
ds = ufl.Measure("dx", subdomain_data=facet_tags, subdomain_id=5)

intergal = dolfinx.fem.form(u_test * ds)
print(facets,dolfinx.fem.assemble_scalar(intergal))


A line integral in 3D is not an integral over facets, but edges. This integration measure has not been fully implemented in Dolfinx (one would have to either use @jpdean’s Mixed dimensional branch).
or interpolate from then parent mesh to a sub mesh of the edges in question.