[newfield_revision] Nearest Neighbors, Reduction Functions
Jeffrey Oldham
oldham at codesourcery.com
Thu Aug 2 22:51:39 UTC 2001
The attached patch revises FieldOffset, eliminating the Type template
parameter and adding miscellaneous functionality such as operator<<
and operator==. The eliminated Type parameter indicated whether the
underlying field had one subfield or multiple subfields, which are
accessed differently. Stephen, Scott, and I all repeatedly tripped
over the template parameter, which propagated itself through much of
the code changes. Instead, of a compile-time parameter, we use a
run-time computation to decide how to access a field. The associated
test program is also revised.
More importantly, this patch also adds FieldOffsetList,
nearestNeighbors(), NearestNeighbors, associated reduction functions,
and associated test programs FieldReductions and NearestNeighbors. A
FieldOffsetList is a list of FieldOffsets. Given input and output
centerings, nearestNeighbors() returns a std::vector of
FieldOffsetLists. Each FieldOffsetList specifies the locations of the
input field values closest (using the Manhattan norm) to the
corresponding output field value. Many of the Caramana hydrodynamics
and Chevron kernel operations use nearest neighbors since physics
operates locally. FieldReductions, taking an input field, a
FieldOffsetList, and an output cell location, apply a reduction to the
specified input cells.
All these changes were added to the newfield_revision branch, which is
being developed by Scott Haney, Jeffrey D. Oldham, and Stephen Smith.
The branch code is not guaranteed to compile, be correct, or stable.
2001-08-02 Jeffrey D. Oldham <oldham at codesourcery.com>
* FieldCentering.cmpl.cpp
(CanonicalCentering::CanonicalCentering): Fix continuity label for
canonical discontinuous vertex centering.
* Field.h: Include NearestNeighbors.h
(View2<..., FieldOffset<...> >): New class merging previous
FieldOffset<...,[01]> code.
(View2<..., FieldOffset<...,1> >): Removed.
(View2<..., FieldOffset<...,0> >): Removed.
* FieldOffset.h: (FieldOffset<...>): New class merging previous
FieldOffset<...,[01]> code.
(FieldOffset<...,1>): Removed.
(FieldOffset<...,0>): Removed.
(FieldOffset::FieldOffset()): New function.
(FieldOffset::setSubFieldNumber): Likewise.
(FieldOffset::modifyCellOffset): Likewise.
(operator<<(...,FieldOffset)): Likewise.
(operator==(FieldOffset,...)): Likewise.
(operator!=(FieldOffset,...)): Likewise.
(FieldOffsetList): New class.
(operator<<(...,FieldOffsetList)): New function.
(accumulate): Likewise.
(sum): Likewise.
(average): Likewise.
(minimum): Likewise.
(maximum): Likewise.
* NearestNeighbors.h (NearestNeighborClass): New class.
(nearestNeighbors): New function.
* objfile.mk ($(UNIQUE)_OBJS): Add FieldOffset.cmpl.o.
* tests/FieldOffset.cpp: Removed second template parameter from
all uses.
* tests/FieldReductions.cpp: New testing file.
* tests/NearestNeighbors.cpp: Likewise.
* tests/makefile: Add support for FieldReductions and
NearestNeighbors. Remove support for ScalarAdvection and
ScalarAdvectionXB.
Tested on sequential Linux using gcc3.0 by building Pooma library and NewField executables
also tested using KCC by Stephen Smith
Applied to newfield_revision development branch
Approved by Stephen Smith
Thanks,
Jeffrey D. Oldham
oldham at codesourcery.com
-------------- next part --------------
Index: Field.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/Field.h,v
retrieving revision 1.15.2.3
diff -c -p -r1.15.2.3 Field.h
*** Field.h 2001/07/24 19:50:44 1.15.2.3
--- Field.h 2001/08/02 19:48:50
***************
*** 68,73 ****
--- 68,74 ----
#include "NewField/FieldCreateLeaf.h"
#include "NewField/FieldCentering.h"
#include "NewField/FieldOffset.h"
+ #include "NewField/NearestNeighbors.h"
#include "NewField/PrintField.h"
*************** struct View1<Field<GeometryTag, T, Engin
*** 577,584 ****
};
! // FIXME: Add specializations for FieldOffset<dimensions, 1> and
! // FIXME: FieldOffset<dimensions, 0>.
//-----------------------------------------------------------------------------
// AltView1 avoids an instantiation problem that arises when two
--- 578,584 ----
};
! // FIXME: Add specializations for FieldOffset<dimensions>.
//-----------------------------------------------------------------------------
// AltView1 avoids an instantiation problem that arises when two
*************** struct View2<Field<GeometryTag, T, Engin
*** 740,752 ****
//-----------------------------------------------------------------------------
! // View2<Field, FieldOffset<Dim,1>, Loc<Dim> > specialization for
// indexing a field with a FieldOffset and a Loc.
//-----------------------------------------------------------------------------
template<class GeometryTag, class T, class EngineTag, int Dim>
struct View2<Field<GeometryTag, T, EngineTag>,
! FieldOffset<Dim, 1>,
Loc<Dim> >
{
// Convenience typedef for the thing we're taking a view of.
--- 740,752 ----
//-----------------------------------------------------------------------------
! // View2<Field, FieldOffset<Dim>, Loc<Dim> > specialization for
// indexing a field with a FieldOffset and a Loc.
//-----------------------------------------------------------------------------
template<class GeometryTag, class T, class EngineTag, int Dim>
struct View2<Field<GeometryTag, T, EngineTag>,
! FieldOffset<Dim>,
Loc<Dim> >
{
// Convenience typedef for the thing we're taking a view of.
*************** struct View2<Field<GeometryTag, T, Engin
*** 766,849 ****
inline static
Type_t make(const Subject_t &f,
! const FieldOffset<dimensions, 1> &fo,
const Loc<dimensions> &loc)
{
CTAssert(dimensions == Dim);
- PAssert(f.numSubFields() > 0);
#if POOMA_BOUNDS_CHECK
! PInsist(contains(f.totalDomain(), loc + fo.cellOffset()),
! "Field view bounds error.");
! #endif
! return f[fo.subFieldNumber()].engine()(loc + fo.cellOffset());
! }
!
! inline static
! ReadType_t makeRead(const Subject_t &f,
! const FieldOffset<dimensions, 1> &fo,
! const Loc<dimensions> &loc)
! {
! PAssert(f.numSubFields() > 0);
!
! #if POOMA_BOUNDS_CHECK
! PInsist(contains(f.totalDomain(), loc + fo.cellOffset()),
! "Field view bounds error.");
#endif
! return f[fo.subFieldNumber()].engine().read(loc + fo.cellOffset());
! }
! };
!
!
! //-----------------------------------------------------------------------------
! // View2<Field, FieldOffset<Dim,0>, Loc<Dim> > specialization for
! // indexing a field with a FieldOffset and a Loc.
! //-----------------------------------------------------------------------------
!
! template<class GeometryTag, class T, class EngineTag, int Dim>
! struct View2<Field<GeometryTag, T, EngineTag>,
! FieldOffset<Dim, 0>,
! Loc<Dim> >
! {
! // Convenience typedef for the thing we're taking a view of.
!
! typedef Field<GeometryTag, T, EngineTag> Subject_t;
!
! // The field's dimension (i.e., the number of indices required to select a point).
!
! enum { dimensions = Subject_t::dimensions };
!
! // The return types.
!
! typedef typename Subject_t::Element_t ReadType_t;
! typedef typename Subject_t::ElementRef_t Type_t;
!
! // The functions that do the indexing.
!
! inline static
! Type_t make(const Subject_t &f,
! const FieldOffset<dimensions, 0> &fo,
! const Loc<dimensions> &loc)
! {
! CTAssert(dimensions == Dim);
!
#if POOMA_BOUNDS_CHECK
! PInsist(contains(f.totalDomain(), loc + fo.cellOffset()),
! "Field view bounds error.");
#endif
! return f.engine()(loc + fo.cellOffset());
}
inline static
ReadType_t makeRead(const Subject_t &f,
! const FieldOffset<dimensions, 0> &fo,
const Loc<dimensions> &loc)
{
#if POOMA_BOUNDS_CHECK
! PInsist(contains(f.totalDomain(), loc + fo.cellOffset()),
! "Field view bounds error.");
#endif
! return f.engine().read(loc + fo.cellOffset());
}
};
--- 766,813 ----
inline static
Type_t make(const Subject_t &f,
! const FieldOffset<dimensions> &fo,
const Loc<dimensions> &loc)
{
CTAssert(dimensions == Dim);
+ if (f.numSubFields() > 0) {
#if POOMA_BOUNDS_CHECK
! PInsist(contains(f[fo.subFieldNumber()].totalDomain(),
! loc + fo.cellOffset()),
! "Field view bounds error.");
#endif
! return f[fo.subFieldNumber()].engine()(loc + fo.cellOffset());
! }
! else {
#if POOMA_BOUNDS_CHECK
! PInsist(contains(f.totalDomain(), loc + fo.cellOffset()),
! "Field view bounds error.");
! return f.engine()(loc + fo.cellOffset());
#endif
! }
}
inline static
ReadType_t makeRead(const Subject_t &f,
! const FieldOffset<dimensions> &fo,
const Loc<dimensions> &loc)
{
+ if (f.numSubFields() > 0) {
#if POOMA_BOUNDS_CHECK
! PInsist(contains(f[fo.subFieldNumber()].totalDomain(),
! loc + fo.cellOffset()),
! "Field view bounds error.");
! #endif
! return f[fo.subFieldNumber()].engine().read(loc + fo.cellOffset());
! }
! else {
! #if POOMA_BOUNDS_CHECK
! PInsist(contains(f.totalDomain(), loc + fo.cellOffset()),
! "Field view bounds error.");
#endif
! return f.engine().read(loc + fo.cellOffset());
! }
}
};
Index: FieldCentering.cmpl.cpp
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/Attic/FieldCentering.cmpl.cpp,v
retrieving revision 1.1.2.4
diff -c -p -r1.1.2.4 FieldCentering.cmpl.cpp
*** FieldCentering.cmpl.cpp 2001/07/17 23:16:10 1.1.2.4
--- FieldCentering.cmpl.cpp 2001/08/02 19:48:50
*************** CanonicalCentering<Dim>::CanonicalCenter
*** 285,290 ****
--- 285,295 ----
centering.addValue(orientation, position);
centering_table_m[VertexType][Continuous][AllDim%(1<<Dim)] =
centering;
+
+ centering = Centering<Dim>(VertexType, Discontinuous);
+ orientation = 0;
+ position = 0.0;
+ centering.addValue(orientation, position);
position(0) = 1.0; centering.addValue(orientation, position);
if (Dim > 1) {
position(1) = 1.0; centering.addValue(orientation, position);
Index: FieldOffset.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/Attic/FieldOffset.h,v
retrieving revision 1.1.2.1
diff -c -p -r1.1.2.1 FieldOffset.h
*** FieldOffset.h 2001/07/24 19:50:44 1.1.2.1
--- FieldOffset.h 2001/08/02 19:48:51
***************
*** 29,34 ****
--- 29,41 ----
//-----------------------------------------------------------------------------
// Classes:
// FieldOffset
+ // FieldOffsetList
+ // Functions:
+ // accumulate
+ // sum
+ // average
+ // minimum
+ // maximum
//-----------------------------------------------------------------------------
#ifndef POOMA_NEWFIELD_OFFSET_H
***************
*** 39,44 ****
--- 46,55 ----
//
// FieldOffset
// - specifies a relative cell offset and subfield number
+ // FieldOffsetList
+ // - is a sequence of FieldOffset's
+ // FieldOffsetList reductions
+ // - computations using the entries in a FieldOffsetList
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
***************
*** 46,73 ****
//-----------------------------------------------------------------------------
#include "Domain/Loc.h"
//-----------------------------------------------------------------------------
// Forward declarations:
//-----------------------------------------------------------------------------
! template <int Dim, int Type=1>
class FieldOffset;
//-----------------------------------------------------------------------------
! // Full Description of FieldOffset:
//
// Given a field f, a Loc loc, and a field offset (offset,num), a
// field value can be obtained. Since each value specified by the
// field's centering is stored in a separate subfield, the notation
! // f[num](loc + offset) yields the value.
//
! // Accessing values for fields with exactly one value per cell differs
! // from accessing fields with multiple subfields. If a field has
! // exactly one value per cell, use FieldOffset<Dim, 0>, which does not
! // store a subfield number. If a field has multiple subfields, use
! // FieldOffset<Dim, 1>, which stores a subfield number.
//
//-----------------------------------------------------------------------------
--- 57,134 ----
//-----------------------------------------------------------------------------
#include "Domain/Loc.h"
+ #include <iostream>
+ #include <vector>
+ #include <functional>
+ #include <algorithm>
+ #include <iterator>
+
//-----------------------------------------------------------------------------
// Forward declarations:
//-----------------------------------------------------------------------------
! template <int Dim>
class FieldOffset;
+ template <int Dim>
+ class FieldOffsetList;
+
//-----------------------------------------------------------------------------
! // Full Description of FieldOffset<Dim>:
//
// Given a field f, a Loc loc, and a field offset (offset,num), a
// field value can be obtained. Since each value specified by the
// field's centering is stored in a separate subfield, the notation
! //
! //-----------------------------------------------------------------------------
!
!
! //-----------------------------------------------------------------------------
! // Full Description of FieldOffsetList<Dim>:
! //
! // A set of FieldOffset's can be stored in a FieldOffsetList. The
! // list has a fixed size. The following member operations are
! // supported:
! // size_type size() const - return the number of FieldOffset's.
! // const_reference operator[](size_type n) - return the nth
! // FieldOffset.
! //
! //-----------------------------------------------------------------------------
!
!
! //-----------------------------------------------------------------------------
! // Full Description of FieldOffsetList<Dim> Reductions:
! //
! // FIXME: What do we do for data-parallel statements?
! //
! // Many computations using `FieldOffsetList's perform similar
! // computations. The provided reduction functions support these
! // computations. All take a field, a FieldOffsetList, and a Loc. The
! // Loc specifies a cell within the given field. Together with a
! // FieldOffsetList and a field, a field value is specified. The
! // reduction function combines all the specified field values
! // corresponding to `FieldOffset's in the list. The list must not be
! // empty.
! //
! // accumulate: Sequentially applies the given binary function to all
! // field offsets in the list.
! // sum: Adds the values indicated by the field offset list.
! // average: Averages the values indicated by the field offset
! // list. Note the division is computed using the element
! // type, e.g., integral or floating point division.
! // minimum: Returns the smallest of the indicated values.
! // maximum: Returns the largest of the indicated values.
//
! //-----------------------------------------------------------------------------
!
! //-----------------------------------------------------------------------------
! // Full Description of findFieldOffsetList():
//
+ // Given an input centering and an output centering,
+ // findFieldOffsetList() returns a std::vector<FieldOffsetList<Dim> >.
+ // Each vector entry corresponds to one output value in the output centering.
+ //
//-----------------------------------------------------------------------------
*************** class FieldOffset;
*** 76,89 ****
//-----------------------------------------------------------------------------
template <int Dim>
! class FieldOffset<Dim, 1> {
public:
//---------------------------------------------------------------------------
// User-callable constructors. These ctors are meant to be called by users.
FieldOffset(const Loc<Dim> &loc, const int subFieldNumber = 0)
! : cell_offset_m(loc), subfield_number_m (subFieldNumber)
{
#if POOMA_BOUNDS_CHECK
PInsist(subfield_number_m >= 0, "Erroneous FieldOffset subfield number.");
--- 137,150 ----
//-----------------------------------------------------------------------------
template <int Dim>
! class FieldOffset {
public:
//---------------------------------------------------------------------------
// User-callable constructors. These ctors are meant to be called by users.
FieldOffset(const Loc<Dim> &loc, const int subFieldNumber = 0)
! : cell_offset_m(loc), subfield_number_m(subFieldNumber)
{
#if POOMA_BOUNDS_CHECK
PInsist(subfield_number_m >= 0, "Erroneous FieldOffset subfield number.");
*************** public:
*** 92,97 ****
--- 153,177 ----
}
//---------------------------------------------------------------------------
+ // Internal POOMA constructor. These operations are used internally
+ // by POOMA. They are not really meant to be called by users.
+
+ FieldOffset()
+ : cell_offset_m(Loc<Dim>()), subfield_number_m(0)
+ {}
+
+ inline void setSubFieldNumber(int subFieldNumber)
+ {
+ subfield_number_m = subFieldNumber;
+ return;
+ }
+
+ inline Loc<Dim> &modifyCellOffset()
+ {
+ return cell_offset_m;
+ }
+
+ //---------------------------------------------------------------------------
// Accessors.
inline const Loc<Dim> &cellOffset() const
*************** private:
*** 114,153 ****
};
template <int Dim>
! class FieldOffset<Dim, 0> {
public:
! //---------------------------------------------------------------------------
! // User-callable constructors. These ctors are meant to be called by users.
! FieldOffset(const Loc<Dim> &loc, const int subFieldNumber = 0)
! : cell_offset_m(loc), subfield_number_m (-1)
! { }
//---------------------------------------------------------------------------
! // Accessors.
! inline const Loc<Dim> &cellOffset() const
! {
! return cell_offset_m;
! }
private:
! // The cell offset.
! Loc<Dim> cell_offset_m;
! // The subfield number, if appropriate.
! int subfield_number_m;
};
#endif // POOMA_NEWFIELD_OFFSET_H
// ACL:rcsinfo
// ----------------------------------------------------------------------
! // $RCSfile: FieldOffset.h,v $ $Author: oldham $
! // $Revision: 1.1.2.1 $ $Date: 2001/07/24 19:50:44 $
// ----------------------------------------------------------------------
// ACL:rcsinfo
--- 194,438 ----
};
+ //-----------------------------------------------------------------------------
+ // Overload the << operator to print a FieldOffset to a stream.
+ //-----------------------------------------------------------------------------
+
+ template <int Dim>
+ std::ostream &operator<<(std::ostream &o,
+ const FieldOffset<Dim> &fieldOffset)
+ {
+ return o << "FieldOffset: (" << fieldOffset.cellOffset()
+ << ", " << fieldOffset.subFieldNumber() << ")";
+ }
+
+ //-----------------------------------------------------------------------------
+ // Define equality and inequality operators for FieldOffsets/
+ //-----------------------------------------------------------------------------
+
+ template <int Dim>
+ inline bool
+ operator==(const FieldOffset<Dim> &fieldOffset1,
+ const FieldOffset<Dim> &fieldOffset2)
+ {
+ return
+ fieldOffset1.cellOffset() == fieldOffset2.cellOffset() &&
+ fieldOffset1.subFieldNumber() == fieldOffset2.subFieldNumber();
+ }
+
+ template <int Dim>
+ inline bool
+ operator!=(const FieldOffset<Dim> &fieldOffset1,
+ const FieldOffset<Dim> &fieldOffset2)
+ {
+ return !(fieldOffset1 == fieldOffset2);
+ }
+
+
+ //-----------------------------------------------------------------------------
+ // FieldOffsetList.
+ //-----------------------------------------------------------------------------
+
template <int Dim>
! class FieldOffsetList
! {
public:
! // Exported typedefs.
!
! typedef size_t size_type;
! typedef FieldOffset<Dim>& reference;
! typedef const FieldOffset<Dim>& const_reference;
! // Return the number of FieldOffset's.
!
! size_type size() const
! {
! return v_m.size();
! }
!
! // Return a FieldOffset.
!
! const_reference operator[](const size_type n) const
! {
! #if POOMA_BOUNDS_CHECK
! PInsist(n < size(), "Erroneous FieldOffsetList index.");
! #endif
! return v_m[n];
! }
//---------------------------------------------------------------------------
! // Internal POOMA operators. These operations are used internally
! // by POOMA. They are not really meant to be called by users.
! // Create an empty list. This is used for arrays or std::vectors.
!
! FieldOffsetList() {}
!
! // Create a list that can hold the specified number of entries.
!
! FieldOffsetList(const size_type sz)
! {
! #if POOMA_BOUNDS_CHECK
! PInsist(sz > 0, "Erroneous FieldOffsetList size.");
! #endif
! v_m.reserve(sz);
! }
+ // Copy a vector's entries to this FieldOffsetList.
+
+ FieldOffsetList &operator=(const std::vector<FieldOffset<Dim> > &v)
+ {
+ v_m.resize(v.size());
+ std::copy(v.begin(), v.end(), v_m.begin());
+ return *this;
+ }
+
+ // Permit adding the specified entry.
+
+ reference operator[](const size_type n)
+ {
+ #if POOMA_BOUNDS_CHECK
+ PInsist(n < size(), "Erroneous FieldOffsetList index.");
+ #endif
+ return v_m[n];
+ }
+
private:
! std::vector<FieldOffset<Dim> > v_m;
! };
! //-----------------------------------------------------------------------------
! // Overload the << operator to print a FieldOffsetList to a stream.
! //-----------------------------------------------------------------------------
!
! template <int Dim>
! std::ostream &operator<<(std::ostream &o,
! const FieldOffsetList<Dim> &fieldOffsetList)
! {
! o << "FieldOffsetList:\n";
! for (int index = 0; index < fieldOffsetList.size(); ++index)
! o << fieldOffsetList[index] << std::endl;
! return o;
! }
!
!
! //-----------------------------------------------------------------------------
! // FieldOffsetList Reductions.
! //-----------------------------------------------------------------------------
!
! // Accumulate all the specified field locations using the supplied STL
! // binary function. Then, for each FieldOffset fo in the list, result
! // = binary_op(result, fv), where fv is the corresponding field value.
!
! // FIXME: Add data-parallel code.
!
! template<class GeometryTag, class T, class Expr, int Dim,
! class BinaryFunction>
! inline
! typename Field<GeometryTag, T, Expr>::T_t
! accumulate(BinaryFunction binary_op,
! const Field<GeometryTag, T, Expr>& field,
! const FieldOffsetList<Dim> &lst,
! const Loc<Dim> &loc)
! {
! typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
! typedef typename FieldOffsetList<Dim>::size_type size_type;
! CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));
!
! const size_type lstLength = lst.size();
! PInsist(lstLength > 0, "accumulate must be given a nonempty list.");
! T_t init = field(lst[0], loc);
! for (size_type i = 1; i < lstLength ; ++i)
! init = binary_op(init, field(lst[i], loc));
! return init;
! }
!
!
! // Sum all the values at the field locations.
!
! template<class GeometryTag, class T, class Expr, int Dim>
! inline
! typename Field<GeometryTag, T, Expr>::T_t
! sum(const Field<GeometryTag, T, Expr>& field,
! const FieldOffsetList<Dim> &lst,
! const Loc<Dim> &loc)
! {
! typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
! CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));
! return accumulate(std::plus<T_t>(), field, lst, loc);
! }
!
!
! // Average all the values at the field locations. Note the return
! // value has the same type as the field types so integer division may
! // be used.
!
! template<class GeometryTag, class T, class Expr, int Dim>
! inline
! typename Field<GeometryTag, T, Expr>::T_t
! average(const Field<GeometryTag, T, Expr>& field,
! const FieldOffsetList<Dim> &lst,
! const Loc<Dim> &loc)
! {
! typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
! CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));
! return sum(field, lst, loc) / lst.size();
! }
!
!
! // Return the minimum value of the field locations.
!
! template <class T>
! struct fomin : public std::binary_function<T, T, T>
! {
! T operator()(const T &op1, const T &op2) const {
! return std::min(op1, op2);
! }
! };
!
! template<class GeometryTag, class T, class Expr, int Dim>
! inline
! typename Field<GeometryTag, T, Expr>::T_t
! minimum(const Field<GeometryTag, T, Expr>& field,
! const FieldOffsetList<Dim> &lst,
! const Loc<Dim> &loc)
! {
! typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
! CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));
! return accumulate(fomin<T_t>(), field, lst, loc);
! }
!
!
! // Return the maximum value of the field locations.
!
! template <class T>
! struct fomax : public std::binary_function<T, T, T>
! {
! T operator()(const T &op1, const T &op2) const {
! return std::max(op1, op2);
! }
};
+ template<class GeometryTag, class T, class Expr, int Dim>
+ inline
+ typename Field<GeometryTag, T, Expr>::T_t
+ maximum(const Field<GeometryTag, T, Expr>& field,
+ const FieldOffsetList<Dim> &lst,
+ const Loc<Dim> &loc)
+ {
+ typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
+ CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));
+ return accumulate(fomax<T_t>(), field, lst, loc);
+ }
+
#endif // POOMA_NEWFIELD_OFFSET_H
// ACL:rcsinfo
// ----------------------------------------------------------------------
! // $RCSfile: FieldCentering.h,v $ $Author: oldham $
! // $Revision: 1.1.2.1 $ $Date: 2001/07/16 20:44:59 $
// ----------------------------------------------------------------------
// ACL:rcsinfo
Index: NearestNeighbors.h
===================================================================
RCS file: NearestNeighbors.h
diff -N NearestNeighbors.h
*** /dev/null Tue May 5 14:32:27 1998
--- NearestNeighbors.h Thu Aug 2 13:48:51 2001
***************
*** 0 ****
--- 1,401 ----
+ // -*- 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
+
+ //-----------------------------------------------------------------------------
+ // Function:
+ // nearestNeighbors
+ // Class:
+ // NearestNeighborClass
+ //-----------------------------------------------------------------------------
+
+ #ifndef POOMA_NEWFIELD_NEAREST_NEIGHBORS_H
+ #define POOMA_NEWFIELD_NEAREST_NEIGHBORS_H
+
+ //-----------------------------------------------------------------------------
+ // Overview:
+ //
+ // nearestNeighbors
+ // - yields FieldOffsetLists corresponding to the specified centerings
+ // NearestNeighborClass
+ // - the class performing the work
+ //-----------------------------------------------------------------------------
+
+ //-----------------------------------------------------------------------------
+ // Includes:
+ //-----------------------------------------------------------------------------
+
+ #include "Domain/Loc.h"
+ #include "NewField/FieldCentering.h"
+ #include <vector>
+ #include <numeric>
+ #include <utility>
+ #include <functional>
+ #include <algorithm>
+
+
+ //-----------------------------------------------------------------------------
+ // Forward declarations:
+ //-----------------------------------------------------------------------------
+
+ //-----------------------------------------------------------------------------
+ // Full Description of nearestNeighbors():
+ //
+ // Given input and output centerings, this function computes the
+ // "first shell" of nearest neighbors for each output value. That is,
+ // for each output value, it computes the FieldOffsetList containing
+ // the input values that are closest, with respect to the Manhattan
+ // norm (l_1 norm). The FieldOffsetLists are stored in a std::vector
+ // in the same order as the output values occur in the output
+ // centering.
+ //
+ // If only values within the corresponding input cell are desired,
+ // specify a third parameter.
+ //
+ //
+ // Full Description of NearestNeighborClass:
+ //
+ // This class implements the work of nearestNeighbors().
+ //
+ //-----------------------------------------------------------------------------
+
+ //----------------------------------------------------------------------
+ // NearestNeighborClass.
+ //----------------------------------------------------------------------
+
+ template <int Dim, bool IntraCellOnly = false>
+ class NearestNeighborClass {
+ public:
+
+ typedef FieldOffset<Dim> FieldOffset_t;
+ typedef FieldOffsetList<Dim> FieldOffsetList_t;
+ typedef std::vector<FieldOffset_t> FieldOffset_vt;
+ typedef std::vector<FieldOffsetList_t> Answer_t;
+ typedef Centering<Dim> Center;
+ typedef Center::Positions Positions;
+ typedef Center::Position Position;
+
+ // To compute the set of input values, we maintain a set of input
+ // values and their differences from the output value. In fact, we
+ // maintain the input value index and its difference.
+
+ typedef std::pair<int, Position> MinimumPair;
+ typedef std::vector<MinimumPair> MinimumSet;
+
+
+ NearestNeighborClass(const Center &inputCentering,
+ const Center &outputCentering)
+ : inputCentering_m(inputCentering),
+ outputCentering_m(outputCentering)
+ {
+ PInsist(inputCentering_m.positions().size() > 0,
+ "The input centering must be non-empty.");
+ }
+
+ inline
+ Answer_t operator()()
+ {
+ Answer_t answer;
+ answer.resize(outputCentering_m.size());
+ Positions inputPositions = inputCentering_m.positions();
+ Positions outputPositions = outputCentering_m.positions();
+ Position positionDifference;
+ double minimumDistance;
+
+ // Determine nearest neighbors for each output value.
+
+ for (Answer_t::size_type outputIndex = 0;
+ outputIndex < outputCentering_m.size();
+ ++outputIndex)
+ {
+ // Compute all input values in the first shell.
+
+ MinimumSet minimumSet; // all input values in first shell
+ const Position outputValue = outputPositions[outputIndex];
+ typename Positions::size_type inputIndex = 0;
+
+ // Use the first input value to start computing the minimum.
+
+ positionDifference = inputPositions[inputIndex] - outputValue;
+ minimumDistance =
+ (IntraCellOnly ?
+ manhattanDistance<Manhattan>(positionDifference) :
+ manhattanDistance<ManhattanGrid>(positionDifference));
+ minimumSet.push_back(std::make_pair(inputIndex, positionDifference));
+
+ // Compute the minimum over the rest of the input values.
+
+ for (++inputIndex;
+ inputIndex < inputPositions.size();
+ ++inputIndex) {
+ positionDifference = inputPositions[inputIndex] - outputValue;
+ const double distance =
+ (IntraCellOnly ?
+ manhattanDistance<Manhattan>(positionDifference) :
+ manhattanDistance<ManhattanGrid>(positionDifference));
+ if (distance < minimumDistance + epsilon) {
+ if (distance < minimumDistance) {
+ minimumSet.clear();
+ minimumDistance = distance;
+ }
+ minimumSet.push_back(std::make_pair(inputIndex,
+ positionDifference));
+ }
+ }
+
+
+ // Convert the minimum set to a set of FieldOffsets.
+ // minimumSet has all the minimum distance locations.
+
+ FieldOffset_vt answerHolder;
+ if (IntraCellOnly) {
+ for (MinimumSet::size_type minIndex = 0;
+ minIndex < minimumSet.size();
+ ++minIndex)
+ answerHolder.push_back(FieldOffset_t(Loc<Dim>(0),
+ minimumSet[minIndex].first));
+ }
+ else {
+ FieldOffset_vt partialAnswer;
+ for (MinimumSet::size_type minIndex = 0;
+ minIndex < minimumSet.size();
+ ++minIndex)
+ {
+ // Compute the cell offsets, appending to the set of answers.
+
+ partialAnswer = computeCellOffsets(minimumSet[minIndex].first,
+ minimumSet[minIndex].second);
+ answerHolder.insert(answerHolder.end(),
+ partialAnswer.begin(), partialAnswer.end());
+ }
+
+ // Remove all duplicates from the answer set.
+
+ std::sort(answerHolder.begin(), answerHolder.end(),
+ CompareFieldOffset());
+ answerHolder.erase(std::unique(answerHolder.begin(),
+ answerHolder.end(),
+ EqualFieldOffset()),
+ answerHolder.end());
+ }
+
+ // Store the answer.
+
+ answer[outputIndex] = answerHolder;
+ }
+
+ return answer;
+ }
+
+ private:
+
+ // Given a difference between two positions in logical coordinate
+ // space, return the Manhattan norm distance taking into account
+ // that input values are repeated in every grid cell.
+
+ struct ManhattanGrid : public std::binary_function<double, double, double>
+ {
+ double operator()(const double totalSoFar, double coordinate) const {
+ const double absCoordinate = std::abs(coordinate);
+ return totalSoFar + std::min(absCoordinate, 1-absCoordinate);
+ }
+ };
+
+ // Given a difference between two positions in logical coordinate
+ // space, return the Manhattan norm distance not taking into account
+ // that input values are repeated in every grid cell.
+
+ struct Manhattan : public std::binary_function<double, double, double>
+ {
+ double operator()(const double totalSoFar, double coordinate) const {
+ return totalSoFar + std::abs(coordinate);
+ }
+ };
+
+ template <typename Distance>
+ inline static
+ double manhattanDistance(const Position &difference)
+ {
+ double answer;
+ for (int coordinate = Dim-1; coordinate >= 0; --coordinate)
+ answer = Distance()(answer, difference(coordinate));
+ return answer;
+ }
+
+ // Given an input value in the first shell and its position
+ // difference from the given output value, return a vector of
+ // FieldOffsets of input values in the first shell, taking into
+ // account the repetition of input values throughout the grid. This
+ // is non-trivial because
+ // 1) input values are replicated and multiple values may
+ // be the same distance away
+ // 2) the closest input location may be in a different cell
+
+ inline static const FieldOffset_vt
+ computeCellOffsets(const int inputValueIndex, const Position &difference)
+ {
+ // Start with one empty tuple.
+
+ FieldOffset_vt answer(1);
+ int numTuples = 1;
+
+ // Store the cell offsets for input values.
+
+ int cellOffsetCoordinates[2]; // For our problem, there can
+ // be at most two cell offsets for
+ // each dimension.
+ int numOffsets; // number of cell offsets in
+ // cellOffsetCoordinates
+
+ for (int dimension = 0; dimension < Dim; ++dimension) {
+ numOffsets =
+ convertDifferenceToCellOffsets(difference(dimension),
+ cellOffsetCoordinates);
+ PInsist(numOffsets >= 1 && numOffsets <= 2,
+ "Incorrect number of cell offsets");
+
+ if (numOffsets == 2)
+ // Duplicate the tuples.
+ answer.insert(answer.end(), answer.begin(), answer.end());
+
+ for (int coc = 0; coc < numOffsets; ++coc)
+ for (int tuple = 0; tuple < numTuples; ++tuple)
+ answer[numTuples * coc + tuple].modifyCellOffset()[dimension] =
+ cellOffsetCoordinates[coc];
+
+ numTuples *= numOffsets;
+ }
+
+ // Set the subField numbers.
+
+ for (int i = numTuples-1; i >= 0; --i)
+ answer[i].setSubFieldNumber(inputValueIndex);
+
+ return answer;
+ }
+
+ // Given one coordinate of a difference between two coordinates,
+ // return the corresponding cell offset(s), either one or two.
+
+ inline static
+ int convertDifferenceToCellOffsets(const double difference,
+ int cellOffsetCoordinate[])
+ {
+ if (difference < -0.5 - epsilon) {
+ cellOffsetCoordinate[0] = +1;
+ return 1;
+ }
+ else if (std::abs(difference + 0.5) < epsilon) {
+ cellOffsetCoordinate[0] = +1;
+ cellOffsetCoordinate[1] = 0;
+ return 2;
+ }
+ else if (difference <= 0.5 - epsilon) {
+ cellOffsetCoordinate[0] = 0;
+ return 1;
+ }
+ else if (std::abs(difference - 0.5) < epsilon) {
+ cellOffsetCoordinate[0] = 0;
+ cellOffsetCoordinate[1] = -1;
+ return 2;
+ }
+ else if (difference < 1.0 + epsilon) {
+ cellOffsetCoordinate[0] = -1;
+ return 1;
+ }
+ else
+ PInsist(0, "Out of range difference");
+ }
+
+ // Specify a partial order for FieldOffsets to use when removing
+ // duplicates.
+
+ struct CompareFieldOffset :
+ public std::binary_function<FieldOffset_t, FieldOffset_t, bool> {
+ bool operator()(const FieldOffset_t &op1, const FieldOffset_t &op2) {
+ return (op1.cellOffset() < op2.cellOffset()) ||
+ (op1.cellOffset() == op2.cellOffset() &&
+ op1.subFieldNumber() < op2.subFieldNumber());
+ }
+ };
+
+ // Specify an equality operator for FieldOffsets to use when removing
+ // duplicates.
+
+ struct EqualFieldOffset :
+ public std::binary_function<FieldOffset_t, FieldOffset_t, bool> {
+ bool operator()(const FieldOffset_t &op1, const FieldOffset_t &op2) {
+ return op1.cellOffset() == op2.cellOffset() &&
+ op1.subFieldNumber() == op2.subFieldNumber();
+ }
+ };
+
+ const Center &inputCentering_m;
+ const Center &outputCentering_m;
+
+ // Use epsilon when comparing floating-point numbers, which cannot
+ // be represented precisely.
+
+ static const double epsilon;
+
+ };
+
+
+ template <int Dim, bool IntraCellOnly>
+ const double
+ NearestNeighborClass<Dim, IntraCellOnly>::epsilon = 1.0e-08;
+
+
+ //-----------------------------------------------------------------------------
+ // nearestNeighbors.
+ //-----------------------------------------------------------------------------
+
+ template <int Dim>
+ std::vector<FieldOffsetList<Dim> >
+ nearestNeighbors(const Centering<Dim> &inputCentering,
+ const Centering<Dim> &outputCentering)
+ {
+ return NearestNeighborClass<Dim>(inputCentering, outputCentering)();
+ }
+
+ template <int Dim>
+ std::vector<FieldOffsetList<Dim> >
+ nearestNeighbors(const Centering<Dim> &inputCentering,
+ const Centering<Dim> &outputCentering,
+ const bool)
+ {
+ return NearestNeighborClass<Dim, true>(inputCentering, outputCentering)();
+ }
+
+ #endif // POOMA_NEWFIELD_NEAREST_NEIGHBORS_H
+
+ // ACL:rcsinfo
+ // ----------------------------------------------------------------------
+ // $RCSfile: NearestNeighbors.h,v $ $Author: oldham $
+ // $Revision: 1.1.2.1 $ $Date: 2001/07/16 20:44:59 $
+ // ----------------------------------------------------------------------
+ // ACL:rcsinfo
Index: tests/FieldOffset.cpp
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/tests/Attic/FieldOffset.cpp,v
retrieving revision 1.1.2.1
diff -c -p -r1.1.2.1 FieldOffset.cpp
*** tests/FieldOffset.cpp 2001/07/24 19:50:45 1.1.2.1
--- tests/FieldOffset.cpp 2001/08/02 19:48:51
*************** int main(int argc, char *argv[])
*** 56,115 ****
// Test a field with subfields.
tester.check("f[0](0,0)",
! f(FieldOffset<Dim,1>(Loc<Dim>(0), 0), Loc<Dim>(0)),
-1.0, 1.0e-8);
tester.check("f[0](0,0)",
! f(FieldOffset<Dim,1>(Loc<Dim>(2,1), 0), Loc<Dim>(-2,-1)),
-1.0, 1.0e-8);
tester.check("f[0](2,1)",
! f(FieldOffset<Dim,1>(Loc<Dim>(2,1), 0), Loc<Dim>(0)),
-1.0, 1.0e-8);
tester.check("f[1](0,0)",
! f(FieldOffset<Dim,1>(Loc<Dim>(0), 1), Loc<Dim>(0)),
-2.0, 1.0e-8);
tester.check("f[1](1,2)",
! f(FieldOffset<Dim,1>(Loc<Dim>(1,2), 1), Loc<Dim>(0)),
-2.0, 1.0e-8);
! f(FieldOffset<Dim,1>(Loc<Dim>(3,2), 0), Loc<Dim>(-1,-1)) = 1.3;
! f(FieldOffset<Dim,1>(Loc<Dim>(3,2), 1), Loc<Dim>(-1,-1)) = 10.3;
tester.check("f[0](2,1)",
! f(FieldOffset<Dim,1>(Loc<Dim>(2,1), 0), Loc<Dim>(0)),
1.3, 1.0e-08);
tester.check("f[1](2,1)",
! f(FieldOffset<Dim,1>(Loc<Dim>(2,1), 1), Loc<Dim>(0)),
10.3, 1.0e-08);
tester.check("f[0].read(2,1)",
! f.read(FieldOffset<Dim,1>(Loc<Dim>(2,1), 0), Loc<Dim>(0)),
1.3, 1.0e-08);
tester.check("f[1].read(2,1)",
! f.read(FieldOffset<Dim,1>(Loc<Dim>(2,1), 1), Loc<Dim>(0)),
10.3, 1.0e-08);
// Test a field with no subfields.
Field_t h(canonicalCentering<Dim>(CellType, Continuous, AllDim),
layout, Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
! h(FieldOffset<Dim,0>(Loc<Dim>(0,0)), Loc<Dim>(0,0)) = 1.3;
! h(FieldOffset<Dim,0>(Loc<Dim>(0,0)), Loc<Dim>(0,1)) = 2.3;
! h(FieldOffset<Dim,0>(Loc<Dim>(0,0)), Loc<Dim>(1,0)) = 2.8;
! h(FieldOffset<Dim,0>(Loc<Dim>(1,0)), Loc<Dim>(0,1)) = 3.3;
tester.check("h(0,0)",
! h(FieldOffset<Dim,0>(Loc<Dim>(-1,-1)), Loc<Dim>(1,1)),
1.3, 1.0e-08);
tester.check("h(0,1)",
! h(FieldOffset<Dim,0>(Loc<Dim>(0,1)), Loc<Dim>(0,0)),
2.3, 1.0e-08);
tester.check("h(1,0)",
! h(FieldOffset<Dim,0>(Loc<Dim>(0,1)), Loc<Dim>(1,-1)),
2.8, 1.0e-08);
tester.check("h(1,1)",
! h(FieldOffset<Dim,0>(Loc<Dim>(0,0)), Loc<Dim>(1,1)),
3.3, 1.0e-08);
tester.check("h.read(1,0)",
! h.read(FieldOffset<Dim,0>(Loc<Dim>(0,1)), Loc<Dim>(1,-1)),
2.8, 1.0e-08);
tester.check("h.read(1,1)",
! h.read(FieldOffset<Dim,0>(Loc<Dim>(0,0)), Loc<Dim>(1,1)),
3.3, 1.0e-08);
int ret = tester.results("FieldOffset");
--- 56,115 ----
// Test a field with subfields.
tester.check("f[0](0,0)",
! f(FieldOffset<Dim>(Loc<Dim>(0), 0), Loc<Dim>(0)),
-1.0, 1.0e-8);
tester.check("f[0](0,0)",
! f(FieldOffset<Dim>(Loc<Dim>(2,1), 0), Loc<Dim>(-2,-1)),
-1.0, 1.0e-8);
tester.check("f[0](2,1)",
! f(FieldOffset<Dim>(Loc<Dim>(2,1), 0), Loc<Dim>(0)),
-1.0, 1.0e-8);
tester.check("f[1](0,0)",
! f(FieldOffset<Dim>(Loc<Dim>(0), 1), Loc<Dim>(0)),
-2.0, 1.0e-8);
tester.check("f[1](1,2)",
! f(FieldOffset<Dim>(Loc<Dim>(1,2), 1), Loc<Dim>(0)),
-2.0, 1.0e-8);
! f(FieldOffset<Dim>(Loc<Dim>(3,2), 0), Loc<Dim>(-1,-1)) = 1.3;
! f(FieldOffset<Dim>(Loc<Dim>(3,2), 1), Loc<Dim>(-1,-1)) = 10.3;
tester.check("f[0](2,1)",
! f(FieldOffset<Dim>(Loc<Dim>(2,1), 0), Loc<Dim>(0)),
1.3, 1.0e-08);
tester.check("f[1](2,1)",
! f(FieldOffset<Dim>(Loc<Dim>(2,1), 1), Loc<Dim>(0)),
10.3, 1.0e-08);
tester.check("f[0].read(2,1)",
! f.read(FieldOffset<Dim>(Loc<Dim>(2,1), 0), Loc<Dim>(0)),
1.3, 1.0e-08);
tester.check("f[1].read(2,1)",
! f.read(FieldOffset<Dim>(Loc<Dim>(2,1), 1), Loc<Dim>(0)),
10.3, 1.0e-08);
// Test a field with no subfields.
Field_t h(canonicalCentering<Dim>(CellType, Continuous, AllDim),
layout, Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
! h(FieldOffset<Dim>(Loc<Dim>(0,0)), Loc<Dim>(0,0)) = 1.3;
! h(FieldOffset<Dim>(Loc<Dim>(0,0)), Loc<Dim>(0,1)) = 2.3;
! h(FieldOffset<Dim>(Loc<Dim>(0,0)), Loc<Dim>(1,0)) = 2.8;
! h(FieldOffset<Dim>(Loc<Dim>(1,0)), Loc<Dim>(0,1)) = 3.3;
tester.check("h(0,0)",
! h(FieldOffset<Dim>(Loc<Dim>(-1,-1)), Loc<Dim>(1,1)),
1.3, 1.0e-08);
tester.check("h(0,1)",
! h(FieldOffset<Dim>(Loc<Dim>(0,1)), Loc<Dim>(0,0)),
2.3, 1.0e-08);
tester.check("h(1,0)",
! h(FieldOffset<Dim>(Loc<Dim>(0,1)), Loc<Dim>(1,-1)),
2.8, 1.0e-08);
tester.check("h(1,1)",
! h(FieldOffset<Dim>(Loc<Dim>(0,0)), Loc<Dim>(1,1)),
3.3, 1.0e-08);
tester.check("h.read(1,0)",
! h.read(FieldOffset<Dim>(Loc<Dim>(0,1)), Loc<Dim>(1,-1)),
2.8, 1.0e-08);
tester.check("h.read(1,1)",
! h.read(FieldOffset<Dim>(Loc<Dim>(0,0)), Loc<Dim>(1,1)),
3.3, 1.0e-08);
int ret = tester.results("FieldOffset");
Index: tests/FieldReductions.cpp
===================================================================
RCS file: FieldReductions.cpp
diff -N FieldReductions.cpp
*** /dev/null Tue May 5 14:32:27 1998
--- FieldReductions.cpp Thu Aug 2 13:48:51 2001
***************
*** 0 ****
--- 1,151 ----
+ // -*- 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
+ //-----------------------------------------------------------------------------
+ // Test field reductions.
+ //-----------------------------------------------------------------------------
+
+
+ #include "Pooma/NewFields.h"
+ #include "Utilities/Tester.h"
+ #include <iostream>
+
+
+ // Check the sum, average, minimum, and maximum functions for a
+ // specified position.
+
+ template <class Geometry, class T, class Engine, int Dim>
+ inline bool
+ checkFieldPosition(const Field<Geometry, T, Engine> &f,
+ const FieldOffsetList<Dim> &fol,
+ const Loc<Dim> &loc,
+ const T sumAnswer, const T avAnswer,
+ const T minAnswer, const T maxAnswer,
+ const double tolerance)
+ {
+ return
+ std::abs(sum(f, fol, loc) - sumAnswer) < tolerance &&
+ std::abs(average(f, fol, loc) - avAnswer) < tolerance &&
+ std::abs(minimum(f, fol, loc) - minAnswer) < tolerance &&
+ std::abs(maximum(f, fol, loc) - maxAnswer) < tolerance;
+ }
+
+
+ int main(int argc, char *argv[])
+ {
+ Pooma::initialize(argc, argv);
+ Pooma::Tester tester(argc, argv);
+
+ const double eps = 1.0e-08; // checking tolerance
+ const int Dim = 2;
+ std::vector<FieldOffsetList<Dim> > nn;
+ std::vector<FieldOffsetList<3> > nn3;
+
+ Centering<2> inputCenteringTwo, outputCenteringTwo;
+ Centering<3> inputCenteringThree, outputCenteringThree;
+
+ Interval<Dim> physicalVertexDomain(4, 4);
+ DomainLayout<Dim> layout(physicalVertexDomain, GuardLayers<Dim>(1));
+ typedef Field<UniformRectilinear<Dim>, double, Brick> Field_t;
+
+
+ // Test 2D Discontinuous Vertex -> Continuous Vertex.
+
+ inputCenteringTwo = canonicalCentering<2>(VertexType, Discontinuous, AllDim);
+ outputCenteringTwo = canonicalCentering<2>(VertexType, Continuous, AllDim);
+ nn = nearestNeighbors(inputCenteringTwo, outputCenteringTwo);
+ Field_t g(inputCenteringTwo, layout,
+ Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+
+ g.all() = 2.0;
+ g = -1.0;
+ g(FieldOffset<Dim>(Loc<Dim>(1,1), 0), Loc<Dim>(0,0)) = 17.0;
+ tester.check("discontinuous vertex->continuous vertex",
+ checkFieldPosition(g, nn[0], Loc<Dim>(1),
+ 14.0, 3.5, -1.0, 17.0, eps));
+
+
+ // Test 2D Continuous Cell -> Continuous Cell.
+
+ inputCenteringTwo = canonicalCentering<2>(CellType, Continuous, AllDim);
+ outputCenteringTwo = canonicalCentering<2>(CellType, Continuous, AllDim);
+ nn = nearestNeighbors(inputCenteringTwo, outputCenteringTwo);
+ Field_t f(inputCenteringTwo, layout,
+ Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+
+ f.all() = 2.0;
+ f = -1.0;
+ f(FieldOffset<Dim>(Loc<Dim>(1,1), 0), Loc<Dim>(0,0)) = 17.0;
+ tester.check("cell->cell",
+ checkFieldPosition(f, nn[0], Loc<Dim>(1,1),
+ 17.0, 17.0, 17.0, 17.0, eps));
+
+
+ // Test 2D Discontinuous Face -> Continuous Edge.
+
+ inputCenteringTwo = canonicalCentering<2>(FaceType, Discontinuous, AllDim);
+ outputCenteringTwo = canonicalCentering<2>(EdgeType, Continuous, AllDim);
+ nn = nearestNeighbors(inputCenteringTwo, outputCenteringTwo);
+ Field_t h(inputCenteringTwo, layout,
+ Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+
+ h.all() = 2.0;
+ h = -1.0;
+ h(FieldOffset<Dim>(Loc<Dim>(1), 0), Loc<Dim>(0)) = 17.0;
+ tester.check("discontinuous face->edge",
+ checkFieldPosition(h, nn[0], Loc<Dim>(1),
+ -2.0, -1.0, -1.0, -1.0, eps));
+
+ // Test 3D Discontinuous Vertex -> Continuous Cell.
+
+ inputCenteringThree = canonicalCentering<3>(VertexType, Discontinuous, AllDim);
+ outputCenteringThree = canonicalCentering<3>(CellType, Continuous, AllDim);
+ nn3 = nearestNeighbors(inputCenteringThree, outputCenteringThree);
+
+ Interval<3> physicalVertexDomain3(4, 4, 4);
+ DomainLayout<3> layout3(physicalVertexDomain3, GuardLayers<3>(1));
+ Field<UniformRectilinear<3>, double, Brick>
+ G(inputCenteringThree, layout3, Vector<3>(0.0), Vector<3>(1.0, 2.0, 0.0));
+
+ G.all() = 2.0;
+ G = -1.0;
+ G(FieldOffset<3>(Loc<3>(1), 0), Loc<3>(0)) = 17.0;
+ tester.check("discontinuous vertex->cell",
+ checkFieldPosition(G, nn3[0], Loc<3>(1),
+ -46.0, -46.0/64.0, -1.0, 17.0, eps));
+
+ int ret = tester.results("FieldReductions");
+ Pooma::finalize();
+ return ret;
+ }
+
+ // ACL:rcsinfo
+ // ----------------------------------------------------------------------
+ // $RCSfile: FieldOffset.cpp,v $ $Author: oldham $
+ // $Revision: 1.1.2.1 $ $Date: 2001/07/24 19:50:45 $
+ // ----------------------------------------------------------------------
+ // ACL:rcsinfo
Index: tests/NearestNeighbors.cpp
===================================================================
RCS file: NearestNeighbors.cpp
diff -N NearestNeighbors.cpp
*** /dev/null Tue May 5 14:32:27 1998
--- NearestNeighbors.cpp Thu Aug 2 13:48:52 2001
***************
*** 0 ****
--- 1,325 ----
+ // -*- 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
+ //-----------------------------------------------------------------------------
+ // Test computing the set of nearest neighbors.
+ //-----------------------------------------------------------------------------
+
+ #include "Pooma/NewFields.h"
+ #include "Utilities/Tester.h"
+ #include <vector>
+ #include <algorithm>
+
+
+ // Instead of checking all the nearest neighbor's FieldOffset values,
+ // we check "random" values.
+
+ // Check if a FieldOffset is in the FieldOffsetList.
+ template <int Dim>
+ inline bool
+ checkForFieldOffset(const FieldOffsetList<Dim> &lst,
+ const FieldOffset<Dim> &offset)
+ {
+ for (FieldOffsetList<Dim>::size_type index = 0;
+ index < lst.size();
+ ++index)
+ if (lst[index] == offset)
+ return true;
+ return false;
+ }
+
+
+ // Check for a particular FieldOffset within a vector of FieldOffsetLists.
+ // The arguments should be:
+ // a tester object,
+ // a C-string describing the test(s),
+ // a vector of FieldOffsetLists,
+ // the vector's correct size,
+ // the index of the particular FieldOffsetList to check,
+ // the FieldOffsetList correct size,
+ // a FieldOffset that should be (or not be) present in the list,
+ // whether the FieldOffset should be present.
+
+ template <int Dim>
+ inline bool
+ checkFieldOffset(Pooma::Tester &tester,
+ const char *testExplanation,
+ const std::vector<FieldOffsetList<Dim> > &nn,
+ const std::vector<FieldOffsetList<Dim> >::size_type nnSize,
+ const std::vector<FieldOffsetList<Dim> >::size_type listNum,
+ const FieldOffsetList<Dim>::size_type listSize,
+ const FieldOffset<Dim> &offset,
+ const bool offsetPresent = true)
+ {
+ PInsist(listNum >= 0 && listNum < nnSize,
+ "Incorrect std::vector<FieldOffsetList> index.");
+
+ return
+ tester.check(testExplanation, nn.size() == nnSize) &&
+ tester.check(testExplanation, nn[listNum].size() == listSize) &&
+ tester.check(testExplanation,
+ checkForFieldOffset(nn[listNum], offset) == offsetPresent);
+ }
+
+
+ // Check that the distances are the same for all input values.
+
+ // Compute the Manhattan distance of a difference between positions.
+
+ template <int Dim>
+ inline double
+ manhattanDistance(const Vector<Dim> &difference)
+ {
+ double answer;
+ for (int coordinate = Dim-1; coordinate >= 0; --coordinate)
+ answer += std::abs(difference(coordinate));
+ return answer;
+ }
+
+ // Compute the Manhattan distance between an input centering's value
+ // shifted by an FieldOffset and an output centering's value.
+
+ template <int Dim>
+ inline double
+ manhattanDistance(const Centering<Dim> &inputCentering,
+ const FieldOffset<Dim> &offset,
+ const Centering<Dim> &outputCentering,
+ const int outputIndex)
+ {
+ // Compute the actual input position.
+ Loc<Dim> cellOffset = offset.cellOffset();
+ Vector<Dim> input =
+ inputCentering.positions()[offset.subFieldNumber()];
+ for (int index = Dim-1; index >= 0; --index)
+ input(index) += cellOffset[index].first();
+
+ return manhattanDistance(outputCentering.positions()[outputIndex] - input);
+ }
+
+ // Check that the distance between the input and output values are the
+ // same for all the input values.
+
+ template <int Dim>
+ inline bool
+ sameDistances(const std::vector<FieldOffsetList<Dim> > &nn,
+ const Centering<Dim> &inputCentering,
+ const Centering<Dim> &outputCentering)
+ {
+ typedef std::vector<FieldOffsetList<Dim> >::size_type nn_size_type;
+ typedef FieldOffsetList<Dim>::size_type fol_size_type;
+ PInsist(nn.size() == outputCentering.size(),
+ "Nearest neighbors and output centering must have the same length.");
+
+ for (nn_size_type outputIndex = 0; outputIndex < nn.size(); ++outputIndex)
+ {
+ const double distance =
+ manhattanDistance(inputCentering, nn[outputIndex][0],
+ outputCentering, outputIndex);
+ for (fol_size_type inputIndex = 1;
+ inputIndex < nn[outputIndex].size();
+ ++inputIndex)
+ if (std::abs(distance -
+ manhattanDistance(inputCentering,
+ nn[outputIndex][inputIndex],
+ outputCentering, outputIndex))
+ > 1.0e-08)
+ return false;
+ }
+
+ return true;
+ }
+
+
+ int main(int argc, char *argv[])
+ {
+ Pooma::initialize(argc, argv);
+ Pooma::Tester tester(argc, argv);
+
+ Centering<2> inputCenteringTwo, outputCenteringTwo;
+ Centering<3> inputCenteringThree, outputCenteringThree;
+
+ // Test 2D Continuous Cell -> Continuous Cell.
+
+ inputCenteringTwo = canonicalCentering<2>(CellType, Continuous, AllDim);
+ outputCenteringTwo = canonicalCentering<2>(CellType, Continuous, AllDim);
+ checkFieldOffset(tester, "cell->cell intracell",
+ nearestNeighbors(inputCenteringTwo, outputCenteringTwo,
+ true),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 1,
+ FieldOffset<2>(Loc<2>(0)));
+
+ checkFieldOffset(tester, "cell->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 1,
+ FieldOffset<2>(Loc<2>(0)));
+
+ // Test 2D Continuous Vertex -> Continuous Cell.
+
+ inputCenteringTwo = canonicalCentering<2>(VertexType, Continuous, AllDim);
+ outputCenteringTwo = canonicalCentering<2>(CellType, Continuous, AllDim);
+ checkFieldOffset(tester, "vertex->cell intracell",
+ nearestNeighbors(inputCenteringTwo, outputCenteringTwo,
+ true),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 1,
+ FieldOffset<2>(Loc<2>(0)));
+ tester.check("vertex->cell intracell distances",
+ sameDistances(nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo,
+ true),
+ inputCenteringTwo,
+ outputCenteringTwo));
+
+ checkFieldOffset(tester, "vertex->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<2>(Loc<2>(0)));
+ checkFieldOffset(tester, "vertex->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<2>(Loc<2>(1,1)));
+ tester.check("vertex->cell intercell distances",
+ sameDistances(nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ inputCenteringTwo,
+ outputCenteringTwo));
+
+
+ // Test 2D Discontinuous Vertex -> Continuous Cell.
+
+ inputCenteringTwo = canonicalCentering<2>(VertexType, Discontinuous, AllDim);
+ outputCenteringTwo = canonicalCentering<2>(CellType, Continuous, AllDim);
+ checkFieldOffset(tester, "discontinuous vertex->cell intracell",
+ nearestNeighbors(inputCenteringTwo, outputCenteringTwo,
+ true),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<2>(Loc<2>(0), 0));
+ checkFieldOffset(tester, "discontinuous vertex->cell intracell",
+ nearestNeighbors(inputCenteringTwo, outputCenteringTwo,
+ true),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<2>(Loc<2>(0), 3));
+ tester.check("discontinuous vertex->cell intracell distances",
+ sameDistances(nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo,
+ true),
+ inputCenteringTwo,
+ outputCenteringTwo));
+
+ checkFieldOffset(tester, "discontinuous vertex->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 16,
+ FieldOffset<2>(Loc<2>(0), 0));
+ checkFieldOffset(tester, "discontinuous vertex->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 16,
+ FieldOffset<2>(Loc<2>(0), 3));
+ checkFieldOffset(tester, "discontinuous vertex->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 16,
+ FieldOffset<2>(Loc<2>(0), 3));
+ checkFieldOffset(tester, "discontinuous vertex->cell intercell",
+ nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ /* vector: */ 1, 0,
+ /* FieldOffsetList: */ 16,
+ FieldOffset<2>(Loc<2>(-1,0), 3), false);
+ tester.check("discontinuous vertex->cell intercell distances",
+ sameDistances(nearestNeighbors(inputCenteringTwo,
+ outputCenteringTwo),
+ inputCenteringTwo,
+ outputCenteringTwo));
+
+
+ // Test 3D Continuous Face -> Continuous Edge.
+
+ inputCenteringThree = canonicalCentering<3>(FaceType, Continuous, AllDim);
+ outputCenteringThree = canonicalCentering<3>(EdgeType, Continuous, AllDim);
+ checkFieldOffset(tester, "face->edge intracell",
+ nearestNeighbors(inputCenteringThree,
+ outputCenteringThree,
+ true),
+ /* vector: */ 3, 1,
+ /* FieldOffsetList: */ 2,
+ FieldOffset<3>(Loc<3>(0), 2));
+ tester.check("face->edge intracell distances",
+ sameDistances(nearestNeighbors(inputCenteringThree,
+ outputCenteringThree,
+ true),
+ inputCenteringThree,
+ outputCenteringThree));
+
+ checkFieldOffset(tester, "face->edge intercell",
+ nearestNeighbors(inputCenteringThree,
+ outputCenteringThree),
+ /* vector: */ 3, 1,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<3>(Loc<3>(-1,0,0), 2));
+ checkFieldOffset(tester, "face->edge intercell",
+ nearestNeighbors(inputCenteringThree,
+ outputCenteringThree),
+ /* vector: */ 3, 2,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<3>(Loc<3>(-1,0,0), 1));
+ checkFieldOffset(tester, "face->edge intercell",
+ nearestNeighbors(inputCenteringThree,
+ outputCenteringThree),
+ /* vector: */ 3, 2,
+ /* FieldOffsetList: */ 4,
+ FieldOffset<3>(Loc<3>(-1,-1,-1), 1), false);
+ tester.check("face->edge intercell distances",
+ sameDistances(nearestNeighbors(inputCenteringThree,
+ outputCenteringThree),
+ inputCenteringThree,
+ outputCenteringThree));
+
+ int ret = tester.results("NearestNeighbors");
+ Pooma::finalize();
+ return ret;
+ }
+
+ // ACL:rcsinfo
+ // ----------------------------------------------------------------------
+ // $RCSfile: FieldOffset.cpp,v $ $Author: oldham $
+ // $Revision: 1.1.2.1 $ $Date: 2001/07/24 19:50:45 $
+ // ----------------------------------------------------------------------
+ // ACL:rcsinfo
Index: tests/makefile
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/tests/makefile,v
retrieving revision 1.11.2.3
diff -c -p -r1.11.2.3 makefile
*** tests/makefile 2001/07/24 19:50:45 1.11.2.3
--- tests/makefile 2001/08/02 19:48:52
*************** field_tests:: $(ODIR)/BasicTest1 $(ODIR)
*** 61,69 ****
$(ODIR)/WhereTest $(ODIR)/VectorTest \
$(ODIR)/ScalarCode $(ODIR)/StencilTests \
$(ODIR)/ExpressionTest $(ODIR)/Centerings \
! $(ODIR)/FieldOffset
-
###########################
.PHONY: BasicTest1
--- 61,69 ----
$(ODIR)/WhereTest $(ODIR)/VectorTest \
$(ODIR)/ScalarCode $(ODIR)/StencilTests \
$(ODIR)/ExpressionTest $(ODIR)/Centerings \
! $(ODIR)/FieldOffset $(ODIR)/FieldReductions \
! $(ODIR)/NearestNeighbors
###########################
.PHONY: BasicTest1
*************** FieldOffset: $(ODIR)/FieldOffset
*** 160,184 ****
$(ODIR)/FieldOffset: $(ODIR)/FieldOffset.o
$(LinkToSuite)
! .PHONY: FieldTour3
! FieldTour3: $(ODIR)/FieldTour3
! $(ODIR)/FieldTour3: $(ODIR)/FieldTour3.o
$(LinkToSuite)
! .PHONY: ScalarAdvection
! ScalarAdvection: $(ODIR)/ScalarAdvection
! $(ODIR)/ScalarAdvection: $(ODIR)/ScalarAdvection.o
$(LinkToSuite)
! .PHONY: ScalarAdvectionXB
! ScalarAdvectionXB: $(ODIR)/ScalarAdvectionXB
! $(ODIR)/ScalarAdvectionXB: $(ODIR)/ScalarAdvectionXB.o
$(LinkToSuite)
--- 160,184 ----
$(ODIR)/FieldOffset: $(ODIR)/FieldOffset.o
$(LinkToSuite)
! .PHONY: FieldReductions
! FieldReductions: $(ODIR)/FieldReductions
! $(ODIR)/FieldReductions: $(ODIR)/FieldReductions.o
$(LinkToSuite)
! .PHONY: FieldTour3
! FieldTour3: $(ODIR)/FieldTour3
! $(ODIR)/FieldTour3: $(ODIR)/FieldTour3.o
$(LinkToSuite)
! .PHONY: NearestNeighbors
! NearestNeighbors: $(ODIR)/NearestNeighbors
! $(ODIR)/NearestNeighbors: $(ODIR)/NearestNeighbors.o
$(LinkToSuite)
More information about the pooma-dev
mailing list