Evaluate an expression from the solution

Question

I calculate a solution to two different (sets of coupled) equations. I know that the solution to the second equation is close to the solution of the first equation plus some given, analytical function. What I do is I first compute the solution to the first equation, i.e.:

solution = Function(solutionspace)
solve(F0 == 0, solution, [bc1,bc2,...])

After this, I solve the second equation. Since the solution to the first equation is now stored in the solution variable, using it as an initial guess can be done by just solving the second equation to the same solution variable, i.e.:

solve(F == 0, solution, [bc1,bc2,...])

What I would like to do is to add a constant function to the solution between these two steps. I tried:

solution.sub(0).assign(solution.sub(0) + f)

where f is a given function defined with Expression(...). Similarly for other three subcomponents of the solution. This did not work. What is the correct way to do this?

Answer 1

I think you need to interpolate f before adding. I am not familiar with .sub command but what you intend to do can be done with these two lines.

ff = interpolate(f , solutionspace)
solution.vector()[:] = solution.vector() + ff.vector()

Comments

In my problem, solution is a four component solution (i.e. includes four solutions to four four coupled equations), but I want to add f to only one of them. Your solution looks like I add f to each one of them.

---

You might already have tried this. Does it work?

ff = interpolate(f , solutionspace)
solution.sub(0).vector()[:] = solution.sub(0).vector() + ff.vector()

---

Doesn't work. I get:
Error: Unable to interpolate function into function space.
Reason: Rank of function (0) does not match rank of function space (1).
Where: This error was encountered inside FunctionSpace.cpp.

I also tried interpolating to solutionspace.sub(0) instead of solutionspace. That gives:
Error: Unable to create function.
Reason: Cannot be created from subspace. Consider collapsing the function space.
Where: This error was encountered inside Function.cpp.

---

OK. How about creating another scalar valued function space, then interpolating f on that, and then adding.

V1 = FunctionSpace(mesh, "CG", 1)
ff = interpolate(f , V1)
solution.sub(0).vector()[:] = solution.sub(0).vector() + ff.vector()

---

I get the error:
Error: Unable to access vector of degrees of freedom.
Reason: Cannot access a non-const vector from a subfunction.

I also tried using split() to access one of the subfunctions, but I get the same error. I have to say that the way Fenics represents multi-component solutions to coupled equations feels very counterintuitive. It seems there is no easy way to access or modify just one of the solution components.

---

Did you try vector valued expression? Something like

f = Expression(("1.0", "0.0", "0.0", "0.0"))
ff = interpolate(f , solutionspace)
solution.vector()[:] = solution.vector() + ff.vector()

Example in this link might be more helpful.

http://fenicsproject.org/documentation/dolfin/1.4.0/python/demo/documented/hyperelasticity/python/documentation.html#index-0

---

This seems to actually work. Thanks a lot!

---

Wonderful! Thanks for your vote. Feels good to get points. It does not matter what they mean. :)