dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

Here's a list of various ways that Function::interpolate can fail.
There are many variants of failures:
- too strict assertion
- missing error checks leading to crash
- missing error checks leading to wrong results
- unexpected behaviour where I expected correct results
Seems like this is not just a quick fix, and I have other things to
do, so I'm bailing out here...

(The code below uses SFC-generated function spaces, but those lines
should be easy to switch with FFC-generated code.)



#include <dolfin.h>
#include "generated_code/Elements.h"

using namespace dolfin;

class Vector2D: public Function
{
public:
  Vector2D(const FunctionSpace & V): Function(V) {}
  Vector2D(): Function() {}

  void eval(double* values, const double* x) const
  {
    values[0] = 3.0 + 0.1*x[0] + 0.01*x[1];
    values[1] = 5.0 + 0.1*x[0] + 0.01*x[1];
    values[0] = x[0] + x[1];
    values[1] = x[0] + x[1];
  }

  unsigned rank() const
  { return 1; }

  unsigned dim(unsigned i) const
  { return 2; }
};

int main()
{
  UnitSquare mesh(1, 1);

  // from a VectorElement("CG", "triangle", 1)
  ElementsForm_av2D::CoefficientSpace_v1_2D V1(mesh);
  // from a VectorElement("CG", "triangle", 2)
  ElementsForm_av2D::CoefficientSpace_v2_2D V2(mesh);


  // Ideally, I'd like to say "make f the interpolant of s" where f
already has a function space attached:
  /*
  Vector2D s;
  Function f(V2);
  f.interpolate(s);
  */


  // But with the currently available signatures, I wanted to do this:
  /*
  Vector2D s; // No function space necessary
  Function f(V2);
  s.interpolate(f.vector(), V2);
  */
  // But this gives:
  /*
terminate called after throwing an instance of 'std::runtime_error'
  what():  *** Assertion (_function_space) [at
dolfin/function/Function.cpp:134 in function_space()]
  */
  // For some reason there's an assertion for a function space in s


  // This works: (sqrt(2^2+2^2) = 2.83 in the corner)
  /*
  Vector2D s(V2);
  Function f(V2);
  s.interpolate(f.vector(), V2);
  */


  // This error should be detected, but it only produces wrong results:
  /*
  Vector2D s(V1);
  Function f(V1);
  s.interpolate(f.vector(), V2);
  */


  // This might be ok (not sure how to interpret s(V1)), but produces
wrong results:
  /*
  Vector2D s(V1);
  Function f(V2);
  s.interpolate(f.vector(), V2);
  */


  // This should fail with an error message, but crashes:
  /*
  Vector2D s(V2);
  Function f(V2);
  s.interpolate(f.vector(), V1);
  */
  /*
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
probably memory access out of range
  */


  File file("f.pvd");
  file << f;

  return 0;
}



Martin