matrix from `project`, evaluation at many points at once

4 months ago by
I have an FEM function `u` of which I change the coefficients a large number of times, and for each of the instances I'll have to evaluate the gradient and a fixed set of points:
from dolfin import *
import numpy

XY = numpy.random.rand(77, 2)

mesh = UnitSquareMesh(10, 10)

V = FunctionSpace(mesh, 'CG', 1)
Vgrad = VectorFunctionSpace(mesh, 'DG', 0)

u = Function(V)
for _ in range(1000):
    # set u to something

    g = project(grad(u), Vgrad)
    jac = numpy.array([g(x, y) for x, y in XY])
Now, I would like to speed this up. Two questions come to mind:

  • Can the `project` operator be assembled into a matrix?
  • Can `g` be evaluated at all points XY at once?
Community: FEniCS Project
Do you want a matrix that projects onto your DG 0 space or onto any space?
written 4 months ago by Miguel  
DG0 is fine for now, but the operator really only needs the gradient only at a handful of points (namely `XY`).
written 4 months ago by Nico Schlömer  

1 Answer

4 months ago by
from dolfin import *
import numpy

mesh = UnitSquareMesh(10,10)
Vgrad = VectorFunctionSpace(mesh, 'DG', 0)
V = FunctionSpace(mesh, 'CG', 1)

class InvMassGradient(object):
    def __init__(self, Vgrad, V):
        u_elect = TrialFunction(V)
        v_vec = TestFunction(Vgrad)
        u_vec = TrialFunction(Vgrad)

        # Define the operator to get the gradient
        L_operator = inner(v_vec, grad(u_elect))*dx
        # Mass matrix
        a_vec = inner(v_vec, u_vec)*dx

        A = assemble(a_vec)
        L_vec = assemble(L_operator)
        diag_vector = Function(Vgrad)
        self.L_vec = L_vec
        self.diag_vector = diag_vector
        self.new_grad_func = Function(Vgrad)
    def __rmul__(self,x):
        return NotImplemented
    def __mul__(self,x):
        new_gradient = self.L_vec * x.vector()
        new_gradient_array = new_gradient.array()
        new_gradient_array = new_gradient_array / self.diag_vector.vector().array()
        return self.new_grad_func

grad_operator = InvMassGradient(Vgrad, V)
u = Function(V)
u.interpolate(Expression("x[0]", degree=1))
g = grad_operator * u

g is the value of the gradient at the centroid of the elements.

Please login to add an answer/comment or follow this question.

Similar posts:
Search »