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

176

views

1

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
u.vector().set_local(numpy.random.rand(121))
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

6

```
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)
A.get_diagonal(diag_vector.vector())
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()
self.new_grad_func.vector().set_local(new_gradient_array)
self.new_grad_func.vector().apply("insert")
return self.new_grad_func
grad_operator = InvMassGradient(Vgrad, V)
u = Function(V)
#u.vector().set_local(numpy.random.rand(121))
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.