# Point Sources (Redux)

As a prelude to something more complex, I’d like to solve

-\nabla^2 u = \delta(\mathbf{x}- \mathbf{x}_0) \tag{1}

on the unit circle. Given a finite-dimensional space of test functions V_h = \mathrm{span}(\{\phi_i| i = 1,...,n\}), it seems like the most accurate way to solve (1) would be to take its weak-form of and write the rhs as follows -

\int_{\Omega}\nabla u\cdot\nabla \phi_i dx = \int_{\Omega}\delta \phi_idx \\ = \phi_i(\mathbf{x}_0).

From what I’ve seen, older versions of dolfin had a built-in PointSource’ function that performed just this (or something exactly analogous). Has any such functionality been added to dolfinx or is my best bet simply to use a Gaussian approximation? I worry that in using the latter, numerical quadrature of \int_{\Omega} \delta \phi_i\,dx over the appropriate cell may not yield something sufficiently-close to \phi_i(x).

If my only hope is a Gaussian approximation, is there any standard way to choose the width parameter (commonly referred to as \sigma)?

1 Like

Have you used the search button? This question was discused serveral times, and as of December 2023 there is a solution. I’ll leave you to find it, but once you do please post a reply with the link.

For anyone who may have a similar question, a solution was given here - Delta Implementation - by dokken (post 4). Given a rhs function \delta(\mathbf{x} - \mathbf{x}_0), the solution provided allows one to evaluate the basis functions \phi_i st \mathbf{x}_0 \in \mathrm{supp}(\phi_i).

As an addendum, would anyone mind pointing me towards some method for inserting values into a rhs vector corresponding to specific mesh cells?

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

Copied and pasted program and get this error
AttributeError: ‘dolfinx.cpp.mesh.Geometry_float64’ object has no attribute ‘cmaps’. Did you mean: ‘cmap’?
running dolfinx in fenicsx-env on a Mac.
DOLFINx version: 0.8.0 based on GIT commit: of GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment

Hope you had a great meeting this year!

The api slightly changed. Try to change it ss the error message suggests.

For dolfinx 0.8.0, simply use cmap = mesh.geometry.cmap.

For dolfinx main branch (future 0.9.0), a few further changes in the call to dolfinx.cpp.geometry.determine_point_ownership are required due to Make standalone function for interpolation over non-matching meshes (… · FEniCS/dolfinx@ab79530 · GitHub. I copy the updated code below

# 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
point_ownership_data = dolfinx.cpp.geometry.determine_point_ownership(
mesh._cpp_object, points, 1e-6)
owning_points = np.asarray(point_ownership_data.dest_points).reshape(-1, 3)
cells = point_ownership_data.dest_cells

# Pull owning points back to reference cell
mesh_nodes = mesh.geometry.x
cmap = mesh.geometry.cmap
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)
`
2 Likes

mesh.geometry.cmap worked! Thank you!
I am still making sure I understand how code works, so I may be back later!