[PATCH] Add RectilinearMesh

[Richard Guenther]

Richard Guenther wrote:
> This moves the particle tests to new Field.  RectilinearMesh still
> missing for 2 and 4, though.  Hm, I remembered submitting that long time
> ago.

Indeed I have.  This is mainly RectilinearMesh from v2.1 fixed for v2.4.
particle_test[24] will then work, too (after some minor fix).

Ok?

Richard.


2004Jul11  Richard Guenther <address@hidden>

        * src/Field/Mesh/RectilinearMesh.h: new.
        src/Pooma/Fields.h: include RectilinearMesh.h.

Index: Pooma/Fields.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/Pooma/Fields.h,v
retrieving revision 1.16
diff -u -u -r1.16 Fields.h
--- Pooma/Fields.h      21 Nov 2003 17:36:10 -0000      1.16
+++ Pooma/Fields.h      11 Jul 2004 17:05:25 -0000
@@ -55,6 +55,7 @@
 
 #include "Field/Mesh/NoMesh.h"
 #include "Field/Mesh/UniformRectilinearMesh.h"
+#include "Field/Mesh/RectilinearMesh.h"
 #include "Field/Mesh/MeshFunctions.h"
 #include "Field/Mesh/PositionFunctions.h"
 
Index: Field/Mesh/RectilinearMesh.h
===================================================================
RCS file: Field/Mesh/RectilinearMesh.h
diff -N Field/Mesh/RectilinearMesh.h
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ Field/Mesh/RectilinearMesh.h        11 Jul 2004 17:05:27 -0000
@@ -0,0 +1,904 @@
+// -*- C++ -*-
+// ACL:license
+// ----------------------------------------------------------------------
+// This software and ancillary information (herein called "SOFTWARE")
+// called POOMA (Parallel Object-Oriented Methods and Applications) is
+// made available under the terms described here.  The SOFTWARE has been
+// approved for release with associated LA-CC Number LA-CC-98-65.
+// 
+// Unless otherwise indicated, this SOFTWARE has been authored by an
+// employee or employees of the University of California, operator of the
+// Los Alamos National Laboratory under Contract No. W-7405-ENG-36 with
+// the U.S. Department of Energy.  The U.S. Government has rights to use,
+// reproduce, and distribute this SOFTWARE. The public may copy, distribute,
+// prepare derivative works and publicly display this SOFTWARE without 
+// charge, provided that this Notice and any statement of authorship are 
+// reproduced on all copies.  Neither the Government nor the University 
+// makes any warranty, express or implied, or assumes any liability or 
+// responsibility for the use of this SOFTWARE.
+// 
+// If SOFTWARE is modified to produce derivative works, such modified
+// SOFTWARE should be clearly marked, so as not to confuse it with the
+// version available from LANL.
+// 
+// For more information about POOMA, send e-mail to address@hidden,
+// or visit the POOMA web page at http://www.acl.lanl.gov/pooma/.
+// ----------------------------------------------------------------------
+// ACL:license
+
+/** @file
+ * @ingroup Mesh
+ * @brief
+ * A rectilinear mesh without uniform spacing between vertices.
+ */
+
+#ifndef POOMA_FIELD_MESH_RECTILINEARMESH_H
+#define POOMA_FIELD_MESH_RECTILINEARMESH_H
+
+//-----------------------------------------------------------------------------
+// Includes:
+//-----------------------------------------------------------------------------
+
+#include <iostream>
+
+#include "Engine/ConstantFunctionEngine.h"         // Used in functors
+#include "Engine/IndexFunctionEngine.h"            // Used in functors
+#include "Layout/INode.h"                          // Used in ctors
+#include "Field/FieldEngine/FieldEnginePatch.h" // Used in ctors
+#include "Field/Mesh/NoMesh.h"                  // Base class
+#include "Field/FieldCentering.h"               // Centering<Dim> inline
+#include "Tiny/Vector.h"                        // Class member
+#include "Field/Mesh/MeshTraits.h"              // Template parameter
+
+
+//-----------------------------------------------------------------------------
+/// Holds the data for a rectilinear mesh. That class has a ref-counted
+/// instance of this class.
+//-----------------------------------------------------------------------------
+
+template <class MeshTraits>
+class RectilinearMeshData : public NoMeshData<MeshTraits::dimensions>
+{
+public:
+
+  // shortcuts to MeshTraits types
+  typedef typename MeshTraits::Domain_t Domain_t;
+  typedef typename MeshTraits::MeshData_t MeshData_t;
+  typedef typename MeshTraits::Scalar_t Scalar_t;
+  typedef typename MeshTraits::PointType_t PointType_t;
+  typedef typename MeshTraits::VectorType_t VectorType_t;
+  typedef typename MeshTraits::SpacingsType_t SpacingsType_t;
+  typedef typename MeshTraits::PositionsType_t PositionsType_t;
+
+  enum { dimensions = MeshTraits::dimensions };
+
+
+  //---------------------------------------------------------------------------
+  // Constructors.
+
+  /// We provide a default constructor that creates the object with empty
+  /// domains. To be useful, this object must be replaced by another 
+  /// version via assignment.
+  
+  RectilinearMeshData()
+    { 
+      // This space intentionally left blank.
+    }
+
+  /// This constructor fully constructs the object. It uses the layout to
+  /// compute domains and also initializes the origin and the spacings in each
+  /// coordinate direction. The indices in the layout refer to VERTEX
+  /// positions.
+
+  template<class Layout, class EngineTag>
+  RectilinearMeshData(
+    const Layout &layout,
+    const PointType_t &origin,
+    const Vector<dimensions, Array<1, Scalar_t, EngineTag> > &spacings)
+  : NoMeshData<dimensions>(layout), 
+    origin_m(origin)
+    //spacings_m(spacings)
+    {
+      for (int i=0; i<dimensions; i++) {
+       spacings_m(i).engine() = spacings(i).engine(); // init
+       spacings_m(i).engine().makeOwnCopy(); // FIXME? Do we want this?
+       Interval<1> I(layout.domain()[i]);
+       positions_m(i).engine() = Engine<1, Scalar_t, Brick>(I);
+       positions_m(i)(0) = origin_m(i);
+       // initialize from origin downward the ghost cells
+       for (int j=-1; j>=I.min(); j--)
+         positions_m(i)(j) = positions_m(i).read(j+1) - spacings_m(i).read(j);
+       // initialize from origin upward
+       for (int j=1; j<=I.max(); j++)
+         positions_m(i)(j) = positions_m(i).read(j-1) + 
spacings_m(i).read(j-1);
+      }
+    }
+
+  /// Constructor for constructing evenly spaced rectilinear meshes just
+  /// like UniformRectilinearMesh does.
+
+  template<class Layout>
+  RectilinearMeshData(
+    const Layout &layout,
+    const PointType_t &origin,
+    const VectorType_t &spacings)
+  : NoMeshData<dimensions>(layout), 
+    origin_m(origin)
+    {
+      // for each dimension we allocate engines for spacings & positions
+      // and initialize them according to origin/spacings
+      for (int i=0; i<dimensions; i++) {
+       Interval<1> I(layout.domain()[i]);
+       // allocate and assign spacings
+       spacings_m(i).engine() = Engine<1, Scalar_t, Brick>(I);
+       spacings_m(i)(I) = spacings(i); // no Array.all()
+       Pooma::blockAndEvaluate();
+       // allocate positions, assign origin
+       positions_m(i).engine() = Engine<1, Scalar_t, Brick>(I);
+       positions_m(i)(0) = origin_m(i);
+       // initialize from origin downward the ghost cells
+       for (int j=-1; j>=I.min(); j--)
+         positions_m(i)(j) = positions_m(i).read(j+1) - spacings_m(i).read(j);
+       // initialize from origin upward
+       for (int j=1; j<=I.max(); j++)
+         positions_m(i)(j) = positions_m(i).read(j-1) + 
spacings_m(i).read(j-1);
+      }
+    }
+    
+  /// Copy constructor.
+
+  RectilinearMeshData(const MeshData_t &model)
+  : NoMeshData<dimensions>(model), 
+    origin_m(model.origin_m)
+    //spacings_m(model.spacings_m),
+    //positions_m(model.positions_m)
+    {
+      for (int i=0; i<dimensions; i++) {
+       spacings_m(i).engine() = model.spacings_m(i).engine();
+       positions_m(i).engine() = model.positions_m(i).engine();
+      }
+      // This space intentionally left blank.
+    } 
+    
+  /// @name View constructors.
+  //@{
+  
+  /// Interval view. This means that we simply need to adjust the
+  /// origin by the amount the view is offset from the model's physical
+  /// cell domain. We rely on the base class to do the heavy lifting
+  /// with respect to figuring out the domains correctly.
+  ///
+  /// The Interval supplied must refer to VERTEX positions.
+  
+  RectilinearMeshData(const MeshData_t &model, 
+                     const Interval<dimensions> &d)
+  : NoMeshData<dimensions>(d)
+    {
+      for (int i = 0; i < dimensions; i++) {
+       // FIXME: Wheeee ;) (we cant store a BrickView...
+       // and still dont want to copy)
+       spacings_m(i).engine() = Engine<1, Scalar_t, 
Brick>(&model.spacings_m(i)(d[i])(0), d[i]);
+       positions_m(i).engine() = Engine<1, Scalar_t, 
Brick>(&model.positions_m(i)(d[i])(0), d[i]);
+       origin_m(i) = positions_m(i)(d[i].min());
+      }
+    }
+#if 0  
+  /// FieldEnginePatch view. We don't fiddle with the origin because we are not
+  /// making the domain zero-based.
+  ///
+  /// The domain supplied by the FieldEnginePatch must refer to VERTEX
+  /// positions.
+  
+  RectilinearMeshData(const MeshData_t &model, 
+                     const FieldEnginePatch<dimensions> &p)
+  : NoMeshData<dimensions>(model, p),
+    origin_m(model.origin_m),
+    spacings_m(model.spacings_m),
+    positions_m(model.spacings_m)
+    {
+      std::cerr << "RectilinearMeshData(FieldEnginePatch) constructor called" 
<< std::endl;
+      abort();
+      // FIXME: what does FieldEnginePatch do???
+      for (int i=0; i<dimensions; i++) {
+       spacings_m(i).engine() = model.spacings_m(i).engine();
+       positions_m(i).engine() = model.positions_m(i).engine();
+      }
+    }
+#endif
+  //@}
+
+  //---------------------------------------------------------------------------
+  /// Copy assignment operator.
+  
+  MeshData_t &
+  operator=(const MeshData_t &rhs)
+    {
+      if (this != &rhs)
+        {
+          NoMeshData<dimensions>::operator=(rhs);
+          origin_m = rhs.origin_m;
+         for (int i=0; i<dimensions; i++) {
+           spacings_m(i).engine() = rhs.spacings_m(i).engine();
+           positions_m(i).engine() = rhs.positions_m(i).engine();
+         }
+        }
+        
+      return *this;
+    }
+
+  //---------------------------------------------------------------------------
+  /// Empty destructor is fine. Note, however, that NoMeshData does not have
+  /// a virtual destructor. We must be careful to delete these puppies as
+  /// RectilinearMeshData.
+
+  ~RectilinearMeshData() { }
+
+  //---------------------------------------------------------------------------
+  /// @name General accessors.
+  //@{
+
+  /// The mesh spacing.
+  
+  inline const SpacingsType_t &spacings() const 
+    { 
+      return spacings_m; 
+    }
+
+  /// The mesh vertex positions.
+  
+  inline const PositionsType_t &positions() const 
+    { 
+      return positions_m; 
+    }
+
+  /// The mesh origin.
+
+  inline const PointType_t &origin() const 
+    { 
+      return origin_m; 
+    }
+
+  //@}
+
+private:
+
+  /// Origin of mesh (coordinate vector of first vertex).
+
+  PointType_t origin_m;
+
+  /// Spacings between vertices.
+
+  SpacingsType_t spacings_m;
+
+  /// Vertex positions.
+
+  PositionsType_t positions_m;
+
+};
+
+
+
+
+///
+/// RectilinearMesh is a rectilinear mesh sometimes called a 
+/// "cartesian product" or "tensor product" mesh. Each dimension has a
+/// spacing value between every pair of vertices along that dimension;
+/// these spacings can all be different.
+///
+template<class MeshTraits>
+class RectilinearMesh : public MeshTraits::CoordinateSystem_t
+{
+public:
+
+  //---------------------------------------------------------------------------
+  // Exported typedefs and enumerations.
+
+  typedef MeshTraits MeshTraits_t;
+
+  /// The type of the mesh class.
+
+  typedef typename MeshTraits::Mesh_t Mesh_t;
+
+  /// The type of the mesh data class.
+
+  typedef typename MeshTraits::MeshData_t MeshData_t;
+
+  /// The type of the coordinate system.
+
+  typedef typename MeshTraits::CoordinateSystem_t CoordinateSystem_t;
+
+  /// The type of domains.
+
+  typedef typename MeshTraits::Domain_t Domain_t;
+
+  /// The type of locations.
+
+  typedef typename MeshTraits::Loc_t Loc_t;
+
+  /// The type T, used to represent, for example, volumes & areas, etc.
+
+  typedef typename MeshTraits::T_t T_t;
+
+  /// The type of scalars.
+
+  typedef typename MeshTraits::Scalar_t Scalar_t;
+
+  /// The type of mesh points.
+    
+  typedef typename MeshTraits::PointType_t PointType_t;
+
+  /// The type of vectors used to represent, for example, normals.
+  
+  typedef typename MeshTraits::VectorType_t VectorType_t;
+
+  /// The type used to store spacings.
+
+  typedef typename MeshTraits::SpacingsType_t SpacingsType_t;
+
+  /// The type used to store positions.
+
+  typedef typename MeshTraits::PositionsType_t PositionsType_t;
+
+  /// The number of indices required to select a point in this mesh.
+
+  enum { dimensions = MeshTraits::dimensions };
+
+  /// The number of components of a position vector inside the mesh.
+
+  enum { coordinateDimensions = MeshTraits::coordinateDimensions };
+
+
+  //---------------------------------------------------------------------------
+  // Constructors.
+  
+  /// We supply a default constructor, but it doesn't generate a useful mesh.
+  /// This is accomplished through assignment.
+  
+  RectilinearMesh() 
+  : data_m(new MeshData_t)
+    { 
+      // This space intentionally left blank.
+    }
+
+  /// This constructor fully constructs the object. It uses the layout to
+  /// compute domains and also initializes the origin and the spacings in each
+  /// coordinate direction.
+  ///
+  /// The Layout supplied must refer to VERTEX positions.
+  
+  template<class Layout, class EngineTag>
+  inline RectilinearMesh(const Layout &layout, 
+                        const PointType_t &origin,
+                        const Vector<coordinateDimensions, Array<1, T_t, 
EngineTag> > &spacings)
+  : data_m(new MeshData_t(layout, origin, spacings))
+    { 
+    }
+
+  /// Constructor compatible to UniformRectilinearMesh.
+
+  template<class Layout>
+  inline RectilinearMesh(const Layout &layout,
+                        const PointType_t &origin,
+                        const PointType_t &spacings)
+  : data_m(new MeshData_t(layout, origin, spacings))
+    { 
+    }
+
+  template<class Layout>
+  inline explicit RectilinearMesh(const Layout &layout)
+  : data_m(new MeshData_t(layout, 
+                                          PointType_t(0), 
+                                          PointType_t(1)))
+    { 
+    }
+
+  /// Copy constructor. 
+  
+  inline RectilinearMesh(const Mesh_t &model)
+  : data_m(model.data_m)
+    {
+    }
+    
+  /// @name View constructors
+  /// These are the only possible views of this
+  /// mesh. Other views will make a NoMesh.
+  //@{ 
+  
+  /// Interval view.
+  ///
+  /// The Interval supplied must refer to VERTEX positions.
+  
+  inline RectilinearMesh(const Mesh_t &model, 
+                        const Domain_t &d)
+  : data_m(new MeshData_t(*model.data_m, d))
+    {
+    }
+  
+  /// INode view.
+  ///
+  /// The INode supplied must refer to VERTEX positions.
+  
+  inline RectilinearMesh(const Mesh_t &model, 
+                        const INode<dimensions> &i)
+  : data_m(new MeshData_t(*model.data_m, i.domain()))
+    {
+    }
+#if 0
+  /// FieldEnginePatch view.
+  ///
+  /// The FieldEnginePatch supplied must refer to VERTEX positions.
+  
+  inline RectilinearMesh(const Mesh_t &model, 
+                        const FieldEnginePatch<dimensions> &p)
+  : data_m(new MeshData_t(*model.data_m, p))
+    {
+    }
+#endif
+  //@}
+
+  //---------------------------------------------------------------------------
+  /// Copy assignment operator.
+  
+  inline Mesh_t &
+  operator=(const Mesh_t &rhs)
+    {
+      if (&rhs != this)
+        {
+          data_m = rhs.data_m;
+        }
+      
+      return *this;
+    }
+
+  //---------------------------------------------------------------------------
+  /// Empty destructor is fine. The pointer to the data is ref-counted so its
+  /// lifetime is correctly managed.
+  
+  ~RectilinearMesh()
+    {
+    }
+
+  /// Data member access.
+  const MeshData_t& data() const
+    {
+      return *data_m;
+    }
+  
+  //---------------------------------------------------------------------------
+  /// @name Domain functions.
+  //@{
+  
+  /// The vertex domain, as the mesh was constructed with.
+
+  inline const Domain_t &physicalVertexDomain() const
+    {
+      return data_m->physicalVertexDomain(); 
+    }
+
+  /// Function that returns a domain adjusted to give the indices of the cells.
+
+  inline const Domain_t &physicalCellDomain() const
+    {
+      return data_m->physicalCellDomain(); 
+    }
+
+  /// The total vertex domain, including mesh guard vertices.
+
+  inline const Domain_t &totalVertexDomain() const
+    {
+      return data_m->totalVertexDomain(); 
+    }
+
+  /// The total cell domain, including mesh guard cells.
+
+  inline const Domain_t &totalCellDomain() const
+    {
+      return data_m->totalCellDomain(); 
+    }
+
+  //@}
+
+  //---------------------------------------------------------------------------
+  /// @name General accessors.
+  //@{
+
+  /// The mesh origin.
+
+  inline const PointType_t &origin() const 
+    { 
+      return data_m->origin();
+    }
+
+  /// The mesh spacing. Return type is dependend on mesh type.
+
+  inline const SpacingsType_t &spacings() const 
+    { 
+      return data_m->spacings();
+    }
+
+  /// The mesh positions. Return type is dependend on mesh type.
+
+  inline const PositionsType_t &positions() const 
+    { 
+      return data_m->positions();
+    }
+
+  /// The cell containing a particular point.
+
+  inline Loc_t cellContaining(const PointType_t &point) const
+    {
+      /// FIXME
+      Loc_t loc((0, Pooma::NoInit())); // Avoid a g++ parse error.
+      for (int i = 0; i < dimensions; i++)
+       {
+         const T_t *start = &positions()(i)(0);
+         const T_t *finish = start + 
positions()(i).physicalDomain()[i].length();
+         const T_t *p = std::lower_bound(start, finish, point(i));
+#if POOMA_BOUNDS_CHECK
+         PInsist(p != finish,
+                 "Rectilinear::cellContaining(): point is outside mesh.");
+#endif
+         // The lower_bound function returns the first element that is not
+         // less than the point we're searching for.
+         int j = static_cast<int>(std::distance(start, p));
+         if (*p == point(i))
+           loc[i] = j;
+         else
+           loc[i] = j-1;
+       }
+
+      return loc;
+    }
+
+  /// The lower-left vertex associated with a given cell location.
+    
+  inline PointType_t vertexPosition(const Loc_t &loc) const
+    {
+      PointType_t point;
+      for (int i = 0; i < dimensions; i++)
+        point(i) = positions()(i)(loc[i]); 
+      return point;
+    }
+
+  inline Scalar_t vertexPosition(int dim, int i) const
+    {
+      return positions()(dim)(i);
+    }
+
+  /// The position of a given cell location for canonical cell centering.
+
+  inline PointType_t cellPosition(const Loc_t &loc) const
+    {
+      PointType_t point;
+
+      for (int i=0; i<dimensions; i++)
+       point(i) = positions()(i)(loc[i]) + 0.5*spacings()(i)(loc[i]);
+
+      return point;
+    }
+
+  inline Scalar_t cellPosition(int dim, int i) const
+    {
+      return positions()(dim)(i) + 0.5*spacings()(dim)(i);
+    }
+
+  /// The vertex spacing for a given cell location.
+
+  inline VectorType_t vertexSpacing(const Loc_t &loc) const
+    {
+      VectorType_t delta;
+
+      for (int i=0; i<dimensions; i++)
+       delta(i) = spacings()(i)(loc[i]);
+
+      return delta;
+    }
+
+  inline Scalar_t vertexSpacing(int dim, int i) const
+    {
+      return spacings()(dim)(i);
+    }
+
+  /// The cell spacing for a given cell location.
+
+  inline VectorType_t cellSpacing(const Loc_t &loc) const
+    {
+      VectorType_t delta;
+
+      for (int i=0; i<dimensions; i++)
+       delta(i) = 0.5 * (spacings()(i)(loc[i]) + spacings()(i)(loc[i]+1));
+
+      return delta;
+    }
+
+  inline Scalar_t cellSpacing(int dim, int i) const
+    {
+      return 0.5*(spacings()(dim)(i) + spacings()(dim)(i+1));
+    }
+
+  //@}
+
+
+private:
+
+  /// Our data, stored as a ref-counted pointer to simplify memory management.
+  
+  RefCountedPtr<MeshData_t> data_m;
+};
+
+
+
+///
+/// GenericRM contains mesh functions related functors that are applicapble
+/// regardless of the coordinate system type.
+///
+template <class MeshTraits>
+struct GenericRM {
+
+  /// The coordinate type.
+
+  typedef typename MeshTraits::T_t T_t;
+
+  /// The mesh data class.
+
+  typedef typename MeshTraits::Mesh_t Mesh_t;
+
+  /// The type to represent points.
+
+  typedef typename MeshTraits::PointType_t PointType_t;
+
+  /// The type to represent vectors.
+
+  typedef typename MeshTraits::VectorType_t VectorType_t;
+
+  /// The type to represent the spacings.
+
+  typedef typename MeshTraits::SpacingsType_t SpacingsType_t;
+
+  /// The type used to store positions.
+
+  typedef typename MeshTraits::PositionsType_t PositionsType_t;
+
+  /// The type of locations.
+
+  typedef typename MeshTraits::Loc_t Loc_t;
+
+  /// Dimensionality of the mesh.
+
+  enum { dimensions = MeshTraits::dimensions };
+  enum { coordinateDimensions = MeshTraits::coordinateDimensions };
+
+
+  //---------------------------------------------------------------------------
+  /// Support for the positions() function. We need to provide a functor for
+  /// use with IndexFunction-engine. We also need to export the
+  /// PositionsEngineTag_t typedef and the positionsFunctor() member function,
+  /// which computes the positions using the centering point positions.
+  /// The indices passed in refer to cells.
+  
+  class PositionsFunctor {
+  public:
+  
+    /// Need to be able to default construct since we fill in the details
+    /// after the fact.
+
+    // WARNING! For Arrays to be initialized (copy constructed, assigned,
+    //          etc.) correctly, even in the case of uninitialized targets
+    //          we need to copy the engines explicitly rather than rely
+    //          on the compiler generating correct copy constructors and
+    //          assignment operators.
+    // FIXME! Technically we either can dump the copy constructor or the
+    //        assignment operator.
+
+    PositionsFunctor() { }
+    
+    PositionsFunctor(const Mesh_t &m, 
+                     const Centering<dimensions> &c)
+      : centering_m(c.position(0))
+      {
+       for (int i=0; i<dimensions; i++) {
+         positions_m(i).engine() = m.positions()(i).engine();
+         spacings_m(i).engine() = m.spacings()(i).engine();
+       }
+      }
+
+    PositionsFunctor(const PositionsFunctor &m)
+      :        centering_m(m.centering_m)
+    {
+      for (int i=0; i<dimensions; i++) {
+       positions_m(i).engine() = m.positions_m(i).engine();
+       spacings_m(i).engine() = m.spacings_m(i).engine();
+      }
+    }
+
+    PositionsFunctor& operator=(const PositionsFunctor &m)
+    {
+      centering_m = m.centering_m;
+      for (int i=0; i<dimensions; i++) {
+       positions_m(i).engine() = m.positions_m(i).engine();
+       spacings_m(i).engine() = m.spacings_m(i).engine();
+      }
+
+      return *this;
+    }
+
+    inline PointType_t operator()(int i0) const
+      {
+        return PointType_t(positions_m(0).read(i0) + 
spacings_m(0).read(i0)*centering_m(0));
+      }
+      
+    inline PointType_t operator()(int i0, int i1) const
+      {
+        return PointType_t(positions_m(0).read(i0) + 
spacings_m(0).read(i0)*centering_m(0),
+                          positions_m(1).read(i1) + 
spacings_m(1).read(i1)*centering_m(1));
+      }
+
+    inline PointType_t operator()(int i0, int i1, int i2) const
+      {
+        return PointType_t(positions_m(0).read(i0) + 
spacings_m(0).read(i0)*centering_m(0),
+                          positions_m(1).read(i1) + 
spacings_m(1).read(i1)*centering_m(1),
+                          positions_m(2).read(i2) + 
spacings_m(2).read(i2)*centering_m(2));
+      }
+
+  private:
+
+    PositionsType_t positions_m;
+    SpacingsType_t spacings_m;
+    typename Centering<dimensions>::Position centering_m;
+
+  };
+  
+  typedef IndexFunction<PositionsFunctor> PositionsEngineTag_t;
+
+  template <class PositionsEngineTag>
+  void initializePositions(
+    Engine<dimensions, PointType_t, PositionsEngineTag> &e, 
+    const Centering<dimensions> &c) const
+    {
+      e.setFunctor(typename PositionsEngineTag::Functor_t(static_cast<const 
Mesh_t&>(*this), c));
+    }
+
+
+  //---------------------------------------------------------------------------
+  /// Support for the spacings() function. We need to provide a functor for
+  /// use with IndexFunction-engine. We also need to export the
+  /// SpacingsEngineTag_t typedef and the spacingsFunctor() member function,
+  /// which computes the spacings using the centering point positions.
+  /// The indices passed in refer to cells.
+  
+  class SpacingsFunctor {
+  public:
+  
+    /// Need to be able to default construct since we fill in the details
+    /// after the fact.
+
+    // WARNING! For Arrays to be initialized (copy constructed, assigned,
+    //          etc.) correctly, even in the case of uninitialized targets
+    //          we need to copy the engines explicitly rather than rely
+    //          on the compiler generating correct copy constructors and
+    //          assignment operators.
+    // FIXME! Technically we either can dump the copy constructor or the
+    //        assignment operator.
+
+    SpacingsFunctor() { }
+    
+    SpacingsFunctor(const Mesh_t &m, 
+                   const Centering<dimensions> &c)
+      : centering_m(c.position(0))
+      {
+       for (int i=0; i<dimensions; i++)
+         spacings_m(i).engine() = m.spacings()(i).engine();
+      }
+
+    SpacingsFunctor(const SpacingsFunctor &m)
+      :        centering_m(m.centering_m)
+    {
+      for (int i=0; i<dimensions; i++)
+       spacings_m(i).engine() = m.spacings_m(i).engine();
+    }
+
+    SpacingsFunctor& operator=(const SpacingsFunctor &m)
+    {
+      centering_m = m.centering_m;
+      for (int i=0; i<dimensions; i++)
+       spacings_m(i).engine() = m.spacings_m(i).engine();
+
+      return *this;
+    }
+
+    /* FIXME: the following may cause an out of bound condition, if
+     *        the spacings is queried for the last cell - spacings
+     *        for non-existing cells are used then.
+     */
+
+    inline VectorType_t operator()(int i0) const
+      {
+        return VectorType_t(spacings_m(0).read(i0)
+                          + 
(spacings_m(0).read(i0+1)-spacings_m(0).read(i0))*centering_m(0));
+      }
+      
+    inline VectorType_t operator()(int i0, int i1) const
+      {
+        return VectorType_t(spacings_m(0).read(i0)
+                          + 
(spacings_m(0).read(i0+1)-spacings_m(0).read(i0))*centering_m(0),
+                          spacings_m(1).read(i1)
+                          + 
(spacings_m(1).read(i1+1)-spacings_m(1).read(i1))*centering_m(1));
+      }
+
+    inline VectorType_t operator()(int i0, int i1, int i2) const
+      {
+        return VectorType_t(spacings_m(0).read(i0)
+                          + 
(spacings_m(0).read(i0+1)-spacings_m(0).read(i0))*centering_m(0),
+                          spacings_m(1).read(i1)
+                          + 
(spacings_m(1).read(i1+1)-spacings_m(1).read(i1))*centering_m(1),
+                          spacings_m(2).read(i2)
+                          + 
(spacings_m(2).read(i2+1)-spacings_m(2).read(i2))*centering_m(2));
+      }
+
+  private:
+
+    SpacingsType_t spacings_m;
+    typename Centering<dimensions>::Position centering_m;
+
+  };
+  
+  typedef IndexFunction<SpacingsFunctor> SpacingsEngineTag_t;
+  
+  void initializeSpacings(
+    Engine<dimensions, PointType_t, SpacingsEngineTag_t> &e, 
+    const Centering<dimensions> &c) const
+    {
+      e.setFunctor(SpacingsFunctor(static_cast<const Mesh_t&>(*this), c));
+    }
+
+  
+  //---------------------------------------------------------------------------
+  /// Support for the outwardNormals() and coordinateNormals() functions. 
+  /// We also need to export the NormalsEngineTag_t typedef and the 
+  /// initializeNormals() member function, which sets the appropriate constant 
+  /// value (since the normals exactly align with the coordinate axes).
+  /// The boolean value passed is true if we are asking for outward normals,
+  /// as opposed to coordinate normals. The indices passed in refer to cells.
+
+  typedef ConstantFunction NormalsEngineTag_t;
+  
+  void initializeNormals(
+    Engine<dimensions, VectorType_t, NormalsEngineTag_t> &e, 
+    const Centering<dimensions> &c,
+    bool outward = true) const
+    {
+      // Check some pre-conditions. We need there to be a single centering
+      // point and it must be face-centered.
+      
+      PAssert(c.size() == 1);
+      PAssert(c.centeringType() == FaceType);
+      
+      // Generate the normals. The coordinate normals are computed from
+      // 1 - orientation. Then, if we are on the near face, indicated by
+      // position == 0.0, we need to multiply by -1.0 if we are doing
+      // outward normals.
+      
+      VectorType_t normal;
+      for (int i = 0; i < dimensions; i++)
+        {
+          normal(i) = static_cast<T_t>(1 - c.orientation(0)[i].first());
+          if (outward && c.position(0)(i) == 0.0)
+            normal(i) *= static_cast<T_t>(-1);
+        }
+        
+      e.setConstant(normal);
+    }
+
+};
+
+
+
+#endif // POOMA_FIELD_MESH_RECTILINEARMESH_H
+
+// ACL:rcsinfo
+// ----------------------------------------------------------------------
+// $RCSfile: RectilinearMesh.h,v $   $Author: oldham $
+// $Revision: 1.4 $   $Date: 2001/12/11 20:43:30 $
+// ----------------------------------------------------------------------
+// ACL:rcsinfo
+