[newfield_revision Patch] Support Field Replication

Jeffrey Oldham oldham at codesourcery.com
Thu Aug 16 18:26:02 UTC 2001


Specifying one field's values by copying values in another field is a
frequent operation supported by this patch.  To support it, we
implement the resulting field as an engine based on the input field
rather than performing the computation when replicate() is called.

2001-08-16  Jeffrey D. Oldham  <oldham at codesourcery.com>

        * FieldOffset.h: Remove unnecessary header files.
        (replicate): New function.
        * DiffOps/FieldShiftEngine.h (FieldShiftSimple::make): Revise
        parameters.  Move error checking from View2<Field, FieldOffset,
        Centering>.
        (FieldShiftSimple::make(Expression, vector<FieldOffsetList>,
        Centering>)): New function analogous to above make() but
        supporting a vector of FieldOffsetList's.
        (View2<Field,FieldOffset,Centering>): Fix initial comment.
        (View2<Field,FieldOffset,Centering>::make): Move error checking to
        FieldShiftSimple::make.
        (View2<Field,FieldOffset,Centering>::makeRead): Likewise.
        (View2<Field,vector<FieldOffset>,Centering>): New class supporting
        a vector of FieldOffsetList's.
        * tests/makefile (run_tests): Add Replicate.
        (field_tests): Likewise.
        Add PHONY Replicate rule.
	* tests/Replicate.cpp: New file to test replicate().

Applied to      newfield_revision branch
Approved by	Stephen Smith
Tested on       sequential Linux using gcc 3.0.1 by compiling Pooma library and NewField tests

Thanks,
Jeffrey D. Oldham
oldham at codesourcery.com
-------------- next part --------------
Index: FieldOffset.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/Attic/FieldOffset.h,v
retrieving revision 1.1.2.4
diff -c -p -r1.1.2.4 FieldOffset.h
*** FieldOffset.h	2001/08/15 18:09:45	1.1.2.4
--- FieldOffset.h	2001/08/16 17:18:00
***************
*** 65,73 ****
  #include "Domain/Loc.h"
  #include <iostream>
  #include <vector>
- #include <functional>
- #include <algorithm>
- #include <iterator>
  
  
  //-----------------------------------------------------------------------------
--- 65,70 ----
*************** class FieldOffsetList;
*** 129,134 ****
--- 126,145 ----
  //-----------------------------------------------------------------------------
  
  //-----------------------------------------------------------------------------
+ // Full Description of replicate:
+ //
+ // Copy field values to the specified locations.  The first field
+ // parameter specifies the field supplying the values to replicate.
+ // The second std::vector<FieldOffsetList> parameter specifies, for
+ // each value in the returned field, which input field value to use.
+ // The vector's length must match the number of values in each output
+ // field's cell.  For example, the output field's first value is
+ // copied from the location specified by the vector's first list.  The
+ // third parameter indicates the returned field's centering.
+ //
+ //-----------------------------------------------------------------------------
+ 
+ //-----------------------------------------------------------------------------
  // Full Description of findFieldOffsetList():
  //
  // Given an input centering and an output centering,
*************** max(const Field<GeometryTag, T, Expr>& f
*** 438,443 ****
--- 449,483 ----
    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);
+ }
+ 
+ 
+ //-----------------------------------------------------------------------------
+ // replicate.
+ //-----------------------------------------------------------------------------
+ 
+ template <class GeometryTag, class T, class Expr, int Dim>
+ inline
+ typename
+ View2<Field<GeometryTag, T, Expr>, std::vector<FieldOffset<Dim> >,
+       Centering<Dim> >::Type_t
+ replicate(const Field<GeometryTag, T, Expr>& field,
+           const std::vector<FieldOffsetList<Dim> > &vec,
+           const Centering<Dim> &centering)
+ {
+   CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));
+   typedef typename std::vector<FieldOffsetList<Dim> >::size_type vsize_type;
+   PInsist(vec.size() > 0, "Cannot replicate no values.");
+   PInsist(vec.size() == centering.size(),
+ 	  "Vector and output centering sizes must match.");
+ 
+   std::vector<FieldOffset<Dim> > vecFO(vec.size());
+   for (vsize_type i = 0; i < vec.size(); ++i) {
+     PInsist(vec[i].size() == 1, "Can replicate only one value.");
+     vecFO[i] = vec[i][0];
+   }
+ 
+   return field(vecFO, centering);
  }
  
  
