How to construct a sparse matrix used for calculating function values at multiple points

dolfinx
,
1

Hi everyone. Recently, I try to convert my program based on FEniCS 2019 to the new version FEniCSx. I meet a problem for transforming the following code to the current version due to lack of some functions, e.g., dolfin_element.evaluate_basis_all

import fenics as fe

def construct_measurement_matrix(xs, V):
    # This function generate measurement matrix 
    # xs: measurement points
    # V:  function space generated by FEniCS
    # Example: 
    # Let V be a function space generated by FEniCS
    # u is a function genrated by FEniCS based on function space V
    # points = np.array([(0,0), (0.5,0.5)])
    # S = construct_measurement_matrix(ponits, V)
    # S@u.vector()[:] is a vector (u[0, 0], u[0.5, 0.5])

    nx, dim = xs.shape
    mesh = V.mesh()
    coords = mesh.coordinates()
    cells = mesh.cells()
    dolfin_element = V.dolfin_element()
    dofmap = V.dofmap()
    bbt = mesh.bounding_box_tree()
    sdim = dolfin_element.space_dimension()
    v = np.zeros(sdim)
    rows = np.zeros(nx*sdim, dtype='int')
    cols = np.zeros(nx*sdim, dtype='int')
    vals = np.zeros(nx*sdim)
    for k in range(nx):
        # Loop over all interpolation points
        x = xs[k, :]
        if dim > 1:
            p = fe.Point(x[0], x[1])
        elif dim == 1:
            p = fe.Point(x)
        # Find cell for the point
        cell_id = bbt.compute_first_entity_collision(p)
        # Vertex coordinates for the cell
        xvert = coords[cells[cell_id, :], :]
        # Evaluate the basis functions for the cell at x
        v = dolfin_element.evaluate_basis_all(x, xvert, cell_id)
        jj = np.arange(sdim*k, sdim*(k+1))
        rows[jj] = k
        # Find the dofs for the cell
        cols[jj] = dofmap.cell_dofs(cell_id)
        vals[jj] = v

    ij = np.concatenate((np.array([rows]), np.array([cols])), axis=0)
    M = sps.csr_matrix((vals, ij), shape=(nx, V.dim()))
    return M

Any help would be great to have!

Junxiong Jia

2

due to lack of some functions, e.g., dolfin_element.evaluate_basis_all

You may want to have a look at Creating and tabulating an element — Basix 0.7.0 documentation

You may also want to show us the parts of the code you already did manage to convert from FEniCS to FEniCSx. You can expect the community to help, but not to do the entire rewrite of the snippet for you.

3

Sorry for not posting my current parts of the code. I am not an expert for programming finite element methods, just use the functionals of FEniCSx to solve some inverse problems. FEniCSx provides a quick and convenient tool to test Bayesian or optimization inverse algorithms without too much attention on solving PDEs. The following code is my temporary solution, but it must be not efficient enough.

import numpy as np
import scipy
import scipy.sparse as sp
from dolfinx import fem, la, geometry
from dolfinx.fem import assemble_matrix, assemble_vector
import ufl
from mpi4py import MPI


def construct_measure_matrix(points, V):
    '''
    Input:
        V: function space produced by dolfinx.fem.FunctionSpace
        points: input points, e.g.,
                points = [[x1_1,x2_1], [x1_2,x2_2], [x1_3,x2_3]] stands for
                there are two measured points with coordinates:
                x1=(x1_1,x1_2,_x1_3), x2=(x2_1,x2_2,x2_3)
    Output:
        S: a scipy.sparse.csr_matrix
    '''
    assert points.shape[0] == 3, "The points should be shape 3xN"
    domain = V.mesh
    bb_tree = geometry.bb_tree(domain, domain.topology.dim)
    cells = []
    points_on_proc = []
    cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
    colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
    for i, point in enumerate(points.T):
        if len(colliding_cells.links(i)) > 0:
            points_on_proc.append(point)
            cells.append(colliding_cells.links(i)[0])
    points_on_proc = np.array(points_on_proc, dtype=np.float64)
    ut = fem.Function(V)
    S = sp.csr_matrix((len(points_on_proc), len(ut.x.array)))
    for i, _ in enumerate(ut.x.array):
        ut.x.array[:] = 0.0
        ut.x.array[i] = 1.0
        S[:, i] = ut.eval(points_on_proc, cells).squeeze()
    return S
4

What you are searching for is how to create a point source.
Following is code for getting the basis functions evaluated at a point in the cell:

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

This gives the list of cells on your process, and the corresponding basis values (per row), that you in turn can insert into your matrix/vector.

3 Likes
5

Thank you very much for your code, I will study your code carefully.

6

Thanks again for your codes. Based on your code, I rewrite the codes as follows, which seems work well.

import dolfinx
from dolfinx import fem, la, geometry
from dolfinx.fem import assemble_matrix, assemble_vector
import ufl
from mpi4py import MPI


def construct_measure_matrix(V, points):
    '''
    Input:
        V: function space produced by dolfinx.fem.FunctionSpace
        points: input points, e.g.,
                points = [[x1_1, x1_2, x1_3], [x2_1, x2_2, x2_3]] stands for
                there are two measured points with coordinates:
                x1=(x1_1,x1_2,_x1_3), x2=(x2_1,x2_2,x2_3)
    Output:
        S: a scipy.sparse.csr_matrix

    Remark:
    The following code is based on the codes provided by Jorgen Dokken at the website:
    https://fenicsproject.discourse.group/t/how-to-construct-a-sparse-matrix-used-for-calculating-function-values-at-multiple-points/13081/4
    '''
    nx, dim = points.shape
    assert dim == 3, "The points should be shape Nx3 (N is the number of points)"
    mesh = V.mesh
    sdim = V.element.space_dimension
    # 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)
    ## The extral parameter 1e-6 seems not allowed
    _, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points)
    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))
        ## the command np.asarry(cells) seems leading to some errors
        values = expr.eval(mesh, 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)

    rows = np.zeros(nx*sdim, dtype='int')
    cols = np.zeros(nx*sdim, dtype='int')
    vals = np.zeros(nx*sdim)
    for k in range(nx):
        jj = np.arange(sdim*k, sdim*(k+1))
        rows[jj] = k
        # Find the dofs for the cell
        cols[jj] = V.dofmap.cell_dofs(cells[k])
        vals[jj] = basis_values[0][sdim*k:sdim*(k+1)]

    ij = np.concatenate((np.array([rows]), np.array([cols])), axis=0)
    S = sp.csr_matrix((vals, ij), shape=(nx, V.tabulate_dof_coordinates().shape[0]))
    return S
1 Like
7

Hello,

How can you get the basis derivatives of the specific cell based on your code?

8

Use u.dx(i) in the Expression for the i-th directional derivative.

9

What’s the most straight-forward way to insert values into a vector (say the rhs vector of a LinearProblem object)? I’ve looked at the documentation for the LinearProblem class and I didn’t see any method allowing the insertion of values corresponding to specific cells.

10

You can access the underlying vector in LinearProblem by:
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/generated/dolfinx.fem.petsc.html#dolfinx.fem.petsc.LinearProblem.b

To insert for specific values in specific cells, you would have to be way more specific, as to how the values relate to the degrees of freedom etc, which, in my personal opinion, deserves a new post.

1 Like