Can't interpolate data from function to a slice plane. Error: Unable to create mapping of degrees of freedom

Question

Hi,

I'm trying to compute average function values on a slice plane.

A python example was described here earlier by mikael-mortensen.
Python example works like a charm.

Link to the example:
http://fenicsproject.org/qa/4517/how-can-i-obtain-plane-averages

I created a cpp version of "compute average function values on a slice plane", but I caught an error,

The error raises (with dolfin::LagrangeInterpolator) when I try to interpolate data from one function (3d mesh) to another function (2d sub mesh).

I use ubuntu 12.04 64-bit, dolfin 1.4 (in the dolfin 1.5 source code looks the same for the dolfin::LagrangeInterpolator class).

The error code:

*** -------------------------------------------------------------------------
*** Error:   Unable to create mapping of degrees of freedom.
*** Reason:  Geometric dimension of the UFC dofmap (dim = 2) and the mesh (dim = 3) do not match.
*** Where:   This error was encountered inside DofMap.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Place where the error happens:

lagrange_interpolator.interpolate(u, u_0);

My cpp code below: 1 cpp file and 3 ufl files. You can convert ufl files to a header files by the following command:

ffc -l dolfin file_name.ufl

fenics_average_plane.cpp.

#include <iostream>

#include <dolfin.h>

#include "average_scalar_2d.h"
#include "scalar_2d.h"
#include "scalar_3d.h"


// Define test Expression.
class MyExpression : public dolfin::Expression {
public:

  void eval(dolfin::Array<double>& values, const dolfin::Array<double>& x) const {
    //values[0] = 10;
    values[0] = x[2];
  }

};

// Define Inlet SubDomain.
class Inlet : public dolfin::SubDomain {
  bool inside(const dolfin::Array<double>& x, bool on_boundary) const {
    return x[2] < DOLFIN_EPS; 
  }
};

int main() {

  dolfin::parameters["allow_extrapolation"] = true;

  // Create a mesh.
  dolfin::UnitCubeMesh mesh(10, 10, 10);

  // Create a function space based on the mesh (3d).
  scalar_3d::FunctionSpace Q(mesh);

  // Create an Expression.
  MyExpression u_0_expr;

  // Create a function which will be interpolated to a slice plane.
  dolfin::Function u_0(Q);

  u_0.interpolate(u_0_expr);

  // Create a BoundaryMesh.
  dolfin::BoundaryMesh boundary_mesh(mesh, "exterior");

  // Create a MeshFunction.
  dolfin::CellFunction<std::size_t> subdomains(boundary_mesh, 0);

  // Create an Inlet subdomain.
  Inlet inlet;

  // Mark the Inlet subdomain.
  inlet.mark(subdomains, 1);

  // Create a sub mesh (a slice plane).
  dolfin::SubMesh sub_mesh(boundary_mesh, subdomains, 1);

  // Get the slice plane coordinates.
  std::vector<double> sub_mesh_coord = sub_mesh.coordinates();

  // Get a sub mesh coordinates size.
  std::size_t sub_mesh_coord_size = sub_mesh_coord.size();

  // Get a geometry dimension.
  std::size_t geom_dim = mesh.geometry().dim();  
  std::cerr << "geom_dim = " << geom_dim << std::endl;

  // Move the slice in the z-direction by changing (redefine) z coordinate in the slice plane points.
  for (std::size_t i = geom_dim - 1; i < sub_mesh_coord_size; i += geom_dim) {
    sub_mesh_coord[i] = 0.5;
  }

  // Create a function space based on the sub mesh.
  scalar_2d::FunctionSpace Q_sub(sub_mesh);
  dolfin::Function u(Q_sub);

  // Create a LagrangeInterpolator to interpolate data from u_0 to u (defined on the sub mesh).
  dolfin::LagrangeInterpolator lagrange_interpolator;
  lagrange_interpolator.interpolate(u, u_0);

  // Other interpolation option (alternative).
  //u.interpolate(u_0);

  // Plot functions to check the results.
  dolfin::plot(u_0, "u_0");
  dolfin::plot(u, "u");  
  dolfin::interactive();

  // Evaluate functional on the sub mesh. 
  // Compute average function value on the slice plane.
  average_scalar_2d::Functional M(sub_mesh, u);
  double u_average = assemble(M);  
  std::cerr << "u_average = " << u_average << std::endl;

  // Output data to a file.
  dolfin::File file_u("output_data/u.pvd");
  file_u << u;  

  std::cerr << "The end" << std::endl;

  // Stop console in the Debug mode.
  std::cin.get();
}

average_scalar_2d.ufl

cell = triangle

element = FiniteElement("Lagrange", cell, 1)

u = Coefficient(element)

M = u*dx

scalar_2d.ufl

cell = triangle

element = FiniteElement("Lagrange", cell, 1)

scalar_3d.ufl

cell = tetrahedron

element = FiniteElement("Lagrange", cell, 1)

Thanks in advance!

Best regards,
Maksim