Interpolated function is horizontally flipped

Question

Hey,

I try to interpolate a given function (as expression and function)

def exact_solution(x,t):
   return exp(t-x)*sin(t+x)
u_expression = Expression("exp(t-x[0])*sin(t+x[0])",t=0)

from $t \in (0,1)$ and $\Omega = [0,1]$. In each time step, I use

u_int = interpolate(u_expression,V)
u_int_array = u_int.vector().array()

and save afterwards the nodal values of the interpolated and exact function in arrays

  for i in range(0,nx):
      u_int_nod_vals.append(u_int_array[i])
      u_ex_nod_vals.append(exact_solution(x_array[i],t))
  sol_exact.append(u_ex_nod_vals)
  sol_int.append(u_int_nod_vals)

However, all the values in "sol_exact" are horizontally flipped. At first I thought my output function is the problem, but then the exact function should be flipped as well. Did anybody else encounter this problem?

Here is my entire source code:

from __future__ import division

from dolfin import *
import numpy as np

import matplotlib as mpl
mpl.use( "agg" )
import matplotlib.pyplot as plt
from matplotlib import animation
from math import exp
from math import sin

# ------------------- functions
def save_solution(sol, xmax, title, filename):
  ymax = max(max(l) for l in sol)
  fig = plt.figure()
  axe = plt.axes(xlim=(0, xmax*1.1), ylim=(0, ymax*1.1))
  line, = axe.plot([], [], lw=2, color = 'g')
  plt.xlabel('Position')
  plt.ylabel('Value')
  plt.title(title)

  def init():
    line.set_data([], [])
    return line,

  def animate(i):
    y = np.array(sol[i])
    line.set_data(x_array, y)
    return line,

  anim = animation.FuncAnimation(fig, animate, init_func=init, frames=len(sol), interval=20, blit=True)
  anim.save(filename+'.mp4', fps=15)
# ------------------- functions

nx = 80
xmax = 1
mesh1D = IntervalMesh(nx, 0, xmax)

# definition of solving method
V = FunctionSpace(mesh1D, 'Lagrange', 1)

dt = 1/nx
T=1
t = dt

def exact_solution(x,t):
    return exp(t-x)*sin(t+x)

u_expression = Expression("exp(t-x[0])*sin(t+x[0])",t=0)

#constants for ploting
x_points = []
for i in range(0,nx):
  x_points.append((xmax*i)/(1.0*nx)) 

x_array = np.array(x_points)

sol_exact = []
sol_int = []

while t <= T:
    u_ex_nod_vals = []

    u_int = []
    u_int_nod_vals = []

    u_int = interpolate(u_expression,V)
    u_int_array = u_int.vector().array()

    for i in range(0,nx):
      u_int_nod_vals.append(u_int_array[i])
      u_ex_nod_vals.append(exact_solution(x_array[i],t))

    sol_exact.append(u_ex_nod_vals)
    sol_int.append(u_int_nod_vals)

    t += dt
    u_expression.t=t


save_solution(sol_int,xmax, 'interpolated solution', 'sol_interpolate')
save_solution(sol_exact,xmax, 'Exact solution', 'sol_exact')

Comments

I think

 u_int_array = u_int.vector().array()

is causing the problem by "rearranging" the array. A simple workaround would be sort/invert it, but I would like to know what's exactly happening, i.e. is this the expected behaviour?

Answer 1

In this special case, you can set

parameter["reorder_dofs_serial"] = False

at the beginning of the code.

In general however, you can not assume anything on the ordering of the dofs (your u_int_array). You can however access the coordinates of the dofs through

V.dofmap().tabulate_all_coordinates(mesh1D)

and work with that for more general cases.

Comments

Ok, I get it, so I would need to implement a function which maps my coordinates to the corresponding index of the dof, and then I could call u_int_array[index] and would get the correct value, am I right?

---

Well, as Maximillian points out below, there is this function vertex_to_dof_map to find this mapping. However, this only works for a FunctionSpace with dofs only on vertices.

---

Thanks a lot, I totally forgot about the data structures and dofs etc.
As I'm going to use my code for an example in 2D and 3D as well I'm going to use

 V.dofmap().tabulate_all_coordinates(mesh1D)

which, if I understand correctly, has, in 2D, the two coordinates of the first dof, then the two coordinates of the next dof and so on, in 3D I would have to extract 3 values, corresponding to the three coordinates. So, after having extracted the coordinates in order of the dofs, I just have to iterate over u_int.vector() and write the value in my array sol_int.
Thank you two a lot, it helped me to understand Fenics a little more :-)

Answer 2

In general, dolfin will reorder the degrees of freedom of your function for efficiency to speed up calculations (especially in parallel). You can avoid this (in serial calculations) by setting

parameters['reorder_dofs_serial'] = False

Alternatively, there is a convenience function which gives you a mapping between the vertex indices and the indices of the degrees of freedom.

Try to replace your original line:

u_int_array = u_int.vector().array()

with the following:

u_int_array = u_int.vector()[vertex_to_dof_map(V)].array()

This should reorder the vector entries into the "expected" order and produce the correct result.

Comments

Thanks a lot Max, but I'm trying to use V.dofmap().tabulate_all_coordinates(mesh1D) as I want to apply my program to more general cases as well.