Index: DiffOps/FieldShiftEngine.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/DiffOps/Attic/FieldShiftEngine.h,v
retrieving revision 1.1.2.2
diff -c -p -r1.1.2.2 FieldShiftEngine.h
*** DiffOps/FieldShiftEngine.h	2001/08/13 21:14:30	1.1.2.2
--- DiffOps/FieldShiftEngine.h	2001/08/16 18:07:48
*************** struct FieldShiftSimple
*** 372,378 ****
    typedef Engine<outputDim, OutputElement_t, OutputEngineTag_t> SEngine_t;
  
    static inline
!   Type_t make(const Expression &f, const Loc<outputDim> &offset,
                const Centering<outputDim> &centering)
    {
      // This should be h(centering, f.mesh(), f.layout())
--- 372,379 ----
    typedef Engine<outputDim, OutputElement_t, OutputEngineTag_t> SEngine_t;
  
    static inline
!   Type_t make(const Expression &f,
! 	      const FieldOffset<outputDim> &s1,
                const Centering<outputDim> &centering)
    {
      // This should be h(centering, f.mesh(), f.layout())
*************** struct FieldShiftSimple
*** 382,400 ****
  
      // Could change this to loop over centerings.
  
!     GuardLayers<outputDim> og(f.fieldEngine().guardLayers());
      for (int d = 0; d < outputDim; d++)
!     {
!       og.lower(d) += offset[d].first();
!       og.upper(d) -= offset[d].first();
!     }
  
      // need to set domain???
      h.fieldEngine().guardLayers() = og;
!     h.fieldEngine().engine() = SEngine_t(f.engine(), offset, f.domain());
  
      return h;
    }
  };
  
  //-----------------------------------------------------------------------------
--- 383,488 ----
  
      // Could change this to loop over centerings.
  
! #   if POOMA_BOUNDS_CHECK
!     if (f.numSubFields() > 0)
!       {
! 	PInsist((s1.subFieldNumber() < f.numSubFields()) &&
! 		(s1.subFieldNumber() >= 0),
! 		"subField bounds error");
! 	PInsist(contains(f[s1.subFieldNumber()].totalDomain(),
! 			 f[s1.subFieldNumber()].domain() + s1.cellOffset()),
! 		"Field operator()(FieldOffset) bounds error");
!       }
!     else
!       {
! 	PInsist(s1.subFieldNumber() == 0,
! 		"subField bounds error");
! 	PInsist(contains(f.totalDomain(), f.domain() + s1.cellOffset()),
! 		"Field operator()(FieldOffset) bounds error");
!       }
! #   endif   
! 
!     Expression fld = 
!       (f.numSubFields() > 0) ? f[s1.subFieldNumber()] : f;
!     const Loc<outputDim> &offset = s1.cellOffset();
! 
!     GuardLayers<outputDim> og(fld.fieldEngine().guardLayers());
      for (int d = 0; d < outputDim; d++)
!       {
! 	og.lower(d) += offset[d].first();
! 	og.upper(d) -= offset[d].first();
!       }
  
      // need to set domain???
      h.fieldEngine().guardLayers() = og;
