OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
IntTSC.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /***************************************************************************
3  *
4  * The IPPL Framework
5  *
6  *
7  * Visit http://people.web.psi.ch/adelmann/ for more details
8  *
9  ***************************************************************************/
10 
11 #ifndef INT_TSC_H
12 #define INT_TSC_H
13 
14 /* IntTSC.h -- Definition of simple class to perform cloud-in-cell
15  interpolation of data for a single particle to or from a IPPL Field.
16  This interpolation scheme is also referred to as linear interpolation,
17  area-weighting (in 2D), or volume-weighting (in 3D). */
18 
19 // include files
20 #include "Particle/Interpolator.h"
21 #include "Field/Field.h"
22 
23 
24 // forward declaration
25 class IntTSC;
26 
27 // specialization of InterpolatorTraits
28 template <class T, unsigned Dim>
29 struct InterpolatorTraits<T,Dim,IntTSC> {
30  typedef CacheDataTSC<T,Dim> Cache_t;
31 };
32 
33 
34 
35 // IntTSCImpl class definition
36 template <unsigned Dim>
37 class IntTSCImpl : public Interpolator {
38 
39 public:
40  // constructor/destructor
41  IntTSCImpl() {}
42  ~IntTSCImpl() {}
43 
44  // gather/scatter functions
45 
46  // scatter particle data into Field using particle position and mesh
47  template <class FT, class M, class C, class PT>
48  static
49  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
50  const Vektor<PT,Dim>& ppos, const M& mesh) {
51  CenteringTag<C> ctag;
52  Vektor<PT,Dim> gpos, dpos, delta;
53  NDIndex<Dim> ngp;
54  CompressedBrickIterator<FT,Dim> fiter;
55  //BENI: offset not needed, since NGP is used as center
56  //int lgpoff[Dim];
57  // find nearest grid point for particle position, store in NDIndex obj
58  ngp = FindNGP(mesh, ppos, ctag);
59  // get position of ngp
60  FindPos(gpos, mesh, ngp, ctag);
61  // get mesh spacings
62  FindDelta(delta, mesh, ngp, ctag);
63  // Try to find ngp in local fields and get iterator
64  fiter = getFieldIter(f,ngp);
65  // Now we find the offset from ngp to next-lowest grip point (lgp)
66  /* //BENI: we take nearest grid point as it is as the center. No modification needed.
67  for (d=0; d<Dim; ++d) {
68  if (gpos(d) > ppos(d)) {
69  lgpoff[d] = -1; // save the offset
70  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
71  }
72  else {
73  lgpoff[d] = 0; // save the offset
74  }
75  }
76  // adjust position of Field iterator to lgp position
77  fiter.moveBy(lgpoff);
78  */
79  // get distance from ppos to gpos
80  dpos = ppos - gpos;
81  // normalize dpos by mesh spacing
82  dpos /= delta;
83  // accumulate into local elements
84  ERRORMSG("IntTSC::scatter: not implemented for Dim>3!!"<<endl);
85  return;
86  }
87 
88  // scatter particle data into Field using particle position and mesh
89  // and cache mesh information for reuse
90  template <class FT, class M, class C, class PT>
91  static
92  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
93  const Vektor<PT,Dim>& ppos, const M& mesh,
94  NDIndex<Dim>& ngp, int lgpoff[Dim], Vektor<PT,Dim>& dpos) {
95  CenteringTag<C> ctag;
96  Vektor<PT,Dim> gpos, delta;
97  CompressedBrickIterator<FT,Dim> fiter;
98  // find nearest grid point for particle position, store in NDIndex obj
99  ngp = FindNGP(mesh, ppos, ctag);
100  // get position of ngp
101  FindPos(gpos, mesh, ngp, ctag);
102  // get mesh spacings
103  FindDelta(delta, mesh, ngp, ctag);
104  // Try to find ngp in local fields and get iterator
105  fiter = getFieldIter(f,ngp);
106  // Now we find the offset from ngp to next-lowest grip point (lgp)
107  /*
108  for (d=0; d<Dim; ++d) {
109  if (gpos(d) > ppos(d)) {
110  lgpoff[d] = -1; // save the offset
111  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
112  }
113  else {
114  lgpoff[d] = 0; // save the offset
115  }
116  }
117  // adjust position of Field iterator to lgp position
118  fiter.moveBy(lgpoff);
119  // get distance from ppos to gpos
120  */
121  dpos = ppos - gpos;
122  // normalize dpos by mesh spacing
123  dpos /= delta;
124  // accumulate into local elements
125  ERRORMSG("IntTSC::scatter: not implemented for Dim>3!!"<<endl);
126  return;
127  }
128 
129  // scatter particle data into Field using cached mesh information
130  template <class FT, class M, class C, class PT>
131  static
132  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
133  const NDIndex<Dim>& ngp, const int lgpoff[Dim],
134  const Vektor<PT,Dim>& dpos) {
135  CompressedBrickIterator<FT,Dim> fiter;
136  // Try to find ngp in local fields and get iterator
137  fiter = getFieldIter(f,ngp);
138  // adjust position of Field iterator to lgp position
139  //fiter.moveBy(lgpoff);
140  // accumulate into local elements
141  ERRORMSG("IntTSC::scatter: not implemented for Dim>3!!"<<endl);
142  return;
143  }
144 
145  // gather particle data from Field using particle position and mesh
146  template <class FT, class M, class C, class PT>
147  static
148  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
149  const Vektor<PT,Dim>& ppos, const M& mesh) {
150  CenteringTag<C> ctag;
151  Vektor<PT,Dim> gpos, dpos, delta;
152  NDIndex<Dim> ngp;
153  CompressedBrickIterator<FT,Dim> fiter;
154  int lgpoff[Dim];
155  // find nearest grid point for particle position, store in NDIndex obj
156  ngp = FindNGP(mesh, ppos, ctag);
157  // get position of ngp
158  FindPos(gpos, mesh, ngp, ctag);
159  // get mesh spacings
160  FindDelta(delta, mesh, ngp, ctag);
161  // Try to find ngp in local fields and get iterator
162  fiter = getFieldIter(f,ngp);
163  /*
164  // Now we find the offset from ngp to next-lowest grip point (lgp)
165  for (d=0; d<Dim; ++d) {
166  if (gpos(d) > ppos(d)) {
167  lgpoff[d] = -1; // save the offset
168  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
169  }
170  else {
171  lgpoff[d] = 0; // save the offset
172  }
173  }
174  // adjust position of Field iterator to lgp position
175  fiter.moveBy(lgpoff);
176  */
177  // get distance from ppos to gpos
178  dpos = ppos - gpos;
179  // normalize dpos by mesh spacing
180  dpos /= delta;
181  // accumulate into particle attrib
182  ERRORMSG("IntTSC::gather: not implemented for Dim>3!!"<<endl);
183  return;
184  }
185 
186  // gather particle data from Field using particle position and mesh
187  // and cache mesh information for reuse
188  template <class FT, class M, class C, class PT>
189  static
190  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
191  const Vektor<PT,Dim>& ppos, const M& mesh,
192  NDIndex<Dim>& ngp, int lgpoff[Dim], Vektor<PT,Dim>& dpos) {
193  CenteringTag<C> ctag;
194  Vektor<PT,Dim> gpos, delta;
195  CompressedBrickIterator<FT,Dim> fiter;
196  // find nearest grid point for particle position, store in NDIndex obj
197  ngp = FindNGP(mesh, ppos, ctag);
198  // get position of ngp
199  FindPos(gpos, mesh, ngp, ctag);
200  // get mesh spacings
201  FindDelta(delta, mesh, ngp, ctag);
202  // Try to find ngp in local fields and get iterator
203  fiter = getFieldIter(f,ngp);
204  /*
205  // Now we find the offset from ngp to next-lowest grip point (lgp)
206  for (d=0; d<Dim; ++d) {
207  if (gpos(d) > ppos(d)) {
208  lgpoff[d] = -1; // save the offset
209  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
210  }
211  else {
212  lgpoff[d] = 0; // save the offset
213  }
214  }
215  // adjust position of Field iterator to lgp position
216  fiter.moveBy(lgpoff);
217  */
218  // get distance from ppos to gpos
219  dpos = ppos - gpos;
220  // normalize dpos by mesh spacing
221  dpos /= delta;
222  // accumulate into particle attrib
223  ERRORMSG("IntTSC::gather: not implemented for Dim>3!!"<<endl);
224  return;
225  }
226 
227  // gather particle data from Field using cached mesh information
228  template <class FT, class M, class C, class PT>
229  static
230  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
231  const NDIndex<Dim>& ngp, const int lgpoff[Dim],
232  const Vektor<PT,Dim>& dpos) {
233  CompressedBrickIterator<FT,Dim> fiter;
234  // Try to find ngp in local fields and get iterator
235  fiter = getFieldIter(f,ngp);
236  // adjust position of Field iterator to lgp position
237  //fiter.moveBy(lgpoff);
238  // accumulate into particle attrib
239  ERRORMSG("IntTSC::gather: not implemented for Dim>3!!"<<endl);
240  return;
241  }
242 
243 };
244 
245 
246 template <>
247 class IntTSCImpl<1U> : public Interpolator {
248 
249 public:
250  // constructor/destructor
251  IntTSCImpl() {}
252  ~IntTSCImpl() {}
253 
254  // gather/scatter functions
255 
256  // scatter particle data into Field using particle position and mesh
257  template <class FT, class M, class C, class PT>
258  static
259  void scatter(const FT& pdata, Field<FT,1U,M,C>& f,
260  const Vektor<PT,1U>& ppos, const M& mesh) {
261  CenteringTag<C> ctag;
262  Vektor<PT,1U> gpos, dpos, delta;
263  NDIndex<1U> ngp;
264  CompressedBrickIterator<FT,1U> fiter;
265  // find nearest grid point for particle position, store in NDIndex obj
266  ngp = FindNGP(mesh, ppos, ctag);
267  // get position of ngp
268  FindPos(gpos, mesh, ngp, ctag);
269  // get mesh spacings
270  FindDelta(delta, mesh, ngp, ctag);
271  // Try to find ngp in local fields and get iterator
272  fiter = getFieldIter(f,ngp);
273  // get distance from ppos to gpos
274  dpos = ppos - gpos;
275  // normalize dpos by mesh spacing
276  dpos /= delta;
277  // accumulate into local elements
278  /*
279  *fiter += .25*(3.-4.dpos(0)*dpos(0)) * pdata;
280  fiter.offset(-1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * pdata;
281  fiter.offset(+1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * pdata;
282  */
283  auto W = [dpos](unsigned p, unsigned i) {
284  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
285  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
286  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
287 
288  for (int p0=-1; p0<=1; ++p0) {
289  fiter.offset(p0) += W(p0,0) * pdata;
290  }
291  return;
292  }
293 
294  // scatter particle data into Field using particle position and mesh
295  // and cache mesh information for reuse
296  template <class FT, class M, class C, class PT>
297  static
298  void scatter(const FT& pdata, Field<FT,1U,M,C>& f,
299  const Vektor<PT,1U>& ppos, const M& mesh,
300  NDIndex<1U>& ngp, int lgpoff[1U], Vektor<PT,1U>& dpos) {
301  CenteringTag<C> ctag;
302  Vektor<PT,1U> gpos, delta;
303  CompressedBrickIterator<FT,1U> fiter;
304  // find nearest grid point for particle position, store in NDIndex obj
305  ngp = FindNGP(mesh, ppos, ctag);
306  // get position of ngp
307  FindPos(gpos, mesh, ngp, ctag);
308  // get mesh spacings
309  FindDelta(delta, mesh, ngp, ctag);
310  // Try to find ngp in local fields and get iterator
311  fiter = getFieldIter(f,ngp);
312  // get distance from ppos to gpos
313  dpos = ppos - gpos;
314  // normalize dpos by mesh spacing
315  dpos /= delta;
316  // accumulate into local elements
317  /*
318  *fiter += .25*(3.-4.dpos(0)*dpos(0)) * pdata;
319  fiter.offset(-1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * pdata;
320  fiter.offset(+1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * pdata;
321  */
322  auto W = [dpos](unsigned p, unsigned i) {
323  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
324  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
325  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
326 
327  for (int p0=-1; p0<=1; ++p0) {
328  fiter.offset(p0) += W(p0,0) * pdata;
329  }
330 
331  return;
332  }
333 
334  // scatter particle data into Field using cached mesh information
335  template <class FT, class M, class C, class PT>
336  static
337  void scatter(const FT& pdata, Field<FT,1U,M,C>& f,
338  const NDIndex<1U>& ngp, const int lgpoff[1U],
339  const Vektor<PT,1U>& dpos) {
340  CompressedBrickIterator<FT,1U> fiter;
341  // Try to find ngp in local fields and get iterator
342  fiter = getFieldIter(f,ngp);
343  // adjust position of Field iterator to lgp position
344  //fiter.moveBy(lgpoff);
345  // accumulate into local elements
346  /*
347  *fiter += .25*(3.-4.dpos(0)*dpos(0)) * pdata;
348  fiter.offset(-1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * pdata;
349  fiter.offset(+1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * pdata;
350  */
351  auto W = [dpos](unsigned p, unsigned i) {
352  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
353  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
354  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
355 
356  for (int p0=-1; p0<=1; ++p0) {
357  fiter.offset(p0) += W(p0,0) * pdata;
358  }
359 
360  return;
361  }
362 
363  // gather particle data from Field using particle position and mesh
364  template <class FT, class M, class C, class PT>
365  static
366  void gather(FT& pdata, const Field<FT,1U,M,C>& f,
367  const Vektor<PT,1U>& ppos, const M& mesh) {
368  CenteringTag<C> ctag;
369  Vektor<PT,1U> gpos, dpos, delta;
370  NDIndex<1U> ngp;
371  CompressedBrickIterator<FT,1U> fiter;
372  // find nearest grid point for particle position, store in NDIndex obj
373  ngp = FindNGP(mesh, ppos, ctag);
374  // get position of ngp
375  FindPos(gpos, mesh, ngp, ctag);
376  // get mesh spacings
377  FindDelta(delta, mesh, ngp, ctag);
378  // Try to find ngp in local fields and get iterator
379  fiter = getFieldIter(f,ngp);
380  // get distance from ppos to gpos
381  dpos = ppos - gpos;
382  // normalize dpos by mesh spacing
383  dpos /= delta;
384  // accumulate into particle attrib
385  /* pdata = (1 - dpos(0)) * (*fiter) +
386  dpos(0) * fiter.offset(1);*/
387  /*
388  pdata = .25*(3.-4.dpos(0)*dpos(0)) * (*fiter) +
389  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * fiter.offset(-1) +
390  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * fiter.offset(+1);
391 */
392  auto W = [dpos](unsigned p, unsigned i) {
393  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
394  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
395  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
396 
397  pdata = 0;
398  for (int p0=-1; p0<=1; ++p0) {
399  pdata += W(p0,0)*fiter.offset(p0);
400  }
401 
402  return;
403  }
404 
405  // gather particle data from Field using particle position and mesh
406  // and cache mesh information for reuse
407  template <class FT, class M, class C, class PT>
408  static
409  void gather(FT& pdata, const Field<FT,1U,M,C>& f,
410  const Vektor<PT,1U>& ppos, const M& mesh,
411  NDIndex<1U>& ngp, int lgpoff[1U], Vektor<PT,1U>& dpos) {
412  CenteringTag<C> ctag;
413  Vektor<PT,1U> gpos, delta;
414  CompressedBrickIterator<FT,1U> fiter;
415  // find nearest grid point for particle position, store in NDIndex obj
416  ngp = FindNGP(mesh, ppos, ctag);
417  // get position of ngp
418  FindPos(gpos, mesh, ngp, ctag);
419  // get mesh spacings
420  FindDelta(delta, mesh, ngp, ctag);
421  // Try to find ngp in local fields and get iterator
422  fiter = getFieldIter(f,ngp);
423  // get distance from ppos to gpos
424  dpos = ppos - gpos;
425  // normalize dpos by mesh spacing
426  dpos /= delta;
427  // accumulate into particle attrib
428  /*
429  pdata = .25*(3.-4.dpos(0)*dpos(0)) * (*fiter) +
430  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * fiter.offset(-1) +
431  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * fiter.offset(+1);
432  */
433  auto W = [dpos](unsigned p, unsigned i) {
434  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
435  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
436  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
437 
438  pdata = 0;
439  for (int p0=-1; p0<=1; ++p0) {
440  pdata += W(p0,0)*fiter.offset(p0);
441  }
442 
443 
444  return;
445  }
446 
447  // gather particle data from Field using cached mesh information
448  template <class FT, class M, class C, class PT>
449  static
450  void gather(FT& pdata, const Field<FT,1U,M,C>& f,
451  const NDIndex<1U>& ngp, const int lgpoff[1U],
452  const Vektor<PT,1U>& dpos) {
453  CompressedBrickIterator<FT,1U> fiter;
454  // Try to find ngp in local fields and get iterator
455  fiter = getFieldIter(f,ngp);
456  // accumulate into particle attrib
457  /*
458  pdata = .25*(3.-4.dpos(0)*dpos(0)) * (*fiter) +
459  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * fiter.offset(-1) +
460  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * fiter.offset(+1);
461  */
462  auto W = [dpos](unsigned p, unsigned i) {
463  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
464  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
465  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
466 
467  pdata = 0;
468  for (int p0=-1; p0<=1; ++p0) {
469  pdata += W(p0,0)*fiter.offset(p0);
470  }
471 
472  return;
473  }
474 
475 };
476 
477 
478 template <>
479 class IntTSCImpl<2U> : public Interpolator {
480 
481 public:
482  // constructor/destructor
483  IntTSCImpl() {}
484  ~IntTSCImpl() {}
485 
486  // gather/scatter functions
487 
488  // scatter particle data into Field using particle position and mesh
489  template <class FT, class M, class C, class PT>
490  static
491  void scatter(const FT& pdata, Field<FT,2U,M,C>& f,
492  const Vektor<PT,2U>& ppos, const M& mesh) {
493  CenteringTag<C> ctag;
494  Vektor<PT,2U> gpos, dpos, delta;
495  NDIndex<2U> ngp;
496  CompressedBrickIterator<FT,2U> fiter;
497  // find nearest grid point for particle position, store in NDIndex obj
498  ngp = FindNGP(mesh, ppos, ctag);
499  // get position of ngp
500  FindPos(gpos, mesh, ngp, ctag);
501  // get mesh spacings
502  FindDelta(delta, mesh, ngp, ctag);
503  // Try to find ngp in local fields and get iterator
504  fiter = getFieldIter(f,ngp);
505  // get distance from ppos to gpos
506  dpos = ppos - gpos;
507  // normalize dpos by mesh spacing
508  dpos /= delta;
509  // accumulate into local elements
510  /*
511  *fiter += .25*(3.-4.dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
512  fiter.offset(-1,-1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
513  fiter.offset(-1,0) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
514  fiter.offset(-1,+1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
515  fiter.offset(0,-1) += .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
516  fiter.offset(0,+1) += .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
517  fiter.offset(+1,-1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
518  fiter.offset(+1,0) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
519  fiter.offset(+1,+1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
520 */
521  auto W = [dpos](unsigned p, unsigned i) {
522  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
523  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
524  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
525 
526  for (int p0=-1; p0<=1; ++p0) {
527  for (int p1=-1; p1<=1; ++p1) {
528  fiter.offset(p0,p1) += W(p0,0) * W(p1,1) * pdata;
529  }
530  }
531  return;
532  }
533 
534  // scatter particle data into Field using particle position and mesh
535  // and cache mesh information for reuse
536  template <class FT, class M, class C, class PT>
537  static
538  void scatter(const FT& pdata, Field<FT,2U,M,C>& f,
539  const Vektor<PT,2U>& ppos, const M& mesh,
540  NDIndex<2U>& ngp, int lgpoff[2U], Vektor<PT,2U>& dpos) {
541  CenteringTag<C> ctag;
542  Vektor<PT,2U> gpos, delta;
543  CompressedBrickIterator<FT,2U> fiter;
544  // find nearest grid point for particle position, store in NDIndex obj
545  ngp = FindNGP(mesh, ppos, ctag);
546  // get position of ngp
547  FindPos(gpos, mesh, ngp, ctag);
548  // get mesh spacings
549  FindDelta(delta, mesh, ngp, ctag);
550  // Try to find ngp in local fields and get iterator
551  fiter = getFieldIter(f,ngp);
552  // get distance from ppos to gpos
553  dpos = ppos - gpos;
554  // normalize dpos by mesh spacing
555  dpos /= delta;
556  // accumulate into local elements
557  /*
558  *fiter += .25*(3.-4.dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
559  fiter.offset(-1,-1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
560  fiter.offset(-1,0) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
561  fiter.offset(-1,+1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
562  fiter.offset(0,-1) += .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
563  fiter.offset(0,+1) += .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
564  fiter.offset(+1,-1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
565  fiter.offset(+1,0) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
566  fiter.offset(+1,+1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
567 */
568 
569  auto W = [dpos](unsigned p, unsigned i) {
570  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
571  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
572  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
573 
574  for (int p0=-1; p0<=1; ++p0) {
575  for (int p1=-1; p1<=1; ++p1) {
576  fiter.offset(p0,p1) += W(p0,0) * W(p1,1) * pdata;
577  }
578  }
579 
580  return;
581  }
582 
583  // scatter particle data into Field using cached mesh information
584  template <class FT, class M, class C, class PT>
585  static
586  void scatter(const FT& pdata, Field<FT,2U,M,C>& f,
587  const NDIndex<2U>& ngp, const int lgpoff[2U],
588  const Vektor<PT,2U>& dpos) {
589  CompressedBrickIterator<FT,2U> fiter;
590  // Try to find ngp in local fields and get iterator
591  fiter = getFieldIter(f,ngp);
592  // adjust position of Field iterator to lgp position
593  // fiter.moveBy(lgpoff);
594  // accumulate into local elements
595  /*
596  *fiter += .25*(3.-4.dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
597  fiter.offset(-1,-1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
598  fiter.offset(-1,0) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
599  fiter.offset(-1,+1) += .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
600  fiter.offset(0,-1) += .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
601  fiter.offset(0,+1) += .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
602  fiter.offset(+1,-1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
603  fiter.offset(+1,0) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * pdata;
604  fiter.offset(+1,+1) += .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * pdata;
605 */
606 
607  auto W = [dpos](unsigned p, unsigned i) {
608  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
609  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
610  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
611 
612  for (int p0=-1; p0<=1; ++p0) {
613  for (int p1=-1; p1<=1; ++p1) {
614  fiter.offset(p0,p1) += W(p0,0) * W(p1,1) * pdata;
615  }
616  }
617 
618  return;
619  }
620 
621  // gather particle data from Field using particle position and mesh
622  template <class FT, class M, class C, class PT>
623  static
624  void gather(FT& pdata, const Field<FT,2U,M,C>& f,
625  const Vektor<PT,2U>& ppos, const M& mesh) {
626  CenteringTag<C> ctag;
627  Vektor<PT,2U> gpos, dpos, delta;
628  NDIndex<2U> ngp;
629  CompressedBrickIterator<FT,2U> fiter;
630  // find nearest grid point for particle position, store in NDIndex obj
631  ngp = FindNGP(mesh, ppos, ctag);
632  // get position of ngp
633  FindPos(gpos, mesh, ngp, ctag);
634  // get mesh spacings
635  FindDelta(delta, mesh, ngp, ctag);
636  // Try to find ngp in local fields and get iterator
637  fiter = getFieldIter(f,ngp);
638  // get distance from ppos to gpos
639  dpos = ppos - gpos;
640  // normalize dpos by mesh spacing
641  dpos /= delta;
642  // accumulate into particle attrib
643  /*
644  pdata = .25*(3.-4.dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * (*fiter) +
645  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(-1,-1) +
646  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * fiter.offset(-1,0) +
647  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(-1,+1) +
648  .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(0,-1) +
649  .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(0,+1) +
650  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) *fiter.offset(+1,-1) +
651  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * fiter.offset(+1,0) +
652  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(+1,+1) ;
653  */
654  pdata = 0;
655  auto W = [dpos](unsigned p, unsigned i) {
656  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
657  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
658  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
659 
660  for (int p0=-1; p0<=1; ++p0) {
661  for (int p1=-1; p1<=1; ++p1) {
662  pdata += W(p0,0) * W(p1,1) * fiter.offset(p0,p1);
663  }
664  }
665 
666  return;
667  }
668 
669  // gather particle data from Field using particle position and mesh
670  // and cache mesh information for reuse
671  template <class FT, class M, class C, class PT>
672  static
673  void gather(FT& pdata, const Field<FT,2U,M,C>& f,
674  const Vektor<PT,2U>& ppos, const M& mesh,
675  NDIndex<2U>& ngp, int lgpoff[2U], Vektor<PT,2U>& dpos) {
676  CenteringTag<C> ctag;
677  Vektor<PT,2U> gpos, delta;
678  CompressedBrickIterator<FT,2U> fiter;
679  // find nearest grid point for particle position, store in NDIndex obj
680  ngp = FindNGP(mesh, ppos, ctag);
681  // get position of ngp
682  FindPos(gpos, mesh, ngp, ctag);
683  // get mesh spacings
684  FindDelta(delta, mesh, ngp, ctag);
685  // Try to find ngp in local fields and get iterator
686  fiter = getFieldIter(f,ngp);
687  // get distance from ppos to gpos
688  dpos = ppos - gpos;
689  // normalize dpos by mesh spacing
690  dpos /= delta;
691  // accumulate into particle attrib
692  /*
693  pdata = .25*(3.-4.dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * (*fiter) +
694  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(-1,-1) +
695  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * fiter.offset(-1,0) +
696  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(-1,+1) +
697  .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(0,-1) +
698  .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(0,+1) +
699  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) *fiter.offset(+1,-1) +
700  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * fiter.offset(+1,0) +
701  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(+1,+1) ;
702  */
703 
704  pdata = 0;
705  auto W = [dpos](unsigned p, unsigned i) {
706  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
707  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
708  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
709 
710  for (int p0=-1; p0<=1; ++p0) {
711  for (int p1=-1; p1<=1; ++p1) {
712  pdata += W(p0,0) * W(p1,1) * fiter.offset(p0,p1);
713  }
714  }
715 
716  return;
717  }
718 
719  // gather particle data from Field using cached mesh information
720  template <class FT, class M, class C, class PT>
721  static
722  void gather(FT& pdata, const Field<FT,2U,M,C>& f,
723  const NDIndex<2U>& ngp, const int lgpoff[2U],
724  const Vektor<PT,2U>& dpos) {
725  CompressedBrickIterator<FT,2U> fiter;
726  // Try to find ngp in local fields and get iterator
727  fiter = getFieldIter(f,ngp);
728  // accumulate into particle attrib
729  /*
730  pdata = .25*(3.-4.dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * (*fiter) +
731  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(-1,-1) +
732  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * fiter.offset(-1,0) +
733  .125*(1.-4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(-1,+1) +
734  .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(0,-1) +
735  .25*(3.-4.dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(0,+1) +
736  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.-4.*dpos(1)+4.*dpos(1)*dpos(1)) *fiter.offset(+1,-1) +
737  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .25*(3.-4.dpos(1)*dpos(1)) * fiter.offset(+1,0) +
738  .125*(1.+4.*dpos(0)+4.*dpos(0)*dpos(0)) * .125*(1.+4.*dpos(1)+4.*dpos(1)*dpos(1)) * fiter.offset(+1,+1) ;
739  */
740 
741  pdata = 0;
742  auto W = [dpos](unsigned p, unsigned i) {
743  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
744  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
745  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
746 
747  for (int p0=-1; p0<=1; ++p0) {
748  for (int p1=-1; p1<=1; ++p1) {
749  pdata += W(p0,0) * W(p1,1) * fiter.offset(p0,p1);
750  }
751  }
752 
753  return;
754  }
755 
756 };
757 
758 
759 template <>
760 class IntTSCImpl<3U> : public Interpolator {
761 
762 public:
763  // constructor/destructor
764  IntTSCImpl() {}
765  ~IntTSCImpl() {}
766 
767  // gather/scatter functions
768 
769  // scatter particle data into Field using particle position and mesh
770  template <class FT, class M, class C, class PT>
771  static
772  void scatter(const FT& pdata, Field<FT,3U,M,C>& f,
773  const Vektor<PT,3U>& ppos, const M& mesh) {
774  CenteringTag<C> ctag;
775  Vektor<PT,3U> gpos, dpos, delta;
776  NDIndex<3U> ngp;
777  CompressedBrickIterator<FT,3U> fiter;
778  //int lgpoff[3U];
779  //unsigned int d;
780  // find nearest grid point for particle position, store in NDIndex obj
781  ngp = FindNGP(mesh, ppos, ctag);
782  // get position of ngp
783  FindPos(gpos, mesh, ngp, ctag);
784  // get mesh spacings
785  FindDelta(delta, mesh, ngp, ctag);
786  // Try to find ngp in local fields and get iterator
787  fiter = getFieldIter(f,ngp);
788  // get distance from ppos to gpos
789  dpos = ppos - gpos;
790  // normalize dpos by mesh spacing
791  dpos /= delta;
792  // accumulate into local elements
793  auto W = [dpos](unsigned p, unsigned i) {
794  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
795  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
796  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
797 
798  //double interpol_tot = 0;
799  for (int p0=-1; p0<=1; ++p0) {
800  for (int p1=-1; p1<=1; ++p1) {
801  for (int p2=-1; p2<=1; ++p2) {
802  fiter.offset(p0,p1,p2) += W(p0,0) * W(p1,1) * W(p2,2) * pdata;
803  //interpol_tot+= W(p0,0) * W(p1,1) * W(p2,2) * pdata;
804  }
805  }
806  }
807  //std::cout << "the total interpolation is = " << interpol_tot << std::endl;
808  return;
809  }
810 
811  // scatter particle data into Field using particle position and mesh
812  // and cache mesh information for reuse
813  template <class FT, class M, class C, class PT>
814  static
815  void scatter(const FT& pdata, Field<FT,3U,M,C>& f,
816  const Vektor<PT,3U>& ppos, const M& mesh,
817  NDIndex<3U>& ngp, int lgpoff[3U], Vektor<PT,3U>& dpos) {
818  CenteringTag<C> ctag;
819  Vektor<PT,3U> gpos, delta;
820  CompressedBrickIterator<FT,3U> fiter;
821  // find nearest grid point for particle position, store in NDIndex obj
822  ngp = FindNGP(mesh, ppos, ctag);
823  // get position of ngp
824  FindPos(gpos, mesh, ngp, ctag);
825  // get mesh spacings
826  FindDelta(delta, mesh, ngp, ctag);
827  // Try to find ngp in local fields and get iterator
828  fiter = getFieldIter(f,ngp);
829  // get distance from ppos to gpos
830  dpos = ppos - gpos;
831  // normalize dpos by mesh spacing
832  dpos /= delta;
833  // accumulate into local elements
834  auto W = [dpos](unsigned p, unsigned i) {
835  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
836  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
837  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
838 
839  //double interpol_tot = 0;
840  for (int p0=-1; p0<=1; ++p0) {
841  for (int p1=-1; p1<=1; ++p1) {
842  for (int p2=-1; p2<=1; ++p2) {
843  fiter.offset(p0,p1,p2) += W(p0,0) * W(p1,1) * W(p2,2) * pdata;
844  //interpol_tot+= W(p0,0) * W(p1,1) * W(p2,2) * pdata;
845  }
846  }
847  }
848 
849 // std::cout << "the total interpolation is = " << interpol_tot << std::endl;
850  return;
851  }
852 
853  // scatter particle data into Field using cached mesh information
854  template <class FT, class M, class C, class PT>
855  static
856  void scatter(const FT& pdata, Field<FT,3U,M,C>& f,
857  const NDIndex<3U>& ngp, const int lgpoff[3U],
858  const Vektor<PT,3U>& dpos) {
859  CompressedBrickIterator<FT,3U> fiter;
860  // Try to find ngp in local fields and get iterator
861  fiter = getFieldIter(f,ngp);
862 
863  auto W = [dpos](unsigned p, unsigned i) {
864  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
865  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
866  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
867 
868  for (int p0=-1; p0<=1; ++p0) {
869  for (int p1=-1; p1<=1; ++p1) {
870  for (int p2=-1; p2<=1; ++p2) {
871  fiter.offset(p0,p1,p2) += W(p0,0) * W(p1,1) * W(p2,2) * pdata;
872  }
873  }
874  }
875 
876 // std::cout << "HELLOOO from scatter" << std::endl;
877 
878 
879 
880  return;
881  }
882 
883  // gather particle data from Field using particle position and mesh
884  template <class FT, class M, class C, class PT>
885  static
886  void gather(FT& pdata, const Field<FT,3U,M,C>& f,
887  const Vektor<PT,3U>& ppos, const M& mesh) {
888  CenteringTag<C> ctag;
889  Vektor<PT,3U> gpos, dpos, delta;
890  NDIndex<3U> ngp;
891  CompressedBrickIterator<FT,3U> fiter;
892  // find nearest grid point for particle position, store in NDIndex obj
893  ngp = FindNGP(mesh, ppos, ctag);
894  // get position of ngp
895  FindPos(gpos, mesh, ngp, ctag);
896  // get mesh spacings
897  FindDelta(delta, mesh, ngp, ctag);
898  // Try to find ngp in local fields and get iterator
899  fiter = getFieldIter(f,ngp);
900  // get distance from ppos to gpos
901  dpos = ppos - gpos;
902  // normalize dpos by mesh spacing
903  dpos /= delta;
904  // accumulate into particle attrib
905 
906  pdata = 0;
907  auto W = [dpos](unsigned p, unsigned i) {
908  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
909  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
910  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
911 
912  for (int p0=-1; p0<=1; ++p0) {
913  for (int p1=-1; p1<=1; ++p1) {
914  for (int p2=-1; p2<=1; ++p2) {
915  pdata += W(p0,0) * W(p1,1) * W(p2,2) * fiter.offset(p0,p1,p2);
916  }
917  }
918  }
919  //std::cout << "GATHERED pdata = " << pdata << std::endl;
920 
921  return;
922  }
923 
924  // gather particle data from Field using particle position and mesh
925  // and cache mesh information for reuse
926  template <class FT, class M, class C, class PT>
927  static
928  void gather(FT& pdata, const Field<FT,3U,M,C>& f,
929  const Vektor<PT,3U>& ppos, const M& mesh,
930  NDIndex<3U>& ngp, int lgpoff[3U], Vektor<PT,3U>& dpos) {
931  CenteringTag<C> ctag;
932  Vektor<PT,3U> gpos, delta;
933  CompressedBrickIterator<FT,3U> fiter;
934  // find nearest grid point for particle position, store in NDIndex obj
935  ngp = FindNGP(mesh, ppos, ctag);
936  // get position of ngp
937  FindPos(gpos, mesh, ngp, ctag);
938  // get mesh spacings
939  FindDelta(delta, mesh, ngp, ctag);
940  // Try to find ngp in local fields and get iterator
941  fiter = getFieldIter(f,ngp);
942  // get distance from ppos to gpos
943  dpos = ppos - gpos;
944  // normalize dpos by mesh spacing
945  dpos /= delta;
946  // accumulate into particle attrib
947 
948  pdata = 0;
949  auto W = [dpos](unsigned p, unsigned i) {
950  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
951  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
952  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
953 
954  for (int p0=-1; p0<=1; ++p0) {
955  for (int p1=-1; p1<=1; ++p1) {
956  for (int p2=-1; p2<=1; ++p2) {
957  pdata += W(p0,0) * W(p1,1) * W(p2,2) * fiter.offset(p0,p1,p2);
958  }
959  }
960  }
961 
962  return;
963  }
964 
965  // gather particle data from Field using cached mesh information
966  template <class FT, class M, class C, class PT>
967  static
968  void gather(FT& pdata, const Field<FT,3U,M,C>& f,
969  const NDIndex<3U>& ngp, const int lgpoff[3U],
970  const Vektor<PT,3U>& dpos) {
971  CompressedBrickIterator<FT,3U> fiter;
972  // Try to find ngp in local fields and get iterator
973  fiter = getFieldIter(f,ngp);
974  // adjust position of Field iterator to lgp position
975  //fiter.moveBy(lgpoff);
976  // accumulate into particle attrib
977  pdata = 0;
978  auto W = [dpos](unsigned p, unsigned i) {
979  if (p==-1) return .125*(1.-4.*dpos(i)+4.*dpos(i)*dpos(i));
980  else if (p==0) return .25*(3.-4.*dpos(i)*dpos(i));
981  else if (p==+1) return .125*(1.+4.*dpos(i)+4.*dpos(i)*dpos(i)); };
982 
983  for (int p0=-1; p0<=1; ++p0) {
984  for (int p1=-1; p1<=1; ++p1) {
985  for (int p2=-1; p2<=1; ++p2) {
986  pdata += W(p0,0) * W(p1,1) * W(p2,2) * fiter.offset(p0,p1,p2);
987  }
988  }
989  }
990 
991  return;
992  }
993 
994 };
995 
996 
997 // IntTSC class -- what the user sees
998 class IntTSC {
999 
1000 public:
1001  // constructor/destructor
1002  IntTSC() {}
1003  ~IntTSC() {}
1004 
1005  // gather/scatter functions
1006 
1007  // scatter particle data into Field using particle position and mesh
1008  template <class FT, unsigned Dim, class M, class C, class PT>
1009  static
1010  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
1011  const Vektor<PT,Dim>& ppos, const M& mesh) {
1012  IntTSCImpl<Dim>::scatter(pdata,f,ppos,mesh);
1013  }
1014 
1015  // scatter particle data into Field using particle position and mesh
1016  // and cache mesh information for reuse
1017  template <class FT, unsigned Dim, class M, class C, class PT>
1018  static
1019  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
1020  const Vektor<PT,Dim>& ppos, const M& mesh,
1021  CacheDataTSC<PT,Dim>& cache) {
1022  IntTSCImpl<Dim>::scatter(pdata,f,ppos,mesh,cache.Index_m,
1023  cache.Offset_m,cache.Delta_m);
1024  }
1025 
1026  // scatter particle data into Field using cached mesh information
1027  template <class FT, unsigned Dim, class M, class C, class PT>
1028  static
1029  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
1030  const CacheDataTSC<PT,Dim>& cache) {
1031  IntTSCImpl<Dim>::scatter(pdata,f,cache.Index_m,
1032  cache.Offset_m,cache.Delta_m);
1033  }
1034 
1035  // gather particle data from Field using particle position and mesh
1036  template <class FT, unsigned Dim, class M, class C, class PT>
1037  static
1038  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
1039  const Vektor<PT,Dim>& ppos, const M& mesh) {
1040  IntTSCImpl<Dim>::gather(pdata,f,ppos,mesh);
1041  }
1042 
1043  // gather particle data from Field using particle position and mesh
1044  // and cache mesh information for reuse
1045  template <class FT, unsigned Dim, class M, class C, class PT>
1046  static
1047  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
1048  const Vektor<PT,Dim>& ppos, const M& mesh,
1049  CacheDataTSC<PT,Dim>& cache) {
1050  IntTSCImpl<Dim>::gather(pdata,f,ppos,mesh,cache.Index_m,
1051  cache.Offset_m,cache.Delta_m);
1052  }
1053 
1054  // gather particle data from Field using cached mesh information
1055  template <class FT, unsigned Dim, class M, class C, class PT>
1056  static
1057  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
1058  const CacheDataTSC<PT,Dim>& cache) {
1059  IntTSCImpl<Dim>::gather(pdata,f,cache.Index_m,
1060  cache.Offset_m,cache.Delta_m);
1061  }
1062 
1063 };
1064 
1065 #endif // INT_TSC_H
1066 
1067 /***************************************************************************
1068  * $RCSfile: IntTSC.h,v $ $Author: adelmann $
1069  * $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:28 $
1070  * IPPL_VERSION_ID: $Id: IntTSC.h,v 1.1.1.1 2003/01/23 07:40:28 adelmann Exp $
1071  ***************************************************************************/
IntTSCImpl< 1U >::scatter
static void scatter(const FT &pdata, Field< FT, 1U, M, C > &f, const NDIndex< 1U > &ngp, const int lgpoff[1U], const Vektor< PT, 1U > &dpos)
Definition: IntTSC.h:337
IntTSCImpl::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh, NDIndex< Dim > &ngp, int lgpoff[Dim], Vektor< PT, Dim > &dpos)
Definition: IntTSC.h:190
CacheDataTSC::Offset_m
int Offset_m[Dim]
Definition: Interpolator.h:166
IntTSC::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntTSC.h:1010
IntTSCImpl< 3U >::scatter
static void scatter(const FT &pdata, Field< FT, 3U, M, C > &f, const Vektor< PT, 3U > &ppos, const M &mesh, NDIndex< 3U > &ngp, int lgpoff[3U], Vektor< PT, 3U > &dpos)
Definition: IntTSC.h:815
Vektor
Definition: TSVMeta.h:24
T
Definition: rbendmap.h:8
IntTSCImpl< 2U >::gather
static void gather(FT &pdata, const Field< FT, 2U, M, C > &f, const NDIndex< 2U > &ngp, const int lgpoff[2U], const Vektor< PT, 2U > &dpos)
Definition: IntTSC.h:722
IntTSCImpl::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const NDIndex< Dim > &ngp, const int lgpoff[Dim], const Vektor< PT, Dim > &dpos)
Definition: IntTSC.h:230
CacheDataTSC::Delta_m
Vektor< T, Dim > Delta_m
Definition: Interpolator.h:167
IntTSCImpl< 3U >::gather
static void gather(FT &pdata, const Field< FT, 3U, M, C > &f, const Vektor< PT, 3U > &ppos, const M &mesh)
Definition: IntTSC.h:886
BrickIterator::offset
T & offset(int i) const
Definition: BrickIterator.h:112
ERRORMSG
#define ERRORMSG(msg)
Definition: IpplInfo.h:399
IntTSC::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntTSC.h:1038
IntTSCImpl< 3U >::gather
static void gather(FT &pdata, const Field< FT, 3U, M, C > &f, const Vektor< PT, 3U > &ppos, const M &mesh, NDIndex< 3U > &ngp, int lgpoff[3U], Vektor< PT, 3U > &dpos)
Definition: IntTSC.h:928
IntTSCImpl< 2U >::scatter
static void scatter(const FT &pdata, Field< FT, 2U, M, C > &f, const Vektor< PT, 2U > &ppos, const M &mesh, NDIndex< 2U > &ngp, int lgpoff[2U], Vektor< PT, 2U > &dpos)
Definition: IntTSC.h:538
CacheDataTSC::Index_m
NDIndex< Dim > Index_m
Definition: Interpolator.h:165
IntTSCImpl< 3U >::gather
static void gather(FT &pdata, const Field< FT, 3U, M, C > &f, const NDIndex< 3U > &ngp, const int lgpoff[3U], const Vektor< PT, 3U > &dpos)
Definition: IntTSC.h:968
InterpolatorTraits< T, Dim, IntTSC >::Cache_t
CacheDataTSC< T, Dim > Cache_t
Definition: IntTSC.h:30
IntTSCImpl< 1U >::scatter
static void scatter(const FT &pdata, Field< FT, 1U, M, C > &f, const Vektor< PT, 1U > &ppos, const M &mesh, NDIndex< 1U > &ngp, int lgpoff[1U], Vektor< PT, 1U > &dpos)
Definition: IntTSC.h:298
IntTSCImpl::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const NDIndex< Dim > &ngp, const int lgpoff[Dim], const Vektor< PT, Dim > &dpos)
Definition: IntTSC.h:132
IntTSCImpl::~IntTSCImpl
~IntTSCImpl()
Definition: IntTSC.h:42
IntTSCImpl::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh, NDIndex< Dim > &ngp, int lgpoff[Dim], Vektor< PT, Dim > &dpos)
Definition: IntTSC.h:92
IntTSC::~IntTSC
~IntTSC()
Definition: IntTSC.h:1003
IntTSC::IntTSC
IntTSC()
Definition: IntTSC.h:1002
IntTSCImpl::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntTSC.h:148
FindDelta
void FindDelta(Vektor< PT, Dim > &delta, const M &mesh, const NDIndex< Dim > &gp, CenteringTag< Cell >)
Definition: Interpolator.h:103
IntTSCImpl< 2U >::scatter
static void scatter(const FT &pdata, Field< FT, 2U, M, C > &f, const Vektor< PT, 2U > &ppos, const M &mesh)
Definition: IntTSC.h:491
FindPos
void FindPos(Vektor< PT, Dim > &pos, const M &mesh, const NDIndex< Dim > &indices, CenteringTag< Cell >)
Definition: Interpolator.h:70
IntTSCImpl< 2U >::scatter
static void scatter(const FT &pdata, Field< FT, 2U, M, C > &f, const NDIndex< 2U > &ngp, const int lgpoff[2U], const Vektor< PT, 2U > &dpos)
Definition: IntTSC.h:586
IntTSC
Definition: IntTSC.h:998
Interpolator::getFieldIter
static CompressedBrickIterator< T, Dim > getFieldIter(const BareField< T, Dim > &f, const NDIndex< Dim > &pt)
Definition: Interpolator.h:191
IntTSCImpl< 1U >::gather
static void gather(FT &pdata, const Field< FT, 1U, M, C > &f, const Vektor< PT, 1U > &ppos, const M &mesh)
Definition: IntTSC.h:366
FindNGP
NDIndex< Dim > FindNGP(const M &mesh, const Vektor< PT, Dim > &ppos, CenteringTag< Cell >)
Definition: Interpolator.h:45
IntTSC::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const CacheDataTSC< PT, Dim > &cache)
Definition: IntTSC.h:1057
IntTSC::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const CacheDataTSC< PT, Dim > &cache)
Definition: IntTSC.h:1029
IntTSC::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh, CacheDataTSC< PT, Dim > &cache)
Definition: IntTSC.h:1019
IntTSC::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh, CacheDataTSC< PT, Dim > &cache)
Definition: IntTSC.h:1047
IntTSCImpl< 1U >::gather
static void gather(FT &pdata, const Field< FT, 1U, M, C > &f, const Vektor< PT, 1U > &ppos, const M &mesh, NDIndex< 1U > &ngp, int lgpoff[1U], Vektor< PT, 1U > &dpos)
Definition: IntTSC.h:409
IntTSCImpl::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntTSC.h:49
IntTSCImpl< 1U >::gather
static void gather(FT &pdata, const Field< FT, 1U, M, C > &f, const NDIndex< 1U > &ngp, const int lgpoff[1U], const Vektor< PT, 1U > &dpos)
Definition: IntTSC.h:450
IntTSCImpl::IntTSCImpl
IntTSCImpl()
Definition: IntTSC.h:41
Dim
const unsigned Dim
Definition: P3MPoissonSolver.h:23
IntTSCImpl< 2U >::gather
static void gather(FT &pdata, const Field< FT, 2U, M, C > &f, const Vektor< PT, 2U > &ppos, const M &mesh)
Definition: IntTSC.h:624
IntTSCImpl< 3U >::scatter
static void scatter(const FT &pdata, Field< FT, 3U, M, C > &f, const Vektor< PT, 3U > &ppos, const M &mesh)
Definition: IntTSC.h:772
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
IntTSCImpl< 2U >::gather
static void gather(FT &pdata, const Field< FT, 2U, M, C > &f, const Vektor< PT, 2U > &ppos, const M &mesh, NDIndex< 2U > &ngp, int lgpoff[2U], Vektor< PT, 2U > &dpos)
Definition: IntTSC.h:673
IntTSCImpl< 3U >::scatter
static void scatter(const FT &pdata, Field< FT, 3U, M, C > &f, const NDIndex< 3U > &ngp, const int lgpoff[3U], const Vektor< PT, 3U > &dpos)
Definition: IntTSC.h:856
IntTSCImpl< 1U >::scatter
static void scatter(const FT &pdata, Field< FT, 1U, M, C > &f, const Vektor< PT, 1U > &ppos, const M &mesh)
Definition: IntTSC.h:259
Generated on Sat Jan 4 2020 14:47:37 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5