#include "dolfin.h"

using namespace dolfin;

// PPDE DOLFIN example #3: mesh representation, input/output

real f(real x, real y)
{
  return x * x + y * y;
}

int main()
{
  dolfin_set("output", "plain text");

  Mesh mesh("mymesh.xml.gz");

  int cellcounter = 0;

  for(CellIterator cell(mesh); !cell.end(); ++cell)
  {
    cellcounter++;
  }

  cout << "Cells: " << cellcounter << endl;

  Vector xi(mesh.noNodes());

  for(NodeIterator n(&mesh); !n.end(); ++n)
  {
    Node& node = *n;

    Point& p = node.coord();
    int id = node.id();

    xi(id) = f(p.x, p.y);
  }

  cout << "xi: " << endl;
  xi.show();

  Function U(mesh, xi);

  File solutionfile("mysolution.m");

  solutionfile << U;
}