SpatialOps
NeboStencilBuilder.h
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 #ifndef NEBO_STENCIL_BUILDER_H
24 #define NEBO_STENCIL_BUILDER_H
25 
26 //To define new stencils use the following four macros.
27 //NEBO_FIRST_POINT and NEBO_FIRST_IJK define a new Nebo stencil point collection.
28 //NEBO_ADD_POINT and NEBO_ADD_IJK add a stencil point to an existing Nebo stencil point collection.
29 //
30 //NEBO_FIRST_POINT and NEBO_ADD_POINT are designed to work together.
31 //Likewise, NEBO_FIRST_IJK and NEBO_ADD_IJK are designed to work together.
32 //They can be mixed, but it is not recommended.
33 //
34 //NEBO_FIRST_POINT and NEBO_ADD_POINT are designed to work with stencil points based off of field types.
35 //For examples, look at Stencil2Collection, Stencil4Collection, and FDStencilCollection definitions in this file.
36 //
37 //NEBO_FIRST_IJK and NEBO_ADD_IJK are designed to work on constant stencil points, whose shape does not change depending on field type.
38 //For examples, look at NullStencilCollection and the seven BoxFilter*StencilCollection definitions in this file.
39 
40 
41 
42 //Define a new stencil from an IndexTriplet
43 #define NEBO_FIRST_POINT_WO_TN(POINT) NeboStencilPointCollection< POINT, NeboNil >
44 #define NEBO_FIRST_POINT(POINT) typename NEBO_FIRST_POINT_WO_TN(POINT)
45 
46 //Add a point (IndexTriplet) to an existing stencil
47 #define NEBO_ADD_POINT(POINT) template AddPoint< POINT >::Result
48 
49 //Define a new stencil from three constant integers
50 #define NEBO_FIRST_IJK(X, Y, Z) NeboStencilPointCollection< IndexTriplet<X,Y,Z>, NeboNil >
51 
52 //Add a point (three constant integers) to an existing stencil
53 #define NEBO_ADD_IJK(X, Y, Z) AddPoint< IndexTriplet<X,Y,Z> >::Result
54 
55 namespace SpatialOps {
56 
61  template< typename List, typename IT >
62  struct ListSubtract
63  {
66  result;
67  };
68 
69  template< typename Point, typename IT >
71  {
73  result;
74  };
75 
85  template<typename OperatorType, typename PntCltnT, typename SrcFieldT, typename DestFieldT>
87  public:
88  typedef OperatorType type;
89  typedef PntCltnT PointCollectionType;
90  typedef SrcFieldT SrcFieldType;
91  typedef DestFieldT DestFieldType;
94 
95  // typedefs for when argument is a Nebo expression
96  template<typename Arg>
97  struct WithArg {
100  };
101 
102  // typedefs for when argument is a field
106 
111  : coefCollection_(coefs)
112  {}
113 
114  ~NeboStencilBuilder() {}
115 
119  CoefCollection const & coefs(void) const { return coefCollection_; };
120 
126  void apply_to_field( const SrcFieldType & src, DestFieldType & dest ) const {
127  dest <<= operator()(src);
128  }
129 
134  inline FieldResult operator ()( const SrcFieldType & src ) const {
135  return FieldResult(FieldStencil(FieldArg(src), coefs()));
136  }
137 
142  template<typename Arg>
143  inline typename WithArg<Arg>::Result
144  operator ()( const NeboExpression<Arg, SrcFieldType> & src ) const {
145  typedef typename WithArg<Arg>::Stencil Stencil;
146  typedef typename WithArg<Arg>::Result Result;
147  return Result(Stencil(src.expr(), coefs()));
148  }
149 
150  private:
151  const CoefCollection coefCollection_;
152  };
153 
154  template<typename OperatorType, typename SrcFieldType, typename DestFieldType>
156  // source field offset
157  typedef typename SrcFieldType::Location::Offset SrcOffset;
158  // destination field offset
159  typedef typename DestFieldType::Location::Offset DestOffset;
160  // low (first) stencil point location (relative to the destination point)
162  // high (second) stencil point location (relative to the destination point)
164  // collection of all stencil points in this stencil
165  typedef NEBO_FIRST_POINT(LowStPt)::NEBO_ADD_POINT(HighStPt) StPtCollection;
166  };
167 
168  template<typename OperatorType, typename SrcFieldType, typename DestFieldType>
170  template<bool Boolean, typename True, typename False>
171  struct TemplateIf;
172 
173  template<typename True, typename False>
174  struct TemplateIf<true, True, False> { True typedef result; };
175 
176  template<typename True, typename False>
177  struct TemplateIf<false, True, False> { False typedef result; };
178 
179  // source field offset
180  typedef typename SrcFieldType::Location::Offset SrcOffset;
181  // destination field offset
182  typedef typename DestFieldType::Location::Offset DestOffset;
183  // unit vectors
187  // first direction (unit vector)
188  typedef typename TemplateIf<((int)(SrcOffset::X) != (int)(DestOffset::X)),
189  XUnit,
190  YUnit>::result FirstDir;
191  // second direction (unit vector)
192  typedef typename TemplateIf<((int)(SrcOffset::Z) != (int)(DestOffset::Z)),
193  ZUnit,
194  YUnit>::result SecondDir;
195  // source offset in the first direction
197  // source offset in the second direction
199  // destination offset in the first direction
201  // destination offset in the second direction
203  // low value in the first direction
204  typedef typename GreaterThan<SrcInFirstDir,
205  DestInFirstDir>::result::Negate LoValInFirstDir;
206  // high value in the first direction
207  typedef typename LessThan<SrcInFirstDir,
208  DestInFirstDir>::result HiValInFirstDir;
209  // low value in the second direction
210  typedef typename GreaterThan<SrcInSecondDir,
211  DestInSecondDir>::result::Negate LoValInSecondDir;
212  // high value in the second direction
213  typedef typename LessThan<SrcInSecondDir,
214  DestInSecondDir>::result HiValInSecondDir;
215  // stencil point locations (relative to the destination point)
220  // collection of all stencil points in this stencil
221  typedef NEBO_FIRST_POINT(StPt1)::NEBO_ADD_POINT(StPt2)
222  ::NEBO_ADD_POINT(StPt3)::NEBO_ADD_POINT(StPt4) StPtCollection;
223  };
224 
225  template<typename OperatorType, typename SrcFieldType, typename DestFieldType>
227  typedef typename OperatorType::DirT DirT;
228  typedef typename UnitTriplet<DirT>::type DirVec;
229  typedef typename DirVec::Negate LowStPt;
230  typedef DirVec HighStPt;
231  typedef NEBO_FIRST_POINT(LowStPt)::NEBO_ADD_POINT(HighStPt) StPtCollection;
232  };
233 
243  template<typename OperatorType, typename PntCltnT, typename SrcFieldT, typename DestFieldT>
245  public:
246  typedef OperatorType type;
247  typedef PntCltnT PointCollectionType;
248  typedef SrcFieldT SrcFieldType;
249  typedef DestFieldT DestFieldType;
252 
253  // typedefs for when argument is a Nebo expression
254  template<typename Arg>
255  struct WithArg {
258  };
259 
260  // typedefs for when argument is a field
264 
269  : coefCollection_(coefs)
270  {}
271 
273 
277  CoefCollection const & coefs(void) const { return coefCollection_; };
278 
284  void apply_to_field( const SrcFieldType & src, DestFieldType & dest ) const {
285  dest <<= operator()(src);
286  }
287 
292  inline FieldResult operator ()( const SrcFieldType & src ) const {
293  return FieldResult(FieldStencil(FieldArg(src), coefs()));
294  }
295 
300  template<typename Arg>
301  inline typename WithArg<Arg>::Result
302  operator ()( const NeboExpression<Arg, SrcFieldType> & src ) const {
303  typedef typename WithArg<Arg>::Stencil Stencil;
304  typedef typename WithArg<Arg>::Result Result;
305  return Result(Stencil(src.expr(), coefs()));
306  }
307 
308  private:
309  const CoefCollection coefCollection_;
310  };
311 
320  template<typename PntCltnT, typename SrcFieldT, typename DestFieldT>
322  public:
323  typedef PntCltnT PointCollectionType;
324  typedef SrcFieldT SrcFieldType;
325  typedef DestFieldT DestFieldType;
326 
327  // typedefs for when argument is a Nebo expression
328  template<typename Arg>
329  struct WithArg {
332  };
333 
334  // typedefs for when argument is a field
338 
343  {}
344 
346 
352  void apply_to_field( const SrcFieldType & src, DestFieldType & dest ) const {
353  dest <<= operator()(src);
354  }
355 
360  inline FieldResult operator ()( const SrcFieldType & src ) const {
361  return FieldResult(FieldStencil(FieldArg(src)));
362  }
363 
368  template<typename Arg>
369  inline typename WithArg<Arg>::Result
370  operator ()( const NeboExpression<Arg, SrcFieldType> & src ) const {
371  typedef typename WithArg<Arg>::Stencil Stencil;
372  typedef typename WithArg<Arg>::Result Result;
373  return Result(Stencil(src.expr()));
374  }
375  };
376 
378  typedef NEBO_FIRST_IJK(0, 0, 0) StPtCollection;
379  };
380 
389  template<typename PntCltnT, typename SrcFieldT, typename DestFieldT>
391  public:
392  typedef PntCltnT PointCollectionType;
393  typedef SrcFieldT SrcFieldType;
394  typedef DestFieldT DestFieldType;
395 
396  // typedefs for when argument is a Nebo expression
397  template<typename Arg>
398  struct WithArg {
403  };
404 
405  // typedefs for when argument is a field
411 
416  {}
417 
419 
425  void apply_to_field( const SrcFieldType & src, DestFieldType & dest ) const {
426  dest <<= operator()(src);
427  }
428 
433  inline FieldResult operator ()( const SrcFieldType & src ) const {
434  return FieldResult(FieldAverage(FieldStencil(FieldArg(src)),
435  FieldScalar(1.0 / double(PointCollectionType::length))));
436  }
437 
442  template<typename Arg>
443  inline typename WithArg<Arg>::Result
444  operator ()( const NeboExpression<Arg, SrcFieldType> & src ) const {
445  typedef typename WithArg<Arg>::Stencil Stencil;
446  typedef typename WithArg<Arg>::Scalar Scalar;
447  typedef typename WithArg<Arg>::Average Average;
448  typedef typename WithArg<Arg>::Result Result;
449  return Result(Average(Stencil(src.expr()),
450  Scalar(1.0 / double(PointCollectionType::length))));
451  }
452  };
453 
455  typedef NEBO_FIRST_IJK(-1,-1,-1)::NEBO_ADD_IJK( 0,-1,-1)::NEBO_ADD_IJK( 1,-1,-1)
456  ::NEBO_ADD_IJK(-1, 0,-1)::NEBO_ADD_IJK( 0, 0,-1)::NEBO_ADD_IJK( 1, 0,-1)
457  ::NEBO_ADD_IJK(-1, 1,-1)::NEBO_ADD_IJK( 0, 1,-1)::NEBO_ADD_IJK( 1, 1,-1)
458  ::NEBO_ADD_IJK(-1,-1, 0)::NEBO_ADD_IJK( 0,-1, 0)::NEBO_ADD_IJK( 1,-1, 0)
459  ::NEBO_ADD_IJK(-1, 0, 0)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 1, 0, 0)
460  ::NEBO_ADD_IJK(-1, 1, 0)::NEBO_ADD_IJK( 0, 1, 0)::NEBO_ADD_IJK( 1, 1, 0)
461  ::NEBO_ADD_IJK(-1,-1, 1)::NEBO_ADD_IJK( 0,-1, 1)::NEBO_ADD_IJK( 1,-1, 1)
462  ::NEBO_ADD_IJK(-1, 0, 1)::NEBO_ADD_IJK( 0, 0, 1)::NEBO_ADD_IJK( 1, 0, 1)
463  ::NEBO_ADD_IJK(-1, 1, 1)::NEBO_ADD_IJK( 0, 1, 1)::NEBO_ADD_IJK( 1, 1, 1)
464  StPtCollection;
465  };
466 
468  typedef NEBO_FIRST_IJK(-1,-1, 0)::NEBO_ADD_IJK( 0,-1, 0)::NEBO_ADD_IJK( 1,-1, 0)
469  ::NEBO_ADD_IJK(-1, 0, 0)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 1, 0, 0)
470  ::NEBO_ADD_IJK(-1, 1, 0)::NEBO_ADD_IJK( 0, 1, 0)::NEBO_ADD_IJK( 1, 1, 0)
471  StPtCollection;
472  };
473 
475  typedef NEBO_FIRST_IJK(-1, 0,-1)::NEBO_ADD_IJK( 0, 0,-1)::NEBO_ADD_IJK( 1, 0,-1)
476  ::NEBO_ADD_IJK(-1, 0, 0)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 1, 0, 0)
477  ::NEBO_ADD_IJK(-1, 0, 1)::NEBO_ADD_IJK( 0, 0, 1)::NEBO_ADD_IJK( 1, 0, 1)
478  StPtCollection;
479  };
480 
482  typedef NEBO_FIRST_IJK( 0,-1,-1)::NEBO_ADD_IJK( 0, 0,-1)::NEBO_ADD_IJK( 0, 1,-1)
483  ::NEBO_ADD_IJK( 0,-1, 0)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 0, 1, 0)
484  ::NEBO_ADD_IJK( 0,-1, 1)::NEBO_ADD_IJK( 0, 0, 1)::NEBO_ADD_IJK( 0, 1, 1)
485  StPtCollection;
486  };
487 
489  typedef NEBO_FIRST_IJK(-1, 0, 0)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 1, 0, 0)
490  StPtCollection;
491  };
492 
494  typedef NEBO_FIRST_IJK( 0,-1, 0)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 0, 1, 0)
495  StPtCollection;
496  };
497 
499  typedef NEBO_FIRST_IJK( 0, 0,-1)::NEBO_ADD_IJK( 0, 0, 0)::NEBO_ADD_IJK( 0, 0, 1)
500  StPtCollection;
501  };
502 
503  template<typename OperatorType, typename SrcFieldType, typename DestFieldType>
505  // source field offset
506  typedef typename SrcFieldType::Location::Offset SrcOffset;
507  // destination field offset
508  typedef typename DestFieldType::Location::Offset DestOffset;
509  // minus-side stencil point location
511  // plus-side stencil point location
512  typedef typename GreaterThan<SrcOffset,
513  DestOffset>::result::Negate PlusPoint;
514  };
515 
516  //must be a finite difference (FD) stencil, pull direction and points from operator:
517  template<typename OperatorType, typename FieldType>
518  struct MaskShiftPoints<OperatorType, FieldType, FieldType> {
519  typedef typename OperatorType::type OpT;
520  // FDStencilCollection:
522  // minus-side stencil point location
523  typedef typename FDStencil::DirVec MinusPoint;
524  // plus-side stencil point location
525  typedef typename FDStencil::DirVec::Negate PlusPoint;
526  };
527 
536  template<typename OperatorT, typename SrcFieldT, typename DestFieldT>
538  public:
539  typedef OperatorT OperatorType;
540  typedef SrcFieldT SrcFieldType;
541  typedef DestFieldT DestFieldType;
542 
545 
548 
550 
553 
556 
561  {}
562 
564 
569  inline MinusResult minus( const SrcMask & src ) const {
570  return MinusResult(MinusShift(Mask(src)));
571  }
572 
577  inline PlusResult plus( const SrcMask & src ) const {
578  return PlusResult(PlusShift(Mask(src)));
579  }
580  };
581 
585  template<typename CurrentMode,typename Arg0, typename Arg1,typename FieldType>
586  struct NeboSwitch;
587  template<typename Arg0, typename Arg1, typename FieldType>
588  struct NeboSwitch<Initial, Arg0, Arg1, FieldType> {
589  FieldType typedef field_type;
591 #ifdef ENABLE_THREADS
593  typedef ResizeType;
594 #endif
595  NeboSwitch(Arg0 const& arg0, Arg1 const& arg1, int which): arg0_(arg0), arg1_(arg1), which_(which){}
596 #define INVOKE_ARG_METHOD(name) (which_ == 0?arg0_.name:arg1_.name)
597  inline GhostData ghosts_with_bc(void) const {
598  if(which_ < 2)
599  return INVOKE_ARG_METHOD(ghosts_with_bc());
600  else
601  return max(arg0_.ghosts_with_bc(), arg1_.ghosts_with_bc());
602  }
603  inline GhostData ghosts_without_bc(void) const {
604  if(which_ < 2)
605  return INVOKE_ARG_METHOD(ghosts_with_bc());
606  else
607  return max(arg0_.ghosts_without_bc(), arg1_.ghosts_without_bc());
608  }
609  inline bool has_extents(void) const { return INVOKE_ARG_METHOD(has_extents()); }
610  inline IntVec extents(void) const { return INVOKE_ARG_METHOD(extents()); }
611  inline IntVec has_bc(void) const { return INVOKE_ARG_METHOD(has_bc());}
612  inline SeqWalkType init(IntVec const & extents, GhostData const & ghosts, IntVec const & hasBC) const {
613  return SeqWalkType(arg0_.init(extents, ghosts, hasBC),
614  arg1_.init(extents, ghosts, hasBC),
615  which_,
616  get_range(arg0_.extents(), arg0_.ghosts_without_bc()),
617  get_range(arg1_.extents(), arg1_.ghosts_without_bc()));
618  }
619 #ifdef ENABLE_THREADS
620  inline ResizeType resize(void) const {
621  return ResizeType(arg0_.resize(),
622  arg1_.resize(),
623  which_,
624  get_range(arg0_.extents(), arg0_.ghosts_without_bc()),
625  get_range(arg1_.extents(), arg1_.ghosts_without_bc()));
626  }
627 #endif
628 #ifdef __CUDACC__
630  inline bool cpu_ready(void) const { return INVOKE_ARG_METHOD(cpu_ready()); }
631  inline bool gpu_ready(int const deviceIndex) const {return INVOKE_ARG_METHOD(gpu_ready(deviceIndex));}
632  inline GPUWalkType gpu_init(IntVec const & extents,GhostData const & ghosts,IntVec const & hasBC, int const deviceIndex, cudaStream_t const & lhsStream) const {
633  return GPUWalkType(arg0_.gpu_init(extents, ghosts, hasBC, deviceIndex, lhsStream),
634  arg1_.gpu_init(extents, ghosts, hasBC, deviceIndex, lhsStream),
635  which_,
636  get_range(arg0_.extents(), arg0_.ghosts_without_bc()),
637  get_range(arg1_.extents(), arg1_.ghosts_without_bc()));
638  }
639  inline void stream_wait_event(cudaEvent_t const & event) const {
640  arg0_.stream_wait_event(event); arg1_.stream_wait_event(event);
641  }
642 # ifdef NEBO_GPU_TEST
643  inline void gpu_prep(int const deviceIndex) const {
644  arg0_.gpu_prep(deviceIndex);
645  arg1_.gpu_prep(deviceIndex);
646  }
647 # endif
648 #endif
649  private:
650  Arg0 const arg0_;
651  Arg1 const arg1_;
652  int which_;
653  static inline GhostData get_range(IntVec const& extents, GhostData const& ghosts)
654  {
655  return GhostData(-ghosts.get_minus(0),
656  extents[0] + ghosts.get_plus(0),
657  -ghosts.get_minus(1),
658  extents[1] + ghosts.get_plus(1),
659  -ghosts.get_minus(2),
660  extents[2] + ghosts.get_plus(2));
661  }
662  };
663 #ifdef ENABLE_THREADS
664  template<typename Arg0, typename Arg1, typename FieldType>
665  struct NeboSwitch<Resize, Arg0, Arg1, FieldType> {
666  FieldType typedef field_type;
668  NeboSwitch(Arg0 const & arg0, Arg1 const & arg1, int which, GhostData const& range0, GhostData const& range1)
669  : arg0_(arg0), arg1_(arg1), range0_(range0), range1_(range1), which_(which)
670  {}
671  inline SeqWalkType init(IntVec const & extents,GhostData const & ghosts, IntVec const & hasBC) const {
672  return SeqWalkType(arg0_.init(extents, ghosts, hasBC),
673  arg1_.init(extents, ghosts, hasBC),
674  which_,
675  _overlap(range0_, get_range(extents, ghosts)),
676  _overlap(range1_, get_range(extents, ghosts)));
677  }
678  private:
679  Arg0 const arg0_;
680  Arg1 const arg1_;
681  GhostData const range0_;
682  GhostData const range1_;
683  int which_;
684  static inline GhostData _overlap(const GhostData& g1, const GhostData& g2)
685  {
686  return GhostData(std::max(g1.get_minus(0), g2.get_minus(0)),
687  std::min(g1.get_plus(0), g2.get_plus(0)),
688  std::max(g1.get_minus(1), g2.get_minus(1)),
689  std::min(g1.get_plus(1), g2.get_plus(1)),
690  std::max(g1.get_minus(2), g2.get_minus(2)),
691  std::min(g1.get_plus(2), g2.get_plus(2)));
692  }
693  static inline GhostData get_range(IntVec const& extents, GhostData const& ghosts)
694  {
695  return GhostData(-ghosts.get_minus(0),
696  extents[0] + ghosts.get_plus(0),
697  -ghosts.get_minus(1),
698  extents[1] + ghosts.get_plus(1),
699  -ghosts.get_minus(2),
700  extents[2] + ghosts.get_plus(2));
701  }
702  };
703 #endif
704  template<typename Arg0, typename Arg1, typename FieldType>
705  struct NeboSwitch<SeqWalk, Arg0, Arg1, FieldType> {
706  FieldType typedef field_type;
707  NeboSwitch(Arg0 const & arg0, Arg1 const & arg1, int which, GhostData const& range0, GhostData const& range1)
708  : arg0_(arg0), arg1_(arg1), range0_(range0), range1_(range1), which_(which)
709  {}
710  template<typename OptionalArgs>
711  inline bool eval(int const x, int const y, int const z) const {
712  if(which_ < 2)
713  {
714  return INVOKE_ARG_METHOD(template eval<OptionalArgs>(x,y,z));
715  }
716  else
717  {
718  return ((range0_.get_minus(0) <= x && x < range0_.get_plus(0) &&
719  range0_.get_minus(1) <= y && y < range0_.get_plus(1) &&
720  range0_.get_minus(2) <= z && z < range0_.get_plus(2))?arg0_.template eval<OptionalArgs>(x,y,z):0) ||
721  ((range1_.get_minus(0) <= x && x < range1_.get_plus(0) &&
722  range1_.get_minus(1) <= y && y < range1_.get_plus(1) &&
723  range1_.get_minus(2) <= z && z < range1_.get_plus(2))?arg1_.template eval<OptionalArgs>(x,y,z):0);
724  }
725  }
726  private:
727  Arg0 const arg0_;
728  Arg1 const arg1_;
729  GhostData const range0_;
730  GhostData const range1_;
731  int which_;
732  };
733 #ifdef __CUDACC__
734  template<typename Arg0, typename Arg1, typename FieldType>
735  struct NeboSwitch<GPUWalk, Arg0, Arg1, FieldType> {
736  FieldType typedef field_type;
737  NeboSwitch(Arg0 const & arg0, Arg1 const & arg1, int which, GhostData const& range0, GhostData const& range1)
738  : arg0_(arg0), arg1_(arg1), which_(which)
739  {
740  for(int i = 0; i < 3; i ++)
741  {
742  minus0[i] = range0.get_minus(i);
743  plus0[i] = range0.get_plus(i);
744  minus1[i] = range1.get_minus(i);
745  plus1[i] = range1.get_plus(i);
746  }
747  }
748  template<typename OptionalArgs>
749  __device__ inline bool eval(int const x, int const y, int const z) const {
750  if(which_ < 2)
751  {
752  return INVOKE_ARG_METHOD(template eval<OptionalArgs>(x,y,z));
753  }
754  else
755  {
756  return ((minus0[0] <= x && x < plus0[0] &&
757  minus0[1] <= y && y < plus0[1] &&
758  minus0[2] <= z && z < plus0[2])?arg0_.template eval<OptionalArgs>(x,y,z):0) ||
759  ((minus1[0] <= x && x < plus1[0] &&
760  minus1[1] <= y && y < plus1[1] &&
761  minus1[2] <= z && z < plus1[2])?arg1_.template eval<OptionalArgs>(x,y,z):0);
762  }
763  }
764  private:
765  Arg0 const arg0_;
766  Arg1 const arg1_;
767  int minus0[3];
768  int plus0[3];
769  int minus1[3];
770  int plus1[3];
771  int which_;
772  };
773 #endif
774 
780  template <typename SrcMaskType, typename DestMaskType, typename Dir = NODIR>
782  public:
783  typedef std::vector<IntVec> Points;
784  typedef SrcMaskType FromType;
785  typedef DestMaskType ToType;
786 
788  struct type{
789  typedef Dir DirT;
790  };
791  };
792 
794  typename FromType::field_type,
795  typename ToType::field_type> Shift;
796 
797  typedef typename Shift::MinusPoint MinusPointT;
798  typedef typename Shift::PlusPoint PlusPointT;
799 
800  typedef typename Shift::MinusShift MinusShiftT;
801  typedef typename Shift::PlusShift PlusShiftT;
802 
804 
806 
807  ConvertedMask(const SrcMaskType& src, BCSide side1, BCSide side2, BCSide side3):src_(src), face_(side1)
808  {
809  if(side1 == MINUS_SIDE) neg_face_ = PLUS_SIDE;
810  else neg_face_ = MINUS_SIDE;
811  if(side2 == NO_SIDE) side2 = PLUS_SIDE;
812  if(side3 == NO_SIDE) side3 = PLUS_SIDE;
813  all_ = ((side2 != -1 && side2 != side1) || (side3 != -1 && side3 != side1));
814  }
815  inline BCSide which_face()
816  {
817  return neg_face_;
818  }
819  inline const Points& points() const
820  {
821  return src_.points();
822  }
823  inline const SrcMaskType& src_mask() const{
824  return src_;
825  }
826  inline Result mask() const{
827  return Result(Selector(MinusShiftT(src_), PlusShiftT(src_), all_?2:face_!= MINUS_SIDE));
828  }
829  /* this function make a spatial mask from a converted mask */
830  inline DestMaskType to_mask(const typename DestMaskType::field_type& field) const{
831  std::vector<IntVec> points;
832  for(std::vector<IntVec>::const_iterator it = src_.points().begin();
833  it != src_.points().end();
834  it ++)
835  {
836  if(all_)
837  {
838  _add_point(*it - MinusPointT::int_vec(), points, field.window_without_ghost(), field.get_ghost_data());
839  _add_point(*it - PlusPointT::int_vec(), points, field.window_without_ghost(), field.get_ghost_data());
840  }
841  else if(face_ == MINUS_SIDE)
842  {
843  _add_point(*it - MinusPointT::int_vec(), points, field.window_without_ghost(), field.get_ghost_data());
844  }
845  else
846  {
847  _add_point(*it - PlusPointT::int_vec(), points, field.window_without_ghost(), field.get_ghost_data());
848  }
849  }
850  return DestMaskType(field, points);
851  }
852  inline bool all() const{ return all_ == 1; }
853  inline operator Result() const { return mask(); }
854  inline IntVec get_shift_vec(int idx) const
855  {
856  if(all_ && idx == 0) return MinusPointT::int_vec();
857  else if(all_ && idx == 1) return PlusPointT::int_vec();
858  else if(face_ == MINUS_SIDE) return MinusPointT::int_vec();
859  else return PlusPointT::int_vec();
860  }
861  template <typename FieldType, typename Point>
863  {
866  return ReturnType(ShiftType(mask().expr()));
867  }
868  private:
869  //Shift shift_;
870  const SrcMaskType& src_;
871  BCSide neg_face_, face_;
872  bool all_;
873  static inline void _add_point(const IntVec& point, std::vector<IntVec>& vec, const MemoryWindow& maskWindow_, const GhostData& ghosts_)
874  {
875 
876  if(!(point[0] >= -ghosts_.get_minus(0)) ||
877  !(point[0] < (signed int)(maskWindow_.extent(0)) + ghosts_.get_plus(0)) ||
878  !(point[1] >= -ghosts_.get_minus(1)) ||
879  !(point[1] < (signed int)(maskWindow_.extent(1)) + ghosts_.get_plus(1)) ||
880  !(point[2] >= -ghosts_.get_minus(2)) ||
881  !(point[2] < (signed int)(maskWindow_.extent(2)) + ghosts_.get_plus(2))) return;
882  vec.push_back(point);
883  }
884  };
885  template<typename OpT, typename SrcFieldT, typename DestFieldT>
886  struct GetOperatorDirImp{typedef NODIR result;};
887  template<typename OpT, typename SrcFieldT>
888  struct GetOperatorDirImp<OpT, SrcFieldT, SrcFieldT>{typedef typename OpT::type::DirT result;};
889  template<typename OpT>
892  };
893 
894  template <typename Dir, typename SrcFieldT, typename DestFieldT>
895  struct GetFDOperatorDir{typedef NODIR result;};
896  template <typename Dir, typename SrcFieldT>
897  struct GetFDOperatorDir<Dir, SrcFieldT, SrcFieldT>{typedef Dir result;};
898 
907  template<typename OperatorT>
909  public:
910  typedef OperatorT OperatorType;
911  typedef typename OperatorType::SrcFieldType PhiFieldType;
912  typedef typename OperatorType::DestFieldType GammaFieldType;
913 
914  typedef typename OperatorType::PointCollectionType PointCollection;
915 
916  typedef typename PointCollection::Last LowPoint;
917  typedef typename PointCollection::AllButLast NonLowPoints;
918  typedef typename PointCollection::First HighPoint;
919  typedef typename PointCollection::AllButFirst NonHighPoints;
920 
921  typedef typename Subtract<IndexTriplet<0,0,0>, LowPoint>:: result LowGammaPoint;
922  typedef typename Subtract<IndexTriplet<0,0,0>, HighPoint>::result HighGammaPoint;
923  typedef typename ListSubtract<NonLowPoints, LowPoint>:: result NonLowSrcPoints;
926 
928 
933 
934  typedef std::vector<IntVec> Points;
935  typedef Points::const_iterator PointIterator;
936 
938 
942  NeboBoundaryConditionBuilder(OperatorType const & op)
943  : lowCoef_(op.coefs().last()),
944  highCoef_(op.coefs().coef()),
945  minusGamma_(1.0),
946  plusGamma_(1.0),
947  minusPhi_(op.coefs().all_but_last()),
948  plusPhi_(op.coefs().others()),
949  shift_()
950  {}
951 
953 
962  template<typename ExprType, typename MaskType>
963  inline void apply_imp(MaskType mask,
964  PhiFieldType& phi,
965  const IntVec& shift,
967  BCSide side) const
968  {
969  if(side == MINUS_SIDE)
970  cpu_apply<ExprType, MinusGammaType, MinusPhiType>(mask.points(),
971  shift,
972  minusGamma_,
973  minusPhi_,
974  phi,
975  gamma,
976  lowCoef_);
977  else
978  cpu_apply<ExprType, PlusGammaType, PlusPhiType>(mask.points(),
979  shift,
980  plusGamma_,
981  plusPhi_,
982  phi,
983  gamma,
984  highCoef_);
985 
986  }
993  template<typename ExprType>
994  inline void operator()(ConvertedPhiMask mask,
995  PhiFieldType& phi,
996  const NeboExpression<ExprType, GammaFieldType>& gamma) const
997  {
998  if(mask.all())
999  throw std::runtime_error("can not apply boundary condition on more than one face");
1000  if(phi.active_device_index() == CPU_INDEX) {
1001  apply_imp(mask,
1002  phi,
1003  IntVec(0,0,0),
1004  gamma,
1005  mask.which_face());
1006  } else {
1007  if(mask.which_face() == MINUS_SIDE)
1008  phi <<= cond(mask.src_mask(), (minusGamma_(gamma) - minusPhi_(phi)) / lowCoef_)
1009  (phi);
1010  else
1011  phi <<= cond(mask.src_mask(), (plusGamma_(gamma) - plusPhi_(phi)) / highCoef_)
1012  (phi);
1013  }
1014  }
1015  template<typename MaskType, typename GammaType>
1016  inline void apply_indirect_conversion(MaskType& mask, PhiFieldType& phi, GammaType gamma, bool minus, int idx) const
1017  {
1018  IntVec conv_shift = mask.get_shift_vec(idx);
1019  if(minus)
1020  apply_imp(mask,
1021  phi,
1022  Shift::MinusPoint::int_vec() + conv_shift,
1023  gamma,
1024  MINUS_SIDE);
1025  else
1026  apply_imp(mask,
1027  phi,
1028  Shift::PlusPoint::int_vec() + conv_shift,
1029  gamma,
1030  PLUS_SIDE);
1031  }
1032  template<typename ExprType, typename SrcMaskType, typename Dir>
1033  inline void operator()(ConvertedMask<const SpatialMask<SrcMaskType>, const SpatialMask<GammaFieldType>, Dir> mask,
1034  PhiFieldType& phi,
1036  bool minus) const
1037  {
1038  if(phi.active_device_index() == CPU_INDEX)
1039  {
1040  apply_indirect_conversion(mask, phi, gamma, minus, 0);
1041  if(mask.all()) apply_indirect_conversion(mask, phi, gamma, minus, 1);
1042  }
1043  else
1044  {
1045  if(minus)
1046  phi <<= cond(mask.template shift<PhiFieldType, typename Shift::MinusPoint>(), (minusGamma_(gamma) - minusPhi_(phi)) / lowCoef_)
1047  (phi);
1048  else
1049  phi <<= cond(mask.template shift<PhiFieldType, typename Shift::PlusPoint>(), (plusGamma_(gamma) - plusPhi_(phi)) / highCoef_)
1050  (phi);
1051  }
1052 
1053  }
1061  template<typename ExprType>
1062  inline void operator()(SpatialMask<GammaFieldType> mask,
1063  PhiFieldType & phi,
1065  bool minus) const {
1066  if(phi.active_device_index() == CPU_INDEX) {
1067  if(minus)
1068  apply_imp(mask,
1069  phi,
1070  Shift::MinusPoint::int_vec(),
1071  gamma,
1072  MINUS_SIDE);
1073  else
1074  apply_imp(mask,
1075  phi,
1076  Shift::PlusPoint::int_vec(),
1077  gamma,
1078  PLUS_SIDE);
1079  }
1080  else {
1081  if(minus)
1082  phi <<= cond(shift_.minus(mask), (minusGamma_(gamma) - minusPhi_(phi)) / lowCoef_)
1083  (phi);
1084  else
1085  phi <<= cond(shift_.plus(mask), (plusGamma_(gamma) - plusPhi_(phi)) / highCoef_)
1086  (phi);
1087  }
1088  }
1089 
1097  inline void operator()(SpatialMask<GammaFieldType> mask,
1098  PhiFieldType & phi,
1099  const GammaFieldType & gamma,
1100  bool minus) const
1101  {
1102  typedef NeboConstField<Initial, GammaFieldType> GammaField;
1104  (*this)(mask, phi, GammaExpr(GammaField(gamma)), minus);
1105  }
1106 
1113  inline void operator()(ConvertedPhiMask mask,
1114  PhiFieldType & phi,
1115  const GammaFieldType & gamma) const
1116  {
1117  typedef NeboConstField<Initial, GammaFieldType> GammaField;
1119  (*this)(mask, phi, GammaExpr(GammaField(gamma)));
1120  }
1121 
1122  template<typename SrcMaskType, typename Dir>
1123  inline void operator()(ConvertedMask<const SpatialMask<SrcMaskType>, const SpatialMask<GammaFieldType>, Dir> mask,
1124  PhiFieldType& phi,
1125  const GammaFieldType& gamma,
1126  bool minus) const
1127  {
1128  typedef NeboConstField<Initial, GammaFieldType> GammaField;
1130  (*this)(mask, phi, GammaExpr(GammaField(gamma)), minus);
1131  }
1139  inline void operator()(const SpatialMask<GammaFieldType> mask,
1140  PhiFieldType & phi,
1141  const double gamma,
1142  bool minus) const
1143  {
1144  typedef NeboScalar<Initial, double> GammaScalar;
1146  (*this)(mask, phi, GammaExpr(GammaScalar(gamma)), minus);
1147  }
1148 
1156  inline void operator()(ConvertedPhiMask mask,
1157  PhiFieldType & phi,
1158  const double gamma) const
1159  {
1160  typedef NeboScalar<Initial, double> GammaScalar;
1162  (*this)(mask, phi, GammaExpr(GammaScalar(gamma)));
1163  }
1164 
1165  template<typename SrcMaskType, typename Dir>
1166  inline void operator()(ConvertedMask<const SpatialMask<SrcMaskType>, const SpatialMask<GammaFieldType>, Dir> mask,
1167  PhiFieldType& phi,
1168  const double gamma,
1169  bool minus) const
1170  {
1171  typedef NeboScalar<Initial, double> GammaScalar;
1173  (*this)(mask, phi, GammaExpr(GammaScalar(gamma)), minus);
1174  }
1175 
1176  private:
1177 
1188  template<typename ExprType, typename ShiftGamma, typename ShiftPhi>
1189  inline void cpu_apply( const Points & points,
1190  const IntVec & shift,
1191  const ShiftGamma & shiftGamma,
1192  const ShiftPhi & shiftPhi,
1193  PhiFieldType & phi,
1195  const double coef ) const
1196  {
1197  typedef NeboField<Initial, PhiFieldType> LhsTypeInit;
1198  typedef typename LhsTypeInit::SeqWalkType LhsType;
1199  LhsTypeInit lhsInit(phi);
1200  LhsType lhs = lhsInit.init();
1201 
1202  typedef typename ShiftGamma::template WithArg<ExprType>::Stencil GammaType;
1203  typedef typename ShiftPhi::FieldStencil PhiType;
1204  typedef DiffOp<Initial, GammaType, PhiType> DiffType;
1205  typedef NeboScalar<Initial, double> NumType;
1206  typedef DivOp<Initial, DiffType, NumType> RhsTypeInit;
1207  typedef typename RhsTypeInit::SeqWalkType RhsType;
1208  //arguments to init() call are meaningless, but they are not used at all...
1209  RhsType rhs = RhsTypeInit(DiffType(shiftGamma(gamma).expr(),
1210  shiftPhi(phi).expr()),
1211  NumType(coef)).init(IntVec(0,0,0),
1212  GhostData(0),
1213  IntVec(0,0,0));
1214 
1215  PointIterator ip = points.begin();
1216  PointIterator const ep = points.end();
1217  for(; ip != ep; ip++) {
1218  const int x = (*ip)[0] - shift[0];
1219  const int y = (*ip)[1] - shift[1];
1220  const int z = (*ip)[2] - shift[2];
1221  lhs.ref(x, y, z) = rhs.template eval</*TODO: Wrap Into Compile Time Optional Args*/void>(x, y, z);
1222  }
1223  }
1224 
1225  const double lowCoef_;
1226  const double highCoef_;
1227  const MinusGammaType minusGamma_;
1228  const PlusGammaType plusGamma_;
1229  const MinusPhiType minusPhi_;
1230  const PlusPhiType plusPhi_;
1231  const Shift shift_;
1232  };
1233  template <typename DestFieldLocation, typename SrcFieldLocation>
1235  template <int x, int y>
1236  struct Equals{
1237  enum {result = (x == y)?1:0};
1238  };
1239  typedef typename DestFieldLocation::Offset Dest_Offset;
1240  typedef typename SrcFieldLocation::Offset Src_Offset;
1244  };
1246  {
1247  if(Shift && face == NO_SIDE)
1248  throw std::runtime_error("trying to make invalid mask conversion");
1249  }
1250  };
1251  template <typename DstFieldType, typename SrcFieldType> struct MaskConversionFieldTypeCheck;
1252 #define ALLOW_MASK_CONVERSION(from, to) template <> struct MaskConversionFieldTypeCheck<to, from>{};
1253  ALLOW_MASK_CONVERSION(SVolField, XVolField)
1254  ALLOW_MASK_CONVERSION(SVolField, YVolField)
1255  ALLOW_MASK_CONVERSION(SVolField, ZVolField)
1256 
1257 #define ALLOW_MASK_CONVERSION_GROUP(dir)\
1258  ALLOW_MASK_CONVERSION(dir##VolField, dir##SurfXField)\
1259  ALLOW_MASK_CONVERSION(dir##VolField, dir##SurfYField)\
1260  ALLOW_MASK_CONVERSION(dir##VolField, dir##SurfZField)\
1261 
1262  ALLOW_MASK_CONVERSION_GROUP(X)
1263  ALLOW_MASK_CONVERSION_GROUP(Y)
1264  ALLOW_MASK_CONVERSION_GROUP(Z)
1265  ALLOW_MASK_CONVERSION_GROUP(S)
1266 
1274  template <typename DstFieldType, typename Dir, typename SrcMskType>
1275  static inline ConvertedMask<const SrcMskType,
1276  const SpatialMask<DstFieldType>,
1277  typename GetFDOperatorDir<Dir,
1278  DstFieldType,
1279  typename SrcMskType::field_type>::result
1280  > convert(const SrcMskType& src, const BCSide face1, const BCSide face2 = (BCSide)-1, const BCSide face3 = (BCSide)-1)
1281  {
1282  typedef typename SrcMskType::field_type SrcFieldType;
1283  typedef SpatialMask<DstFieldType> DstMskType;
1284  typedef typename DstFieldType::Location DstLocation;
1285  typedef typename SrcFieldType::Location SrcLocation;
1287 
1288  MaskConversionCheck<DstLocation, SrcLocation> conversion_check(face1);
1289  return ConvertedMask<const SrcMskType, const DstMskType, FDDir>(src, face1, face2, face3);
1290  }
1291 
1292  template <typename DstFieldType, typename SrcMskType>
1293  static inline ConvertedMask<const SrcMskType,
1295  typename GetFDOperatorDir<NODIR,
1296  DstFieldType,
1297  typename SrcMskType::field_type>::result
1298  > convert(const SrcMskType& src, const BCSide face1, const BCSide face2 = (BCSide)-1, const BCSide face3 = (BCSide)-1)
1299  {
1300  return convert<DstFieldType, NODIR, SrcMskType>(src, face1, face2, face3);
1301  }
1302 
1303  //TODO: move this to other files
1304  template<typename SrcT, typename DestT, typename Dir>
1305  static inline void masked_assign(
1306  const ConvertedMask<const SpatialMask<SrcT>,
1307  const SpatialMask<DestT>,
1308  Dir>& mask,
1309  DestT& field,
1310  const typename DestT::value_type& rhs)
1311  {
1312 #if __CUDACC__
1313  if(field.active_device_index() == CPU_INDEX)
1314 #endif
1315  masked_assign_imp(mask, field, NeboScalar<Initial, typename DestT::value_type>(rhs));
1316 #if __CUDACC__
1317  else
1318  field <<= cond(mask.mask(), rhs)(field);
1319 #endif
1320  }
1321 
1322  template<typename SrcT, typename DestT, typename Dir>
1323  static inline void masked_assign(
1324  const ConvertedMask<const SpatialMask<SrcT>,
1325  const SpatialMask<DestT>,
1326  Dir>& mask,
1327  DestT& field,
1328  const DestT& rhs)
1329  {
1330 #if __CUDACC__
1331  if(field.active_device_index() == CPU_INDEX)
1332 #endif
1333  masked_assign_imp(mask, field, NeboConstField<Initial, DestT>(rhs));
1334 #if __CUDACC__
1335  else
1336  field <<= cond(mask.mask(), rhs)(field);
1337 #endif
1338  }
1339 
1340  template<typename SrcT, typename DestT, typename Dir, typename RSHType, typename Expr>
1341  static inline void masked_assign(
1342  const ConvertedMask<const SpatialMask<SrcT>,
1343  const SpatialMask<DestT>,
1344  Dir>& mask,
1345  DestT& field,
1347  {
1348 #if __CUDACC__
1349  if(field.active_device_index() == CPU_INDEX)
1350 #endif
1351  masked_assign_imp(mask, field, rhs.expr());
1352 #if __CUDACC__
1353  else
1354  field <<= cond(mask.mask(), rhs)(field);
1355 #endif
1356 
1357  }
1358 
1359  template<typename SrcT, typename DestT, typename Dir, typename RHS>
1360  static inline void masked_assign_imp(
1361  const ConvertedMask<const SpatialMask<SrcT>,
1362  const SpatialMask<DestT>,
1363  Dir>& mask,
1364  DestT& field,
1365  const RHS& rhs)
1366  {
1367  typedef NeboField<Initial, DestT> LHSType;
1368  LHSType ilhs(field);
1369 #if __CUDACC__
1370  if(IS_GPU_INDEX(field.active_device_index())) {
1371  std::ostringstream msg;
1372  msg << "Nebo error in " << "Nebo Masked Assignment" << ":\n"
1373  ;
1374  msg << "Left-hand side of masked assignment allocated on ";
1375  msg << "GPU and this backend does not support GPU execution"
1376  ;
1377  msg << "\n";
1378  msg << "\t - " << __FILE__ << " : " << __LINE__;
1379  throw(std::runtime_error(msg.str()));
1380  }
1381  else {
1382  if(IS_CPU_INDEX(field.active_device_index())) {
1383  if(rhs.cpu_ready()) {
1384 #endif
1385  typename LHSType::SeqWalkType lhs = ilhs.init();
1386  typename RHS::SeqWalkType expr = rhs.init(IntVec(0,0,0),
1387  GhostData(0),
1388  IntVec(0,0,0));
1389 
1390 
1391  for(int i = 0; i <= mask.all(); i ++)
1392  {
1393  IntVec shift = mask.get_shift_vec(i);
1394 
1395  std::vector<IntVec>::const_iterator ip = mask.points().begin();
1396 
1397  std::vector<IntVec>::const_iterator const ep = mask.points().end();
1398 
1399  for(; ip != ep; ip++) {
1400  int const x = (*ip)[0] - shift[0];
1401 
1402  int const y = (*ip)[1] - shift[1];
1403 
1404  int const z = (*ip)[2] - shift[2];
1405 
1406  lhs.ref(x, y, z) = expr.template eval</*TODO: Wrap Into Compile Time Optional Args*/void>(x, y, z);
1407  };
1408  }
1409 #if __CUDACC__
1410  }
1411  else {
1412  std::ostringstream msg;
1413  msg << "Nebo error in " << "Nebo Assignment" << ":\n";
1414  msg << "Left-hand side of assignment allocated ";
1415  msg << "on ";
1416  msg << "CPU but right-hand side is not ";
1417  msg << "(completely) accessible on the same CPU";
1418  msg << "\n";
1419  msg << "\t - " << __FILE__ << " : " << __LINE__;
1420  throw(std::runtime_error(msg.str()));
1421  }
1422  }
1423  else {
1424  std::ostringstream msg;
1425  msg << "Nebo error in " << "Nebo Assignment" << ":\n";
1426  msg << "Left-hand side of assignment allocated ";
1427  msg << "on ";
1428  msg << "unknown device - not on CPU or GPU";
1429  msg << "\n";
1430  msg << "\t - " << __FILE__ << " : " << __LINE__;
1431  throw(std::runtime_error(msg.str()));
1432  }
1433  }
1434 #endif
1435  }
1436 
1473  template<typename DirT,
1474  typename SrcFieldT,
1475  typename DestFieldT,
1476  typename StencilBuilderT,
1477  typename NegativeStencilBuildersT,
1478  typename PositiveStencilBuildersT>
1480  public:
1481 
1482  template<typename PListT, typename NListT, typename dummy = void>
1484  typedef StencilListCheckSameField<typename PListT::Collection,
1485  typename NListT::Collection> StencilCollection;
1486  /* Source Field Check */
1487  typedef typename std::enable_if<std::is_same<typename PListT::First::SrcFieldType,
1488  typename StencilCollection::SrcFieldType>::value,
1489  SrcFieldT>::type SrcFieldType1; //Stencil in varying edge has different source field type
1490  typedef typename std::enable_if<std::is_same<typename NListT::First::SrcFieldType,
1491  typename StencilCollection::SrcFieldType>::value,
1492  SrcFieldT>::type SrcFieldType2; //Stencil in varying edge has different source field type
1493  typedef typename std::enable_if<std::is_same<SrcFieldType1, SrcFieldType2>::value,
1494  SrcFieldType1>::type SrcFieldType; //Stencil in varying edge has different source field type
1495  /* Destination Field Check */
1496  typedef typename std::enable_if<std::is_same<typename PListT::First::DestFieldType,
1497  typename StencilCollection::DestFieldType>::value,
1498  DestFieldT>::type DestFieldType1; //Stencil in varying edge has different destination field type
1499  typedef typename std::enable_if<std::is_same<typename NListT::First::DestFieldType,
1500  typename StencilCollection::DestFieldType>::value,
1501  DestFieldT>::type DestFieldType2; //Stencil in varying edge has different destination field type
1502  typedef typename std::enable_if<std::is_same<DestFieldType1, DestFieldType2>::value,
1503  DestFieldType1>::type DestFieldType; //Stencil in varying edge has different destination field type
1504  };
1505  template<typename dummy>
1507  typedef SrcFieldT SrcFieldType;
1508  typedef DestFieldT DestFieldType;
1509  };
1510  typedef typename StencilListCheckSameField<PositiveStencilBuildersT,
1511  NegativeStencilBuildersT>::SrcFieldType
1513  typedef typename StencilListCheckSameField<PositiveStencilBuildersT,
1514  NegativeStencilBuildersT>::DestFieldType
1516 
1517 
1518  //Conversion of lists from builder to stencil
1519  template<typename ListT, typename Arg>
1521  typedef typename ListT::First::template WithArg<Arg>::Stencil CurrentStencil;
1525  Result;
1526  static inline Result GetStencils(ListT const list, typename Arg::field_type const & field)
1527  {
1528  return StencilWithArgCollection::GetStencils(list.others(), field)((list.current()(field)).expr());
1529  }
1530  static inline Result GetStencils(ListT const list, NeboExpression<Arg, SrcFieldType> const & expr)
1531  {
1532  return StencilWithArgCollection::GetStencils(list.others(), expr)((list.current()(expr)).expr());
1533  }
1534  };
1535  template<typename Arg>
1538 
1539  static inline Result GetStencils(NeboGenericEmptyTypeList const list, Arg const & field)
1540  {
1541  return list;
1542  }
1543  static inline Result GetStencils(NeboGenericEmptyTypeList const list, NeboExpression<Arg, SrcFieldType> const & expr)
1544  {
1545  return list;
1546  }
1547  };
1548 
1549  // typedefs for when argument is a Nebo expression
1550  template<typename Arg>
1551  struct WithArg {
1552  typedef typename StencilBuilderT::template WithArg<Arg>::Stencil MainStencil;
1555  typedef NeboVaryingEdgeStencilOneDim<Initial,
1556  MainStencil,
1557  NegativeStencils,
1558  PositiveStencils,
1559  Arg,
1560  DirT,
1563  };
1564 
1565  // typedefs for when argument is a field
1567  typedef typename StencilBuilderT::template WithArg<FieldArg>::Stencil FieldMainStencil;
1570  typedef NeboVaryingEdgeStencilOneDim<Initial,
1571  FieldMainStencil,
1572  FieldNegativeStencils,
1573  FieldPositiveStencils,
1574  FieldArg,
1575  DirT,
1578 
1582  NeboVaryingEdgeStencilBuilder(const StencilBuilderT & mainST,
1583  const NegativeStencilBuildersT & negST,
1584  const PositiveStencilBuildersT & posST)
1585  : mainST_(mainST), negST_(negST), posST_(posST)
1586  {}
1587 
1589 
1595  void apply_to_field( const SrcFieldType & src, DestFieldType & dest ) const {
1596  dest <<= operator()(src);
1597  }
1598 
1603  inline FieldResult
1604  operator ()( const SrcFieldType & src ) const {
1605  typedef FieldResult FieldResult;
1606  typedef FieldStencil FieldStencil;
1609  FieldArg const arg(src);
1610  return FieldResult(FieldStencil(arg,
1611  mainST_(src).expr(),
1612  NegativeConverter::GetStencils(negST_, src),
1613  PositiveConverter::GetStencils(posST_, src)));
1614  }
1615 
1620  template<typename Arg>
1621  inline typename WithArg<Arg>::Result
1622  operator ()( const NeboExpression<Arg, SrcFieldType> & src ) const {
1623  typedef typename WithArg<Arg>::Stencil Stencil;
1624  typedef typename WithArg<Arg>::Result Result;
1625  typedef StencilListWithArg<NegativeStencilBuildersT, Arg> NegativeConverter;
1626  typedef StencilListWithArg<PositiveStencilBuildersT, Arg> PositiveConverter;
1627  return Result(Stencil(src.expr(),
1628  mainST_(src).expr(),
1629  NegativeConverter::GetStencils(negST_, src),
1630  PositiveConverter::GetStencils(posST_, src)));
1631  }
1632 
1633  private:
1634  const StencilBuilderT mainST_;
1635  const NegativeStencilBuildersT negST_;
1636  const PositiveStencilBuildersT posST_;
1637  };
1638 
1639 } // namespace SpatialOps
1640 #endif // NEBO_STENCIL_BUILDER_H
void operator()(SpatialMask< GammaFieldType > mask, PhiFieldType &phi, const GammaFieldType &gamma, bool minus) const
Apply boundary condition with gamma as a field.
NeboStencilBuilder(const CoefCollection &coefs)
construct a stencil with the specified coefficients
Provides tools to index into a sub-block of memory.
Definition: MemoryWindow.h:63
Supports definition of new Nebo stencils that do NOT invalidate ghost cells.
void apply_imp(MaskType mask, PhiFieldType &phi, const IntVec &shift, const NeboExpression< ExprType, GammaFieldType > &gamma, BCSide side) const
Apply boundary condition with gamma as an expression and a converted/non-converted mask...
CoefCollection const & coefs(void) const
Return coefficient collection.
OperatorT OperatorType
operator type
MaskShiftPoints< OperatorType, SrcFieldType, DestFieldType >::PlusPoint PlusPoint
positive face shift for mask
MinusResult minus(const SrcMask &src) const
Compute the minus side shift.
SrcFieldT SrcFieldType
source field type
Perform compile-time subtraction over a list of IndexTriplet types.
Definition: IndexTriplet.h:314
PntCltnT PointCollectionType
collection of stencil points
NeboBooleanExpression< PlusShift, DestFieldType > PlusResult
result type for positive shift
void operator()(SpatialMask< GammaFieldType > mask, PhiFieldType &phi, const NeboExpression< ExprType, GammaFieldType > &gamma, bool minus) const
Apply boundary condition with gamma as an expression.
PointCollection::AllButLast NonLowPoints
stencil offsets for all but low point in phi, assumes gamma is origin
Used for specifying field type traits.
Definition: IndexTriplet.h:121
void apply_to_field(const SrcFieldType &src, DestFieldType &dest) const
Apply this operator to the supplied source field to produce the supplied destination field...
Perform compile-time compare of two IndexTriplet types.
Definition: IndexTriplet.h:377
PntCltnT PointCollectionType
collection of stencil points
SrcFieldT SrcFieldType
source field type
void apply_to_field(const SrcFieldType &src, DestFieldType &dest) const
Apply this operator to the supplied source field to produce the supplied destination field...
Type for a Field Mask after it has been converted.
OperatorType type
operator type (Interpolant, Gradient, Divergence)
void apply_to_field(const SrcFieldType &src, DestFieldType &dest) const
Apply this operator to the supplied source field to produce the supplied destination field...
PlusResult plus(const SrcMask &src) const
Compute the plus side shift.
Defines a type to represent no direction.
OperatorType::PointCollectionType PointCollection
stencil point collection
NeboAverageStencilBuilder()
construct a stencil
NeboStencilCoefCollection< PointCollection::length > CoefCollection
collection of coefficients
NeboMaskShift< Initial, MinusPoint, Mask, DestFieldType > MinusShift
shift type for negative shift
Supports definition of new Nebo sum stencils, which sums given stencil points WITHOUT coefficients...
DestFieldT DestFieldType
destination field type
void apply_to_field(const SrcFieldType &src, DestFieldType &dest) const
Apply this operator to the supplied source field to produce the supplied destination field...
PointCollection::First HighPoint
stencil offset for high point in phi, assumes gamma is origin
SpatialMask< DestFieldType > DestMask
destination mask type
NeboBoundaryConditionBuilder(OperatorType const &op)
construct a boundary condition
NeboEdgelessStencilBuilder(const CoefCollection &coefs)
construct a stencil with the specified coefficients
for wide stencils where we set on a point rather than a face
SrcFieldT SrcFieldType
source field type
NeboMask< Initial, SrcFieldType > Mask
Nebo mask type.
void operator()(ConvertedPhiMask mask, PhiFieldType &phi, const NeboExpression< ExprType, GammaFieldType > &gamma) const
Apply boundary condition with gamma as an expression and a converted mask as boundary mask...
DestFieldT DestFieldType
destination field type
NeboStencilCoefCollection< PointCollectionType::length > CoefCollection
collection of coefficients
DestFieldT DestFieldType
destination field type
NeboSumStencilBuilder()
construct a stencil
OperatorType type
operator type (Interpolant, Gradient, Divergence)
Supports definition of new Nebo stencils.
Abstracts a mask.
Definition: SpatialMask.h:70
NeboVaryingEdgeStencilBuilder(const StencilBuilderT &mainST, const NegativeStencilBuildersT &negST, const PositiveStencilBuildersT &posST)
construct a stencil with the specified coefficients
Perform compile-time compare of two IndexTriplet types.
Definition: IndexTriplet.h:356
OperatorType::DestFieldType GammaFieldType
type of fields in gamma, which this operator reads
NeboBooleanExpression< MinusShift, DestFieldType > MinusResult
result type for negative shift
PointCollection::AllButFirst NonHighPoints
stencil offsets for all but high point in phi, assumes gamma is origin
OperatorType::SrcFieldType PhiFieldType
field type of phi, which this operator modifies
IntVec get_plus() const
obtain the IntVec containing the number of ghost cells on the (+) faces
Definition: GhostData.h:145
PntCltnT PointCollectionType
collection of stencil points
NeboMaskShift< Initial, PlusPoint, Mask, DestFieldType > PlusShift
shift type for positive shift
void operator()(ConvertedPhiMask mask, PhiFieldType &phi, const double gamma) const
Apply boundary condition with gamma as a scalar.
Holds information about the number of ghost cells on each side of the domain.
Definition: GhostData.h:54
The Nebo Operator for switch.
Supports definition of new Nebo mask shift stencils, which converts a mask of one field type into ano...
Supports definition of new Nebo boundary condition.
SpatialMask< SrcFieldType > SrcMask
source mask type
SrcFieldT SrcFieldType
source field type
BCSide
Allows identification of whether we are setting the BC on the right or left side when using an operat...
NeboMaskShiftBuilder()
construct a stencil
PntCltnT PointCollectionType
collection of stencil points
PointCollection::Last LowPoint
stencil offset for low point in phi, assumes gamma is origin
StencilListCheckSameField< PositiveStencilBuildersT, NegativeStencilBuildersT >::DestFieldType DestFieldType
destination field type
void operator()(const SpatialMask< GammaFieldType > mask, PhiFieldType &phi, const double gamma, bool minus) const
Apply boundary condition with gamma as a scalar.
StencilListCheckSameField< PositiveStencilBuildersT, NegativeStencilBuildersT >::SrcFieldType SrcFieldType
source field type
Supports definition of new Nebo stencils that vary as they approach domain edges. ...
Supports definition of new Nebo average stencils, which automatically averages given stencil points...
NeboStencilCoefCollection< PointCollectionType::length > CoefCollection
collection of coefficients
SrcFieldT SrcFieldType
source field type
Abstracts a field.
Definition: SpatialField.h:132
void apply_to_field(const SrcFieldType &src, DestFieldType &dest) const
Apply this operator to the supplied source field to produce the supplied destination field...
IntVec get_minus() const
obtain the IntVec containing the number of ghost cells on the (-) faces
Definition: GhostData.h:135
void operator()(ConvertedPhiMask mask, PhiFieldType &phi, const GammaFieldType &gamma) const
Apply boundary condition with gamma as a field.
DestFieldT DestFieldType
destination field type
OperatorT OperatorType
type of operator to invert
CoefCollection const & coefs(void) const
Return coefficient collection.
MaskShiftPoints< OperatorType, SrcFieldType, DestFieldType >::MinusPoint MinusPoint
negative face shift for mask
DestFieldT DestFieldType
destination field type