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
# SPDX-License-Identifier: MIT

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),
    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)]
        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)
    points = np.zeros((0, 3), dtype=domain.geometry.x.dtype)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
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])

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

ksp = PETSc.KSP().create(domain.comm)

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

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

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)
a = inner(grad(u), grad(v))*dx

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)
    ps = PointSource(V, [])

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