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

// Chevron Kernel Written Using POOMA's Proposed Field Abstraction

#include <iostream>
#include <cstdlib>
#include "Pooma/NewFields.h"

// This program implements "Implementation of a Flux-Continuous Fnite
// Difference Method for Stratigraphic, Hexahedron Grids," by
// S. H. Lee, H. Tchelepi, and L. J. DeChant, \emph{1999 SPE Reservoir
// Simulation Symposium}, SPE (Society of Petroleum Engineers) 51901.

// Preprocessor symbols:
// PSEUDOCODE: Do not define this symbol.  Surrounds desired code to
//	       deal with different granularity fields.
// DEBUG: If defined, print some information about internal program
//	  values.


template<class GeometryTag, class T, class Expr, int Dim>
inline
typename Field<GeometryTag, T, Expr>::T_t
faceWeightedSum(const Field<GeometryTag, T, Expr>& inputField, 
		const FieldOffsetList<Dim,Type> &lst,
		const Field<GeometryTag, T, Expr>& outputField)
{
  typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
  typedef typename FieldOffsetList<Dim,Type>::size_type size_type;
  CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));

  
	
	      

  // HERE
  
  const size_type lstLength = lst.size();
  PInsist(lstLength > 0, "faceWeightedSum() must be given a nonempty list.");
  T_t init = inputField(lst[0], loc);

  // FIXME inputField.mesh().face(arg).normal() returns a normal
  // vector with length equal to the face's area and direction
  // perpendicular to the face.

  for (size_type i = 1; i < lstLength ; ++i)
    init +=
      outputField.mesh().face(WHICH).normal() *
      inputField(lst[i], loc);
  return init;
}


/** THE PROGRAM **/

int main(int argc, char *argv[])
{
  // Set up the Pooma library.
  Pooma::initialize(argc,argv);
#ifdef DEBUG
  std::cout << "Start program." << std::endl;
#endif // DEBUG

  /* DECLARATIONS */

  // Create a simple layout.
  const unsigned Dim = 2;		// Work in a 2D world.
  const unsigned nXs = 5;		// number of horizontal vertices
  const unsigned nYs = 4;		// number of vertical vertices
  Interval<Dim> meshDomain;
  meshDomain[0] = Interval<1>(nXs);
  meshDomain[1] = Interval<1>(nYs);
  DomainLayout<Dim> meshLayout(meshDomain, GuardLayers<Dim>(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;
  typedef Field<Geometry_t, double, Brick> ConstFields_t; // TODO: Change to ConstField when ConstField is available.
  typedef Tensor<Dim,double,Full> Tensor_t;
  typedef Field<Geometry_t, Tensor_t, Brick> FieldT_t;
  typedef Field<Geometry_t, Tensor_t, Brick> ConstFieldT_t; // TODO: Change to ConstField when ConstField is available.
  typedef Field<Geometry_t, Vector<Dim>, Brick> Fieldv_t;
  typedef Field<Geometry_t, Vector<Dim>, Brick> ConstFieldv_t; // TODO: Change to ConstField when ConstField is available.

  // Cell-centered Fields.

  Centering<Dim> cell = canonicalCentering<Dim>(CellType, Continuous);
  ConstFieldT_t permeability	(cell, meshLayout, origin, spacings);
  ConstFields_t pressure	(cell, meshLayout, origin, spacings);
  Fields_t totalFlux		(cell, meshLayout, origin, spacings);

  // Subcell-centered Field.

  typedef Centering<Dim>::Orientation Orientation;
  typedef Centering<Dim>::Position Position;
  Position position;
  Centering<Dim> subcell(CellType, Continuous);
  position(0) = 0.25;
  position(1) = 0.25; subcell.addValue(Orientation(1), position);
  position(1) = 0.75; subcell.addValue(Orientation(1), position);
  position(0) = 0.75; subcell.addValue(Orientation(1), position);
  position(1) = 0.25; subcell.addValue(Orientation(1), position);
  Fields_t pressureGradient	(subcell, meshLayout, origin, spacings);

  // Spoke-centered Field.

  Centering<Dim> spoke(FaceType, Discontinuous);
  // QUESTION: These are supposed to be Discontinuous, right?
  Orientation orientation;
  // NOTE: This code is not dimension-independent.
  for (int zeroFace = 0; zeroFace < 2; ++zeroFace) {
    orientation = 1; orientation[zeroFace] = 0;
    position(zeroFace) = 0.0;
    position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
    position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
    position(zeroFace) = 1.0;
    position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
    position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
  }
  Fields_t spokeFlux		(spoke, meshLayout, origin, spacings);

  // Face-centered.

  Centering<Dim> disFace = canonicalCentering<Dim>(FaceType, Discontinuous);


  /* INITIALIZATION */

#ifdef PSEUDOCODE
  // Initialize tensors.
  // Initialize the pressures.
  // Initialize coordinates.
#endif // PSEUDOCODE


  /* COMPUTATION */

#ifdef PSEUDOCODE
  // Compute pressureGradients by simultaneously solving several
  // linear equations.  The operands have different centerings.
  // FIXME
  pressureGradients =
    linearAlgebra<2>(pressure /* cell-centered */,
		     /* Interpolate from vertex-centered to cell-centered: */
		     interpolate<Cell,Vertex>(coordinates),
		     permeability /* cell-centered */,
		     normals /* face-centered */);
#endif // PSEUDOCODE

  // Compute the spoke fluxes.

  // We must multiply three quantities, each with a different
  // centering, to yield values at a fourth-centering.  permeability
  // is cell-centered.  pressureGradient is subcell-centered.  The
  // normals are face-centered.  The product is spoke-centered.

  spokeFlux = 
    dot(replicate(dot(replicate(permeability, nearestNeighbors(cell, subcell)),
		      pressureGradient),
		  nearestNeighbors(subcell, spoke)),
	replicate(meshLayout.unitCoordinateNormals(),
		  nearestNeighbors(disFace, spoke)));

  // Sum the spoke fluxes into a cell flux.

  // Q = \sum_{faces f of cell} sign_f area_f \sum_{subfaces sf of f} q_{sf}
  // We compute this in three steps:
  // 1. Add together the flux values on each face to form a
  //    face-centered field.
  // 2. Multiply each value by the magnitude of the face's normal.
  // 3. Add together each face's value.

  totalFlux =
    sum(spokeFlux.mesh().normals().signedMagnitude() *
	sum(spokeFlux, nearestNeighbors(spoke, disFace)),
	findFieldOffsetList(disFace, cell));


  /* TERMINATION */

  std::cout << "total flux:\n" << totalFlux << std::endl;
#ifdef DEBUG
  std::cout << "End program." << std::endl;
#endif // DEBUG
  Pooma::finalize();
  return EXIT_SUCCESS;
}