This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Translating (and stitching) functions

0 votes

I'm currently working on a multi-scale method, which works roughly as follows. The initial problem is set over say, $\Omega = [0,1]\times[0,1]$, which is meshed with a coarse structured mesh (built with UnitSquareMesh). In each rectangular cell $C$ of this mesh, one has to solve a sub-problem $P(C)$, which gives a solution $u_C$. This is done on a fine mesh (which meshes $C$). As several quantities/functions are cell independent, I solve all sub-problems on the same mesh, over say $[0, \epsilon] \times [0, \epsilon]$. Then, the $u_C$ s have to be stitched together to recover a function over $\Omega$.

What I'm trying to do at the moment is:

  1. Translate $u_C$ so that they are defined on the correct domain
  2. Extend them by zero to functions of $\Omega$.
  3. (Add them together, disregarding continuity issues)

I don't know how to do 1., more precisely, how can I define v in the code below so that $v(x, y) = u(x-1, y-1) (= x-1)$ ?

mesh = UnitSquareMesh(10, 10)
translated_mesh = UnitSquareMesh(10, 10)
translated_mesh.coordinates()[:] = translated_mesh.coordinates()+1

V = FunctionSpace(mesh, "P", 1)
tV = FunctionSpace(translate_mesh, "P", 1)

u = Function(V)
u.interpolate(Expression('x[0]'))

v = Function(tV)
# v = ?

I have an idea about 2. although I am unable to make @mikael-mortensen's interpolation function work. Have I missed any obvious option/way to extend a function by zero over a non-matching domain, à la parameters['allow_extrapolation'] = True, but without trying to be smart about the unknown values for the new functions ?

Thanks for any advice you might have, otherwise, thanks for reading this far ;)

asked Sep 29, 2014 by gjankowiak FEniCS Novice (880 points)
edited Sep 29, 2014 by gjankowiak

2 Answers

0 votes
 
Best answer

I'm answering as per my comment on MiroK's answer.
For 1.,

v.vector()[:] = u.vector()

is the way to go. This might require some projection if u is not actually a Function, but rather a Sum, Product or similar.

For 2., use Mikael Mortensen's interpolation function.

answered Feb 12, 2015 by gjankowiak FEniCS Novice (880 points)
0 votes

Hi, this should help with (1) on your list

from dolfin import *

code = '''
class MyFunc : public Expression
{
public:
  std::shared_ptr<GenericFunction> u;

  MyFunc() : Expression() { }

  void eval(Array<double>& values, const Array<double>& x) const
  {
    assert(u);

    // Translate
    Array<double> x_translated(2);
    x_translated[0] = x[0] - 1;
    x_translated[1] = x[1] - 1;

    u->eval(values, x_translated);
  }
};
'''

mesh = UnitSquareMesh(10, 10)
u = Expression('x[0]')

translated_mesh = UnitSquareMesh(10, 10)
translated_mesh.coordinates()[:] = translated_mesh.coordinates()+1
v = Expression(code)
v.u = u

plot(u, mesh=mesh, title='u')
plot(v, mesh=translated_mesh, title='translated u')
interactive()
answered Sep 29, 2014 by MiroK FEniCS Expert (80,920 points)

Thanks for the anwser, I was hoping for something simpler though. While keeping searching, I realized I could simply do

v.vector()[:] = u.vector()

or

v.vector()[:] = u.compute_vertex_values()

Of course this needs the same dofmap (?) for the two meshes/function spaces.

...