[PATCH] Add RectilinearMesh
Jeffrey D. Oldham
oldham at codesourcery.com
Mon Jul 12 18:42:46 UTC 2004
Richard Guenther wrote:
> 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?
Yes. I wonder why I am listed as the last author to touch the file
since I do not remember modifying it. It's been a long time.
> Richard.
>
>
> 2004Jul11 Richard Guenther <richard.guenther at uni-tuebingen.de>
>
> * 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 pooma at acl.lanl.gov,
>+// 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
>+
>
>
--
Jeffrey D. Oldham
oldham at codesourcery.com
More information about the pooma-dev
mailing list