Point Sources (Redux)

General
1

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)?

2

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.

4

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).

5

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

6

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),
                     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 = 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]])
b.x.scatter_reverse(dolfinx.la.InsertMode.add)
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)

image
image625×602 21.1 KB

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)
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

image
image598×586 20.9 KB

7

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!

8

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

9

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
# 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
    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 = 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]])
b.x.scatter_reverse(dolfinx.la.InsertMode.add)
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)
10

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

11

Is it possible to adapt this solution for Helmholtz equation? I just executed it with complex fenics (0.8, conda) and got the following error:

{
	"name": "TypeError",
	"message": "apply_lifting(): incompatible function arguments. The following argument types are supported:
    1. apply_lifting(b: ndarray[dtype=float32, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_float32], constants: list[ndarray[dtype=float32, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=float32, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_float32]], x0: list[ndarray[dtype=float32, writable=False, shape=(*), order='C']], scale: float) -> None
    2. apply_lifting(b: ndarray[dtype=float64, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_float64], constants: list[ndarray[dtype=float64, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=float64, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_float64]], x0: list[ndarray[dtype=float64, writable=False, shape=(*), order='C']], scale: float) -> None
    3. apply_lifting(b: ndarray[dtype=complex64, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_complex64], constants: list[ndarray[dtype=complex64, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=complex64, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_complex64]], x0: list[ndarray[dtype=complex64, writable=False, shape=(*), order='C']], scale: complex) -> None
    4. apply_lifting(b: ndarray[dtype=complex128, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_complex128], constants: list[ndarray[dtype=complex128, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=complex128, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_complex128]], x0: list[ndarray[dtype=complex128, writable=False, shape=(*), order='C']], scale: complex) -> None

Invoked with types: ndarray, list, list, list, list, list, float",
	"stack": "---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[1], line 77
     75     dofs = V.dofmap.cell_dofs(cell)
     76     b.x.array[dofs] += basis_value
---> 77 dolfinx.fem.petsc.apply_lifting(b.vector, [a_compiled], [[bc]])
     78 b.x.scatter_reverse(dolfinx.la.InsertMode.add)
     79 dolfinx.fem.petsc.set_bc(b.vector, [bc])

File ~/miniconda3/envs/fenics-complex/lib/python3.12/site-packages/dolfinx/fem/petsc.py:651, in apply_lifting(b, a, bcs, x0, scale, constants, coeffs)
    649 x0_r = [x.array_r for x in x0]
    650 b_local = stack.enter_context(b.localForm())
--> 651 _assemble.apply_lifting(b_local.array_w, a, bcs, x0_r, scale, constants, coeffs)

File ~/miniconda3/envs/fenics-complex/lib/python3.12/site-packages/dolfinx/fem/assemble.py:348, in apply_lifting(b, a, bcs, x0, scale, constants, coeffs)
    346 _a = [None if form is None else form._cpp_object for form in a]
    347 _bcs = [[bc._cpp_object for bc in bcs0] for bcs0 in bcs]
--> 348 _cpp.fem.apply_lifting(b, _a, constants, coeffs, _bcs, x0, scale)

TypeError: apply_lifting(): incompatible function arguments. The following argument types are supported:
    1. apply_lifting(b: ndarray[dtype=float32, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_float32], constants: list[ndarray[dtype=float32, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=float32, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_float32]], x0: list[ndarray[dtype=float32, writable=False, shape=(*), order='C']], scale: float) -> None
    2. apply_lifting(b: ndarray[dtype=float64, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_float64], constants: list[ndarray[dtype=float64, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=float64, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_float64]], x0: list[ndarray[dtype=float64, writable=False, shape=(*), order='C']], scale: float) -> None
    3. apply_lifting(b: ndarray[dtype=complex64, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_complex64], constants: list[ndarray[dtype=complex64, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=complex64, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_complex64]], x0: list[ndarray[dtype=complex64, writable=False, shape=(*), order='C']], scale: complex) -> None
    4. apply_lifting(b: ndarray[dtype=complex128, shape=(*), order='C'], a: list[dolfinx.cpp.fem.Form_complex128], constants: list[ndarray[dtype=complex128, writable=False, shape=(*), order='C']], coeffs: list[dict[tuple[dolfinx.cpp.fem.IntegralType, int], ndarray[dtype=complex128, writable=False, shape=(*, *), order='C']]], bcs1: list[list[dolfinx.cpp.fem.DirichletBC_complex128]], x0: list[ndarray[dtype=complex128, writable=False, shape=(*), order='C']], scale: complex) -> None

Invoked with types: ndarray, list, list, list, list, list, float"
}
12

Simply change

to

u_bc = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(0.0))
13

Hello! I try to run your code in Version0.9 but there are still some errors.

Traceback (most recent call last):
 in <module>
    dolfinx.fem.petsc.apply_lifting(b, [a_compiled], [[bc]])
  in apply_lifting
    b_local = stack.enter_context(b.localForm())
AttributeError: 'Function' object has no attribute 'localForm'
14

The code can be updated with:

# 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
    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 = 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.x.petsc_vec, [a_compiled], [[bc]])
b.x.scatter_reverse(dolfinx.la.InsertMode.add)
dolfinx.fem.petsc.set_bc(b.x.petsc_vec, [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.x.petsc_vec, uh.x.petsc_vec)
uh.x.scatter_forward()

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

however, please note that a tested and improved version of this code is available in scifem: Point sources in DOLFINx — scifem

15

Thanks! The tutorial is fantastic!