!     h.fieldEngine().engine() = SEngine_t(fld.engine(), offset, fld.domain());
! 
!     return h;
!   }
! 
!   static inline
!   Type_t make(const Expression &f,
! 	      const std::vector<FieldOffset<outputDim> > &vs1,
!               const Centering<outputDim> &centering)
!   {
!     typedef std::vector<FieldOffset<outputDim> >::size_type size_type;
! 
!     // This should be h(centering, f.mesh(), f.layout())
!     // (Ideally centering would come out of offset.)
! 
!     Type_t h(f, centering);
! 
!     // Could change this to loop over centerings.
! 
!     PInsist(vs1.size() == centering.size(),
! 	    "The FieldOffset vector's length must match the centering's size.");
! 
!     // This code should simplify when unified access to fields with
!     // one or more subfields is possible.
  
+     for (size_type s1Index = 0; s1Index < vs1.size(); ++s1Index) {
+       const FieldOffset<outputDim> s1 = vs1[s1Index];
+       Type_t hField = (h.numSubFields() > 0) ? h[s1Index] : h;
+ 
+ #   if POOMA_BOUNDS_CHECK
+       if (f.numSubFields() > 0)
+ 	{
+ 	  PInsist((s1.subFieldNumber() < f.numSubFields()) &&
+ 		  (s1.subFieldNumber() >= 0),
+ 		  "subField bounds error");
+ 	  PInsist(contains(f[s1.subFieldNumber()].totalDomain(),
+ 			   f[s1.subFieldNumber()].domain() + s1.cellOffset()),
+ 		  "Field operator()(FieldOffset) bounds error");
+ 	}
+       else
+ 	{
+ 	  PInsist(s1.subFieldNumber() == 0,
+ 		  "subField bounds error");
+ 	  PInsist(contains(f.totalDomain(), f.domain() + s1.cellOffset()),
+ 		  "Field operator()(FieldOffset) bounds error");
+ 	}
+ #   endif   
+ 
+       Expression fld = 
+ 	(f.numSubFields() > 0) ? f[s1.subFieldNumber()] : f;
+       const Loc<outputDim> &offset = s1.cellOffset();
+ 
+       GuardLayers<outputDim> og(fld.fieldEngine().guardLayers());
+       for (int d = 0; d < outputDim; d++)
+ 	{
+ 	  og.lower(d) += offset[d].first();
+ 	  og.upper(d) -= offset[d].first();
+ 	}
+ 
+       // need to set domain???
+       hField.fieldEngine().guardLayers() = og;
+       hField.fieldEngine().engine() =
+ 	SEngine_t(fld.engine(), offset, fld.domain());
+     }
+ 
      return h;
    }
+ 
  };
  
  //-----------------------------------------------------------------------------
*************** struct LeafFunctor<Engine<Dim, T, FieldS
*** 597,608 ****
  };
  
  //-----------------------------------------------------------------------------
! // View1<Field, FieldOffset> specialization for indexing a field with a
! // FieldOffset.
  //-----------------------------------------------------------------------------
  
  template<class GeometryTag, class T, class EngineTag, int Dim>
! struct View2<Field<GeometryTag, T, EngineTag>, FieldOffset<Dim>, Centering<Dim> >
  {
    // Convenience typedef for the thing we're taking a view of.
    
--- 685,697 ----
  };
  
  //-----------------------------------------------------------------------------
! // View2<Field, FieldOffset, Centering> specialization for indexing a
! // field with a FieldOffset.
  //-----------------------------------------------------------------------------
  
  template<class GeometryTag, class T, class EngineTag, int Dim>
! struct View2<Field<GeometryTag, T, EngineTag>, FieldOffset<Dim>,
!              Centering<Dim> >
  {
    // Convenience typedef for the thing we're taking a view of.
    
*************** struct View2<Field<GeometryTag, T, Engin
*** 621,650 ****
                const Centering<Dim> &c)
    {
      CTAssert(Dim == Subject_t::dimensions);
! 
!     if (f.numSubFields() > 0)
!     {
! #   if POOMA_BOUNDS_CHECK
!       PInsist((s1.subFieldNumber() < f.numSubFields()) &&
!               (s1.subFieldNumber() >= 0),
!               "subField bounds error");
!       PInsist(contains(f[s1.subFieldNumber()].totalDomain(),
!                        f[s1.subFieldNumber()].domain() + s1.cellOffset()),
!               "Field operator()(FieldOffset) bounds error");
! #   endif
!       return FieldShiftSimple<Subject_t>::make(f[s1.subFieldNumber()],
!                                                s1.cellOffset(), c);
!     }
!     else
!     {
! #   if POOMA_BOUNDS_CHECK
!       PInsist(s1.subFieldNumber() == 0,
!               "subField bounds error");
!       PInsist(contains(f.totalDomain(), f.domain() + s1.cellOffset()),
!               "Field operator()(FieldOffset) bounds error");
! #   endif
!       return FieldShiftSimple<Subject_t>::make(f, s1.cellOffset(), c);
!     }
    }
  
    inline static
--- 710,716 ----
                const Centering<Dim> &c)
    {
      CTAssert(Dim == Subject_t::dimensions);
!     return FieldShiftSimple<Subject_t>::make(f, s1, c);
    }
  
    inline static
