What is the best way of performing a "pointwise" operation on (one or multiple) functions defined on a given mesh?
For example, given a vector field 'v(r)' and a scalar field 'a(r)' (where 'r' denotes the position), how can I scale each vector v(r) by the corresponding scalar a(r)?
This example is deliberately trivial to keep things simple, but one could imagine more complex operations like computing the cross product of two vector fields, or even rotating each vector by an angle (specified by another scalar field 'theta') about a certain axis (given by yet another vector field). You get the idea.
Here is a code example for the simple example above:
from dolfin import *
import numpy as np
mesh = BoxMesh(0, 0, 0, 100, 100, 2, 50, 50, 1)
V = VectorFunctionSpace(mesh, 'Lagrange', 1, dim=3)
W = FunctionSpace(mesh, 'Lagrange', 1)
v = interpolate(Expression(['1', '2', '3']), V)
a = Function(W)
# Fill 'a' with some random values
n = W.dim()
a.vector().set_local(np.random.random(n))
### How to scale v by a point-wise?!?
Ideally, I'm imagining something along the following lines:
def f(b, w):
# The arguments 'b' and 'w' should be a single scalar
# and a numpy array of length 3, respectively.
return b*w
# Apply the function f point-wise to the dolfin Functions
# a and v defined above
apply_pointwise(f, a, v)
I can certainly resort to doing low-level operations directly on the function vectors, similar to the following:
# The code below assumes a certain ordering
parameters.reorder_dofs_serial = False
# Extract numpy arrays
a_arr = a.vector().array()
v_arr = v.vector().array()
v_arr_3xN = v_arr.reshape(3, -1)
# Compute the vertex-wise result
a_times_v = a_arr * v_arr_3xN
# Convert the result back to a dolfin.Function
result = Function(V)
result.vector().set_local(a_times_v.ravel())
Note that the above code will not work with dof reordering enabled, since the ordering of the entries of v_arr will be different. There are certainly workarounds for that, but I always feel uneasy when dealing with array entries directly because I don't fully understand the implications of dof reordering (yet). For example, one question I've had for a while is: will degrees of freedom belonging to the same "site" (e.g. to a vertex for CG1 elements) always be "kept together" and appear in the expected order in v.vector().array()? For example, with v defined as above, will the entries ['1', '2', '3'] always appear in this order?
In any case, I'm sure there must be a more high-level approach to point-wise operations which avoids dealing with vector elements directly. What is the recommended way of doing this? Any hints would be highly appreciated. Many thanks in advance!
