Reimplement PointSource

[chrisrichardson]

PointSource has been removed and needs to be reimplemented.

[IgorBaratta]

Should we use an interface similar to the older one?

eg:

class PointSource
 {
 public:
   PointSource(std::shared_ptr<const function::FunctionSpace> V,
               const std::vector<std::pair<geometry::Point, PetsScalar>> sources);

   PointSource(std::shared_ptr<const function::FunctionSpace> V0,
               std::shared_ptr<const function::FunctionSpace> V1,
               const std::vector<std::pair<geometry::Point, PetsScalar>> sources);

   ~PointSource() = default;

   void apply(la::PETScVector& b);

   void apply(la::PETScMatrix& A);
   ....
   ....

[chrisrichardson]

It's not clear what is the best approach. Maybe it is better to do it at the form compiler level?

[punyidea]

Has there been any progress on this issue? An implementation would be appreciated from my side.

[chrisrichardson]

@mscroggs - are you working on it?

[mscroggs]

yes, slowly

[jorgensd]

It should in theory be easier to implement this now that we have https://github.com/FEniCS/dolfinx/blob/90dc151f46eba4bc5373a4e91e302ff6e237e1dd/cpp/dolfinx/geometry/utils.h#L131-L148 determine_point_ownership. The strategy would be:

1. Send in points (no duplicates over processes)
2. Determine owner of point
3. Tabulate basis functions or gradients of basis functions at these points.
4. Push basis functions forward to physical element and compute the local element tensor
5. Add them to matrix-vector

[ryansynk]

Just bumping this to say it would be a helpful feature on my end, any progress?

[jorgensd]

I've made a simple way of getting all the info you need for a point source:

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

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl

def compute_cell_contributions(V, points):
    # Determine what process owns a point and what cells it lies within
    _, _, 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))

        # 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

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

if mesh.comm.rank == 0:
    points = np.array([[0.1, 0.2, 0],
                      [0.91, 0.93, 0]], dtype=mesh.geometry.x.dtype)
else:
    points = np.zeros((0, 3), dtype=mesh.geometry.x.dtype)

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

cells, basis_values = compute_cell_contributions(V, points)
print(cells, basis_values)

ref: https://fenicsproject.discourse.group/t/how-to-construct-a-sparse-matrix-used-for-calculating-function-values-at-multiple-points/13081/4

[jorgensd]

A full reimplementation for Python with comparison to legacy dolfin is posted at: https://fenicsproject.discourse.group/t/point-sources-redux/13496/6?u=dokken

[jorgensd]

Do we want to put the re-implementation of point-source into the library, or should we close this issue?