# Point Sources (Redux)

Following my post above, with a minor adaptations, a point source can be implemented as follows:

# Create a point source for Poisson problem
# Author: Jørgen S. Dokken

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx
import dolfinx.fem.petsc
import numpy as np
import ufl

def compute_cell_contributions(V, points):
# Determine what process owns a point and what cells it lies within
mesh = V.mesh
_, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
mesh._cpp_object, points, 1e-6)
owning_points = np.asarray(owning_points).reshape(-1, 3)

# Pull owning points back to reference cell
mesh_nodes = mesh.geometry.x
cmap = mesh.geometry.cmaps[0]
ref_x = np.zeros((len(cells), mesh.geometry.dim),
dtype=mesh.geometry.x.dtype)
for i, (point, cell) in enumerate(zip(owning_points, cells)):
geom_dofs = mesh.geometry.dofmap[cell]
ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

# Create expression evaluating a trial function (i.e. just the basis function)
u = ufl.TrialFunction(V)
num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
if len(cells) > 0:
# NOTE: Expression lives on only this communicator rank
expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))

# Strip out basis function values per cell
basis_values = values[:num_dofs:num_dofs*len(cells)]
else:
basis_values = np.zeros(
(0, num_dofs), dtype=dolfinx.default_scalar_type)
return cells, basis_values

N = 80
domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)
domain.name = "mesh"
domain.topology.create_connectivity(1, 2)

V = dolfinx.fem.functionspace(domain, ("Lagrange", 1))

facets = dolfinx.mesh.exterior_facet_indices(domain.topology)

if domain.comm.rank == 0:
points = np.array([[0.68, 0.36, 0]], dtype=domain.geometry.x.dtype)
else:
points = np.zeros((0, 3), dtype=domain.geometry.x.dtype)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a_compiled = dolfinx.fem.form(a)

dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets)
u_bc = dolfinx.fem.Constant(domain, 0.)
bc = dolfinx.fem.dirichletbc(u_bc, dofs, V)

b = dolfinx.fem.Function(V)
b.x.array[:] = 0
cells, basis_values = compute_cell_contributions(V, points)
for cell, basis_value in zip(cells, basis_values):
dofs = V.dofmap.cell_dofs(cell)
b.x.array[dofs] += basis_value
dolfinx.fem.petsc.apply_lifting(b.vector, [a_compiled], [[bc]])
dolfinx.fem.petsc.set_bc(b.vector, [bc])
b.x.scatter_forward()

A = dolfinx.fem.petsc.assemble_matrix(a_compiled, bcs=[bc])
A.assemble()

ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.PREONLY)
ksp.getPC().setType(PETSc.PC.Type.LU)
ksp.getPC().setFactorSolverType("mumps")

uh = dolfinx.fem.Function(V)
ksp.solve(b.vector, uh.vector)
uh.x.scatter_forward()

with dolfinx.io.VTXWriter(domain.comm, "uh.bp", [uh], engine="BP4") as bp:
bp.write(0.0)


Reference solution made in legacy fenics:

from dolfin import *
N = 80
mesh = UnitSquareMesh(N, N)

V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)

b = assemble(Constant(0.0)*v*dx)
point = Point(0.68, 0.36)
if mesh.mpi_comm().rank == 0:
ps = PointSource(V, point, 1.0)
else:
ps = PointSource(V, [])
ps.apply(b)

bc = DirichletBC(V, Constant(0.0), "on_boundary")
A = assemble(a)

[bc.apply(A, b) for bc in [bc]]

uh = Function(V)
solve(A, uh.vector(), b)

File("uh.pvd") << uh


7 Likes