// Oldham, Jeffrey D.
// 2001Apr20
// Pooma

#include <iostream>
#include <stdlib.h>
#include <cmath>
#include "Pooma/NewFields.h"

int main(int argc, char *argv[])
{
  // Set up the Pooma library.
  Pooma::initialize(argc,argv);

  // Create a simple layout.
  const unsigned Dim = 2;		// Work in a 2D world.
  const unsigned nHorizVerts = 4;	// number of horizontal vertices
  const unsigned nAngles = 3;		// number of angles
  Interval<Dim> vertexDomain;
  vertexDomain[0] = Interval<1>(nHorizVerts);
  vertexDomain[1] = Interval<1>(nAngles);
  DomainLayout<Dim> layout(vertexDomain, GuardLayers<2>(1));

  // Preparation for Field creation.
  Vector<Dim> origin(0.0);
  Vector<Dim> spacings(1.0,1.0);
  typedef UniformRectilinear<Dim, double, Cartesian<Dim> > Geometry_t;
  typedef Field<Geometry_t, double,      Brick> Fields_t;

  Cell cell;
  Vert vert;
  AllFace allFace;
  Fields_t internalEnergy  (cell, layout, origin, spacings);
  Fields_t coordinates	   (vert, layout, origin, spacings);

  // Initialize the Fields.
  for (unsigned x = 0; x < nHorizVerts; ++x)
    for (unsigned y = 0; y < nAngles; ++y)
      coordinates (x,y) = x+y;
  for (int x = 0; x < nHorizVerts-1; ++x) // Removing "-1" yields different results.
    for (int y = 0; y < nAngles-1; ++y) // Removing "-1" yields different results.
      internalEnergy (x,y) = -(x+y);

  std::cout << "initial internal energy:\n" << internalEnergy << std::endl;
  std::cout << "initial coordinates:\n" << coordinates << std::endl;

  coordinates += internalEnergy; // THE ONLY INTERESTING LINE

  std::cout << "final coordinates:\n" << coordinates << std::endl;
  Pooma::finalize();
  return EXIT_SUCCESS;
}