Derivative of a function wrt expansion coefficients

Question

Hello everybody,

I've started using Fenics for my PhD, and after being able to resolve all issues myself so far, I got stuck on a seemingly simple problem. As a rookie, I would appreciate any help or pointer in the right direction.

Background: I have defined a measure for the deviation of the solution (vector valued function) u from a target at a discrete set of points

$$ p = \sum\limits_{k=1}^K \beta_k \left( \left( u_x \left( x_k \right) - \hat{y}_{x,k} \right)^2 + \left( u_y \left( x_k \right) - \hat{y}_{y,k} \right)^2 \right) $$

I need to evaluate the derivative wrt to the N expansion coefficients of u at K values x_k. I need to do so often, and cheaply. To this end I would like to assemble a matrix containg the derivatives of u wrt expansion coefficients U_j at x_k

$$ \frac{\partial}{\partial U_j} u \left( x_k \right) = \frac{\partial}{\partial U_j} \left( \sum\limits_{i=1}^N \phi_i \left( x_k \right) U_{i} \right) = \phi_j \left( x_k \right) \qquad \forall x_k, k=1,...,K $$

for the x and y "subspace" of a Function, separately (excuse the sloppy notation). Or at least I would need to be able to keep account which dof belongs to x and y. This would allow me to get the gradient by a matrix-vector multiplication.

Approaches considered so far:

• Set all but the j-th vector value of a function to zero and evaluate the function at all K points. Issue : prohibitive, as it takes N*K function evaluations (K approx N, i.e. N^2).
• I'm able to evaluate all basis functions of an element in any given cell, but that doesn't help. Looking at the C++ implementation eval() of a Function, there's probably some way to do this
• I was also not able to piece together a solution ussing differentiate() or diff(), or think of a creative way to assemble this KxN matrix

Maybe a little code example:

mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh,'P',2)
F = Function(V)

x_k = (0.5,0.5)
u_k = F(x_k)

=> How does u_k[0] and u_k[1] depend on the expansion coefficients F.vector(), with correct ordering of the dofs?

In case anybody finds the time to help me out, thanks a lot.

Best Regards

Answer 1

Hi, below is an example of how to build the differentiation matrix over scalar space. To do this for vector/tensor-valued or mixed spaces is not much more involved and is a good learning exercise. Note that the code is serial only.

from dolfin import *
from petsc4py import PETSc
import numpy as np


def diff_coeff_matrix(V, points):
    mesh = V.mesh()

    element = V.element()
    basis_values = np.zeros(element.space_dimension())

    nrows = len(points)
    ncols = V.dim()
    mat = PETSc.Mat().createAIJ(size=(nrows, ncols))
    mat.setUp() 
    mat.assemblyBegin()

    dmap = V.dofmap()
    tree = mesh.bounding_box_tree()
    for row, pt in enumerate(points):
        cell_id = tree.compute_first_entity_collision(Point(*pt))

        # Point missing
        if cell_id >= mesh.num_cells(): continue

        cell = Cell(mesh, cell_id)
        vertex_coordinates = cell.get_vertex_coordinates()
        cell_orientation = cell.orientation()

        element.evaluate_basis_all(basis_values, pt, vertex_coordinates, cell_orientation)

        row_indices = [row]
        col_indices = dmap.cell_dofs(cell_id)

        mat.setValues(row_indices, col_indices, basis_values, PETSc.InsertMode.INSERT_VALUES)
    mat.assemblyEnd()
    return PETScMatrix(mat)

# ----------------------------------------------------------------------------------------

if __name__ == '__main__':
    mesh = UnitSquareMesh(3, 3)
    V = FunctionSpace(mesh, 'CG', 8)

    # The matrix should have an identity property if points are
    # dof coordinates
    points = V.tabulate_dof_coordinates().reshape((V.dim(), -1))
    # Let's test it 
    I = diff_coeff_matrix(V, points)
    assert I.size(0) == I.size(1)

    x = as_backend_type(I).mat().createVecRight() 
    y = x.copy()
    x.setRandom()

    x, y = map(PETScVector, (x, y))

    print x.norm('l2'), y.norm('l2')
    I.mult(x, y)
    print (x-y).norm('l2')

Comments

Thank you!

I discarded this approach for some stupid invalid reason. And you're right, the generalisation to vector spaces is pretty simple, as soon as you have figured out that element.evaluate_basis_all() returns an array of length V.dim()*element.space_dimension() and not just element.space_dimension(), where every basis function is already is evaluated for every vector component. Unfortunately, the lower level functions and APIs of fenics / dolfin are scarcely documented.

---

I recommend reading the source code directly. The entire FEniCS stack is mirrored on github which, unlike bitbucket, let's you search the code base.