*************** struct View2<Field<GeometryTag, T, Engin
*** 652,681 ****
                        const Centering<Dim> &c)
    {
      CTAssert(Dim == Subject_t::dimensions);
  
!     if (f.numSubFields() > 0)
!     {
! #   if POOMA_BOUNDS_CHECK
!       PInsist((s1.subFieldNumber() < f.numSubFields()) &&
!               (s1.subFieldNumber() >= 0),
!               "subField bounds error");
!       PInsist(contains(f[s1.subFieldNumber()].totalDomain(),
!                        f[s1.subFieldNumber()].domain() + s1.cellOffset()),
!               "Field operator()(FieldOffset) bounds error");
! #   endif
!       return FieldShiftSimple<Subject_t>::make(f[s1.subFieldNumber()],
!                                                s1.cellOffset(), c);
!     }
!     else
!     {
! #   if POOMA_BOUNDS_CHECK
!       PInsist(s1.subFieldNumber() == 0,
!               "subField bounds error");
!       PInsist(contains(f.totalDomain(), f.domain() + s1.cellOffset()),
!               "Field operator()(FieldOffset) bounds error");
! #   endif
!       return FieldShiftSimple<Subject_t>::make(f, s1.cellOffset(), c);
!     }
    }
  };
  
--- 718,764 ----
                        const Centering<Dim> &c)
    {
      CTAssert(Dim == Subject_t::dimensions);
+     return FieldShiftSimple<Subject_t>::make(f, s1, c);
+   }
+ };
+ 
+ //-----------------------------------------------------------------------------
+ // View2<Field, vector<FieldOffset>, Centering> specialization for indexing a
+ // field with a vector<FieldOffset>.
+ //-----------------------------------------------------------------------------
+ 
+ template<class GeometryTag, class T, class EngineTag, int Dim>
+ struct View2<Field<GeometryTag, T, EngineTag>, std::vector<FieldOffset<Dim> >,
+              Centering<Dim> >
+ {
+   // Convenience typedef for the thing we're taking a view of.
+   
+   typedef Field<GeometryTag, T, EngineTag> Subject_t;
+   typedef typename Subject_t::Engine_t Engine_t;
+ 
+   // The return types.
+ 
+   typedef Field<GeometryTag, T, FieldShift<Engine_t> > ReadType_t;
+   typedef Field<GeometryTag, T, FieldShift<Engine_t> > Type_t;
+ 
+   // The functions that do the indexing.
+ 
+   inline static
+   Type_t make(const Subject_t &f,
+ 	      const std::vector<FieldOffset<Dim> > &s1,
+               const Centering<Dim> &c)
+   {
+     CTAssert(Dim == Subject_t::dimensions);
+     return FieldShiftSimple<Subject_t>::make(f, s1, c);
+   }
  
!   inline static
!   ReadType_t makeRead(const Subject_t &f,
! 		      const std::vector<FieldOffset<Dim> > &s1,
!                       const Centering<Dim> &c)
!   {
!     CTAssert(Dim == Subject_t::dimensions);
!     return FieldShiftSimple<Subject_t>::make(f, s1, c);
    }
  };
  
