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)
return - div(grad(u)) + grad(p)
def create_bilinear_form(V, Q):
u, p = TrialFunction(V), TrialFunction(Q)
v, q = TestFunction(V), TestFunction(Q)
a_uu = inner(grad(u), grad(v)) * dx
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():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD,
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")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
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?
Thanks in advance!!!