SpatialOps
Grid.cpp
1 /*
2  * Copyright (c) 2014-2017 The University of Utah
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to
6  * deal in the Software without restriction, including without limitation the
7  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
8  * sell copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
19  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
20  * IN THE SOFTWARE.
21  */
22 
23 #include "Grid.h"
25 
26 namespace SpatialOps{
27 
28  template< typename DirT, typename FieldT >
29  double shift(){
30  // obtain the coordinate shift for the given direction.
31  // note that (0,0,0) shifting corresponds to a (dx/2,dy/2,dz/2)
32  // offset whereas (-1,-1,-1) shifting corresponds to (0,0,0).
33  return 0.5*double(1+IndexStagger<FieldT,DirT>::value);
34  };
35 
36  Grid::Grid( const IntVec& npts,
37  const DoubleVec length,
38  const DoubleVec origin )
39  : npts_( npts ),
40  length_( length ),
41  spacing_( length_/npts_ ),
42  origin_( origin )
43  {}
44 
45  template<typename CoordT> unsigned int get_dir();
46  template<> unsigned int get_dir<XDIR>(){ return 0; }
47  template<> unsigned int get_dir<YDIR>(){ return 1; }
48  template<> unsigned int get_dir<ZDIR>(){ return 2; }
49 
50  //------------------------------------------------------------------
51 
52  template< typename CoordT >
53  double Grid::spacing() const
54  {
55  return spacing_[ get_dir<CoordT>() ];
56  }
57 
58  //------------------------------------------------------------------
59 
60  template< typename CoordT, typename FieldT >
61  void Grid::set_coord( FieldT& f ) const
62  {
63  const unsigned int dir = get_dir<CoordT>();
64  const double offset = shift<CoordT,FieldT>() * spacing_[dir] + origin_[dir];
65 
66  typedef typename FieldT::iterator FieldIter;
67 
68 # ifdef ENABLE_CUDA
69  // if the field is on GPU, move it to CPU, populate it, then sync it back.
70  // this is slow, but the Grid class isn't used much in production, and this
71  // could be done only during the setup phase rather than repeatedly.
72  const short devIx = f.active_device_index();
73  const bool isCPU = (devIx == CPU_INDEX);
74  if( !isCPU ) {
75  f.add_device( CPU_INDEX );
76  f.set_device_as_active( CPU_INDEX );
77  }
78 # endif
79 
80  FieldIter iter=f.begin();
81 
82  const MemoryWindow& mwInterior = f.window_without_ghost();
83  const MemoryWindow& mw = f.window_with_ghost();
84 
85  const IntVec lo(0,0,0);
86  const IntVec hi( mw.extent() );
87 
88  const int ixOffset = mwInterior.offset(dir);
89  IntVec ix;
90  for( ix[2]=lo[2]; ix[2]<hi[2]; ++ix[2] ){
91  for( ix[1]=lo[1]; ix[1]<hi[1]; ++ix[1] ){
92  for( ix[0]=lo[0]; ix[0]<hi[0]; ++ix[0] ){
93  *iter = spacing_[dir] * (ix[dir]-ixOffset) + offset;
94  ++iter;
95  }
96  }
97  }
98 # ifdef ENABLE_CUDA
99  if( !isCPU ){
100  f.validate_device( devIx );
101  f.set_device_as_active( devIx );
102  }
103 # endif
104  }
105 
106  //==================================================================
107  // Explicit template instantiation
108  //
109 # define DECLARE_COORD( FIELD ) \
110  template void Grid::set_coord< XDIR, FIELD >( FIELD& ) const; \
111  template void Grid::set_coord< YDIR, FIELD >( FIELD& ) const; \
112  template void Grid::set_coord< ZDIR, FIELD >( FIELD& ) const;
113 
114 # define DECLARE( VOL ) \
115  DECLARE_COORD( VOL ); \
116  DECLARE_COORD( FaceTypes<VOL>::XFace ); \
117  DECLARE_COORD( FaceTypes<VOL>::YFace ); \
118  DECLARE_COORD( FaceTypes<VOL>::ZFace );
119 
120  DECLARE( SVolField );
121  DECLARE( XVolField );
122  DECLARE( YVolField );
123  DECLARE( ZVolField );
124  DECLARE_COORD( VertexField );
125 
126  template double Grid::spacing<XDIR>() const;
127  template double Grid::spacing<YDIR>() const;
128  template double Grid::spacing<ZDIR>() const;
129  //
130  //==================================================================
131 
132 } // namespace SpatialOps
double spacing() const
obtain the grid spacing in the requested direction
Definition: Grid.cpp:53
Provides tools to index into a sub-block of memory.
Definition: MemoryWindow.h:63
void set_coord(FieldT &f) const
set the coordinates on the given field.
Definition: Grid.cpp:61
Grid(const IntVec &npts, const DoubleVec length, const DoubleVec origin=DoubleVec(0, 0, 0))
Definition: Grid.cpp:36
Abstracts a field.
Definition: SpatialField.h:132