Index: tests/Replicate.cpp
===================================================================
RCS file: Replicate.cpp
diff -N Replicate.cpp
*** /dev/null	Tue May  5 14:32:27 1998
--- Replicate.cpp	Thu Aug 16 11:18:00 2001
***************
*** 0 ****
--- 1,98 ----
+ // -*- 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 replicating field values.
+ //-----------------------------------------------------------------------------
+ 
+ #include "Pooma/NewFields.h"
+ #include "Utilities/Tester.h"
+ #include <vector>
+ 
+ 
+ int main(int argc, char *argv[])
+ {
+   Pooma::initialize(argc, argv);
+   Pooma::Tester tester(argc, argv);
+ 
+   const double epsilon = 1.0e-08;
+   const int Dim = 2;
+   Centering<Dim> inputCenteringTwo, outputCenteringTwo;
+   Interval<Dim> physicalVertexDomain(4, 4);
+   DomainLayout<Dim> layout(physicalVertexDomain, GuardLayers<Dim>(1));
+   typedef Field<UniformRectilinear<Dim>, double, Brick> Field_t;
+ 
+   // Test 2D Continuous Cell -> Discontinuous Edge.
+ 
+   inputCenteringTwo = canonicalCentering<Dim>(CellType, Continuous, AllDim);
+   outputCenteringTwo = canonicalCentering<Dim>(EdgeType, Discontinuous, AllDim);
+   Field_t g(outputCenteringTwo, layout,
+ 	    Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+   Field_t f(inputCenteringTwo, layout,
+ 	    Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+   f.all() = 2.0;
+   g.all() = 1.0;
+ 
+   g = replicate(f, nearestNeighbors(inputCenteringTwo,
+ 				    outputCenteringTwo, true),
+ 		outputCenteringTwo);
+ 
+   tester.check("cell->discontinuous edge",
+ 	       g(FieldOffset<Dim>(Loc<Dim>(0), 0), Loc<Dim>(0)),
+ 	       2.0, epsilon);
+ 
+   // Test 2D Continuous Vertex -> Discontinuous Vertex.
+ 
+   inputCenteringTwo =
+     canonicalCentering<Dim>(VertexType, Continuous, AllDim);
+   outputCenteringTwo =
+     canonicalCentering<Dim>(VertexType, Discontinuous, AllDim);
+ 
+   Field_t g2(outputCenteringTwo, layout,
+ 	     Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+   Field_t f2(inputCenteringTwo, layout,
+ 	     Vector<Dim>(0.0), Vector<Dim>(1.0, 2.0));
+   f2.all() = 2.0;
+   g2.all() = 1.0;
+ 
+   g2 = replicate(f2, nearestNeighbors(inputCenteringTwo, outputCenteringTwo),
+ 		 outputCenteringTwo);
+   tester.check("vertex->discontinuous vertex",
+ 	       g2(FieldOffset<Dim>(Loc<Dim>(0), 0), Loc<Dim>(0)),
+ 	       2.0, epsilon);
+ 
+   int ret = tester.results("Replicate");
+   Pooma::finalize();
+   return ret; 
+ }
+ 
+ // ACL:rcsinfo
+ // ----------------------------------------------------------------------
+ // $RCSfile: Replicate.cpp,v $   $Author: oldham $
+ // $Revision: 1.1.2.2 $   $Date: 2001/08/14 20:24:18 $
+ // ----------------------------------------------------------------------
+ // ACL:rcsinfo
Index: tests/makefile
===================================================================
RCS file: /home/pooma/Repository/r2/src/NewField/tests/makefile,v
retrieving revision 1.11.2.8
diff -c -p -r1.11.2.8 makefile
*** tests/makefile	2001/08/14 20:24:18	1.11.2.8
--- tests/makefile	2001/08/16 17:18:00
*************** run_tests: tests
*** 57,62 ****
--- 57,63 ----
  	$(ODIR)/Gradient $(TSTOPTS) 1>Gradient.out 2>&1
  	$(ODIR)/MeshTest1 $(TSTOPTS) 1>MeshTest1.out 2>&1
  	$(ODIR)/NearestNeighbors $(TSTOPTS) 1>NearestNeighbors.out 2>&1
+ 	$(ODIR)/Replicate $(TSTOPTS) 1>Replicate.out 2>&1
  	$(ODIR)/ScalarCode $(TSTOPTS) 1>ScalarCode.out 2>&1
  	$(ODIR)/StencilTests $(TSTOPTS) 1>StencilTests.out 2>&1
  	$(ODIR)/VectorTest $(TSTOPTS) 1>VectorTest.out 2>&1
*************** field_tests:: $(ODIR)/BasicTest1 $(ODIR)
*** 68,73 ****
--- 69,75 ----
  	$(ODIR)/FieldTour1 $(ODIR)/FieldTour2 \
  	$(ODIR)/FieldTour3 $(ODIR)/Gradient\
  	$(ODIR)/MeshTest1 $(ODIR)/NearestNeighbors \
+ 	$(ODIR)/Replicate \
  	$(ODIR)/ScalarCode $(ODIR)/StencilTests \
  	$(ODIR)/VectorTest $(ODIR)/WhereTest
  
*************** $(ODIR)/NearestNeighbors: $(ODIR)/Neares
*** 169,174 ****
--- 171,183 ----
  Positions: $(ODIR)/Positions
  
  $(ODIR)/Positions: $(ODIR)/Positions.o
+ 	$(LinkToSuite)
+ 
+ .PHONY: Replicate
+ 
+ Replicate: $(ODIR)/Replicate
+ 
+ $(ODIR)/Replicate: $(ODIR)/Replicate.o
  	$(LinkToSuite)
  
  .PHONY: ScalarCode


More information about the pooma-dev mailing list