Interpolation of Indexed Functions of Mixed Elements

1

Hello,

I again have problems with Indexed functions.
I want to do the quicker interpolate instead of the slower project on Indexed functions.
However I get an error that Indexed does not have an _cpp_object.
The only way that I know to circumvent this is by projecting, but that is what I want to avoid.

Consider the minimal example (of course the use here does not make much sense):

from fenics import *

mesh = UnitSquareMesh(3, 3)

S1 = FiniteElement('CG', mesh.ufl_cell(), 1)
S1S1 = MixedElement([S1, S1])
Y = FunctionSpace(mesh, S1S1)
Y1 = Y.sub(0).collapse()
Y2 = Y.sub(0).collapse()

y = project(Expression(('x[0]*x[1]', 'x[0] + x[1]'), degree=2), Y)
# this works
y1 = project(y[0], Y1)
# this gives:
#    AttributeError: `Indexed` object has no attribute `_cpp_object`
y2 = interpolate(y[0], Y1)

Any help appreciated.

2

Hi, using a FunctionAssigner should work. By first interpolating the expressions to the non-Mixed spaces, then using a FunctionAssigner to assign each of them to the mixed space. Demo at: https://bitbucket.org/fenics-project/dolfin/raw/ec57db53f2b13f1768214829e8b4f80fc32205cf/python/demo/undocumented/sub-function-assignment/demo_sub-function-assignment.py

3

Hi,
in my use case I do not actually have the expression available anymore. The code should just demonstrate the error that I get in my case. I also do not want to interpolate on the same function space that the expression was originally projected on. This was also just done here for simplicity. So as far as I understood your solution it does not really help.

Generally I found it wired that I can call project on an object that I cannot call interpolate on.

Ty, anyway for your time.

4

Hi,
How about:

u,v = y.split()
y2 = interpolate(u, Y1)
5

Ok, that sort of worked.

In fact my situation was even a little more involved, since I had to deal with a whole array of such functions.
See https://fenicsproject.discourse.group/t/problems-with-numpy-arrays-of-functions-over-mixedelement/1171 for a description of the setup (y below refers to arr_Y of the there).

With your help it could now change my from

y1_cmp = numpy.array([project(yt[0], Y1_cmp) for yt in y], dtype=object)
y2_cmp = numpy.array([project(yt[1], Y2_cmp) for yt in y], dtype=object)

to

y_tmp = [yt[0].ufl_operands[0].split() for yt in y]
y1_cmp = numpy.array([interpolate(yt[0], Y1_cmp) for yt in y_tmp])
y2_cmp = numpy.array([interpolate(yt[1], Y2_cmp) for yt in y_tmp])

I must admit though that I am still confused why this is so “complicated”. It uses the first entry in ufl_operands of yt[0] which is relatively confusing to parse in my opinion.
Why is it necessary to use split here instead of just indexing?
Is my expectation wrong that I could use interpolate basically everywhere where I could use project?

6

They key here is that .split() return a dolfin.Function, which then in turn can be either interpolated or projected.
The indexed version, equivalent of calling u,v = split(y) return an ufl.index.Indexed which is not a dolfin object that can only be handled in assembly routines (when its part of a ufl-form, as in projection).
There is a lengthy discussion on the subject here.
So your assumption isn’t wrong, you just have to be careful on what objects you send into interpolate.

7

I do not want to rob any more of your time so do not feel obliged to reply.
I mostly write this down to strengthen my own understanding.

So here is what I understood.

In my example y[0][0] (being the result of __getitem__) is of type ufl.indexed.Indexed, while
y_tmp[0][0] (being the result of .split()) is of type dolfin.function.function.Function.

My experience now is reflected by this table

Indexed Function
project OK OK
interpolate NOT OK OK

That is my assumption is indeed wrong.

Indeed the docstrings of interpolate and project only mention the second column anyway.
According to them both functions should work with either dolfin.function.function.Function or dolfin.functions.expression.Expression (though the docstrings do not really argree on the name of the submodule of dolfin [“functions” vs. “function”].)

So project(Indexed) does work but somehow it “should” not. It breaks the consistency and is not covert by the docstrings.

Anyway ty a lot for your help.

8

The reason for yy=project(y[0],V) working is that it Is equivalent to

u,v =TrialFunction(V), TestFunction(V)
a=inner(u,v)*dx
l=inner(y[0],v)*dx
yy=Function(V)
solve(a==l, yy)
1 Like