Export basis functions in vectorized form

1

I want to export all basis functions evaluated at a grid (different than the FE-mesh).
I came up with two ways to export the basis functions.

Algorithm 1

import fenics as fe
import numpy as np
mesh = fe.UnitSquareMesh(10, 10, 'crossed')
x = np.linspace(0, 1, 60)
y = np.linspace(0, 1, 60)
V = fe.FunctionSpace(mesh, 'CG', 2)
F = fe.Function(V)

t = time.time()
for ib in range(len(F.vector())):
    F.vector()[:] = 0
    F.vector()[ib] = 1
    for i in range(len(x)):
        for j in range(len(y)):
            val = F(x[i], y[j])

Algorithm 2

import fenics as fe
import numpy as np
mesh = fe.UnitSquareMesh(10, 10, 'crossed')
x = np.linspace(0, 1, 60)
y = np.linspace(0, 1, 60)
V = fe.FunctionSpace(mesh, 'CG', 2)
el = V.element()
bbt = fe.BoundingBoxTree()
bbt.build(mesh)
for i in range(len(x)):
    for j in range(len(y)):
        cell_id = bbt.compute_first_entity_collision(fe.Point(x[i], y[j]))
        cell = fe.Cell(mesh, cell_id)
        coordinate_dofs = cell.get_vertex_coordinates()
        val = el.evaluate_basis_all(np.array([x[i], y[j]]), coordinate_dofs, cell_id)

Algorithm 2 is much more efficient than Alg. 1 (800 times faster). But it’s somehow inefficient because of the nested for-loops. Is there a way to evaluate the basis functions/function in vectorized form? In essence some built-in function that takes an array or list of points as input, instead of a single point.

2

If the grid that you’re evaluating basis functions on is (or can be) represented as a DOLFIN Mesh, then you can use the PETScDMCollection functionality to do this. I wrote up a short example of this in the following thread:

Each row of the matrix A gives the evaluations of all basis functions from the space V1 at a node from the space V2.

2 Likes
3

Thank you, Kamensky,
with PETScDMCollection the code is faster and looks a bit more elegant.

The final code looks like this.

import fenics as fe
import numpy as np
import scipy.sparse as sp

mesh1 = fe.UnitSquareMesh(10, 10, 'left/right')
V1 = fe.FunctionSpace(mesh1, 'CG', 2)

# Function required in the minimal example to prevent a PETSc initialize error
F = fe.Function(V1)
# Evaluation grid
y = np.linspace(0, 1, 40)
x = np.linspace(0, 1, 20)
# New mesh where basis functions from V1 are evaluated
mesh2 =  fe.UnitSquareMesh(19, 39, 'left/right')
# I sorted the mesh2 with in y
(mesh2.coordinates()[:, 1]).sort(axis=0)
xf, yf = mesh2.coordinates()[:,  0], mesh2.coordinates()[:, 1]
# Changing the mesh 2 with any grid I want
yf[:] = np.repeat(y, len(x))
# Create second function space
V2 = fe.FunctionSpace(mesh2, 'CG', 1)
# Transfer matrix
A = fe.PETScDMCollection.create_transfer_matrix(V1,V2)

# Hack solution to rearrange the A matrix in y and x
arg_x = np.argsort(V2.tabulate_dof_coordinates()[:, 0])
coord_x = V2.tabulate_dof_coordinates()[arg_x, :]
arg_y = np.argsort(np.reshape(coord_x[:, 1], (len(x), len(y))))
arg_y = (arg_y.T + np.arange(len(arg_y))*(len(y))).T.flatten()
row, col, val = fe.as_backend_type(A).mat().getValuesCSR()
A = sp.csr_matrix((val, col, row))

# Final basis functions from V1 evaluated at y, x
w = (A[arg_x])[arg_y].T

The only issue seems to be when the mesh1 is ‘crossed’, where this solution differs from my previous implementation. With ‘left’, ‘right’, and ‘left/right’ meshes the different basis functions export routine produce the same results.

4

The reason that "crossed" gives a different result is that the set of mesh vertices (and thus the nodal points of the finite element basis for the \text{CG}_1 space V2) is different. With "crossed", there is an extra vertex in the center of each square “cell” of the structured mesh.

EDIT: Okay, I see now you said the discrepancy occurs when mesh1 is "crossed", not mesh2. Is that a typo? If I understand what you’re doing correctly, then I would not expect a difference in behavior, unless there is a bug in at least one of the methods.

5

Thanks Kamensky.

Yes, I meant mesh1, but I actually had a small bug in my other implementation.
Now I double checked it and the PETScDMCollection functionality works for the crossed mesh as well.