Modifying and assigning solution

I/O
1

Is there a way to somehow modify values of computed solution?
Normally we’d do the following:

u_n.assign(u)

But what if I’d like to modify the solution u before assigning it?
Example - multiply all values by 2, which obviously doesn’t work:

u_n.assign(2*u)

More advanced version would be how to perform such modification according to defined MeshFunction?

2

Hey, you can do such things by directly working on the dof values of the Function. For your example, this would look as follows

u_n.vector()[:] = 2*u.vector()[:]

Note that this works dof wise, and gives the correct result for Lagrangian elements. For other element types (or more complex modifications that, e.g., depend on spatial coordinates) you can use Expressions, which would look as follows

u_expr = Expression('2*u', degree=n, u=u)
u_n = interpolate(u_expr, V)

where the degree has to be chosen appropriately. Instead of interpolate you could also chose to project the expression, but I prefer interpolating it, as it gives the nodally exact values for Lagrangian elements (and if the right degree was chosen)

3

I’m doing something similar using version 2019.1.0. I’m replying to this post because I believe there may be a bug in assigning a subspace of a vector function space. I find that assigning a subspace will not work; you need to assign the entire vector function space simultaneously.

I’m finding that using assign to multiply all values by 2 works:

import fenics as fe
import numpy as np
mesh = fe.UnitIntervalMesh(1)

V = fe.FunctionSpace(mesh, 'Lagrange', 1)
u = fe.Function(V)

u = fe.interpolate(fe.Constant(1.), V)
u.assign(2.*u)
assert np.isclose(u.compute_vertex_values(mesh)[0], 2.), 'Incorrect copying' # passes

However, when I do this for a vector function space, and just try to modify one of the subspaces

import fenics as fe
import numpy as np
mesh = fe.UnitIntervalMesh(1)
V = fe.VectorFunctionSpace(mesh, 'CG', degree=1, dim=2)
u = fe.Function(V)
u = fe.interpolate(fe.Constant((1.,1.)), V)
u_copy = u.copy(deepcopy=True)
u.sub(0).assign(2.*u_copy.sub(0)) # also fails if use u.sub(0).vector()[:] = 2.*u_copy.sub(0).vector()
assert np.isclose(u.sub(0).compute_vertex_values(mesh)[0], 2.), 'Incorrect copying for index 0' # passes
assert np.isclose(u.sub(1).compute_vertex_values(mesh)[0], 1.), 'Incorrect copying for index 1' # fails

The second test fails, the whole vector is = 2.
This tells me that there might be a bug in the assign function; perhaps subspace assigns should not be allowed

If I instead assign the whole vector function space at the same time, as below, the test works

import fenics as fe
import numpy as np
mesh = fe.UnitIntervalMesh(1)
V = fe.VectorFunctionSpace(mesh, 'CG', degree=1, dim=2)
u = fe.interpolate(fe.Constant((1.,1.)), V)

u_expr = fe.Expression(('2*u','1'), degree=1, u=u.sub(0))
u = fe.interpolate(u_expr, V)
assert np.isclose(u.sub(0).compute_vertex_values(mesh)[0], 2.), 'Incorrect copying for index 0'
assert np.isclose(u.sub(1).compute_vertex_values(mesh)[0], 1.), 'Incorrect copying for index 1'
4

For vector-spaces, and mixed spaces, you should use FunctionAssigner, see for instance this test on how to assign to a vector space

1 Like