OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
IntCIC.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_CIC_H
12 #define INT_CIC_H
13 
14 /* IntCIC.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 IntCIC;
26 
27 // specialization of InterpolatorTraits
28 template <class T, unsigned Dim>
29 struct InterpolatorTraits<T,Dim,IntCIC> {
30  typedef CacheDataCIC<T,Dim> Cache_t;
31 };
32 
33 
34 
35 // IntCICImpl class definition
36 template <unsigned Dim>
37 class IntCICImpl : public Interpolator {
38 
39 public:
40  // constructor/destructor
41  IntCICImpl() {}
42  ~IntCICImpl() {}
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  int lgpoff[Dim];
56  unsigned int d;
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  for (d=0; d<Dim; ++d) {
67  if (gpos(d) > ppos(d)) {
68  lgpoff[d] = -1; // save the offset
69  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
70  }
71  else {
72  lgpoff[d] = 0; // save the offset
73  }
74  }
75  // adjust position of Field iterator to lgp position
76  fiter.moveBy(lgpoff);
77  // get distance from ppos to gpos
78  dpos = ppos - gpos;
79  // normalize dpos by mesh spacing
80  dpos /= delta;
81  // accumulate into local elements
82  ERRORMSG("IntCIC::scatter: not implemented for Dim>3!!"<<endl);
83  return;
84  }
85 
86  // scatter particle data into Field using particle position and mesh
87  // and cache mesh information for reuse
88  template <class FT, class M, class C, class PT>
89  static
90  void scatter(const FT& /*pdata*/, Field<FT,Dim,M,C>& f,
91  const Vektor<PT,Dim>& ppos, const M& mesh,
92  NDIndex<Dim>& ngp, int lgpoff[Dim], Vektor<PT,Dim>& dpos) {
93  CenteringTag<C> ctag;
94  Vektor<PT,Dim> gpos, delta;
95  CompressedBrickIterator<FT,Dim> fiter;
96  unsigned int d;
97  // find nearest grid point for particle position, store in NDIndex obj
98  ngp = FindNGP(mesh, ppos, ctag);
99  // get position of ngp
100  FindPos(gpos, mesh, ngp, ctag);
101  // get mesh spacings
102  FindDelta(delta, mesh, ngp, ctag);
103  // Try to find ngp in local fields and get iterator
104  fiter = getFieldIter(f,ngp);
105  // Now we find the offset from ngp to next-lowest grip point (lgp)
106  for (d=0; d<Dim; ++d) {
107  if (gpos(d) > ppos(d)) {
108  lgpoff[d] = -1; // save the offset
109  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
110  }
111  else {
112  lgpoff[d] = 0; // save the offset
113  }
114  }
115  // adjust position of Field iterator to lgp position
116  fiter.moveBy(lgpoff);
117  // get distance from ppos to gpos
118  dpos = ppos - gpos;
119  // normalize dpos by mesh spacing
120  dpos /= delta;
121  // accumulate into local elements
122  ERRORMSG("IntCIC::scatter: not implemented for Dim>3!!"<<endl);
123  return;
124  }
125 
126  // scatter particle data into Field using cached mesh information
127  template <class FT, class M, class C, class PT>
128  static
129  void scatter(const FT& /*pdata*/, Field<FT,Dim,M,C>& f,
130  const NDIndex<Dim>& ngp, const int lgpoff[Dim],
131  const Vektor<PT,Dim>& /*dpos*/) {
132  CompressedBrickIterator<FT,Dim> fiter;
133  // Try to find ngp in local fields and get iterator
134  fiter = getFieldIter(f,ngp);
135  // adjust position of Field iterator to lgp position
136  fiter.moveBy(lgpoff);
137  // accumulate into local elements
138  ERRORMSG("IntCIC::scatter: not implemented for Dim>3!!"<<endl);
139  return;
140  }
141 
142  // gather particle data from Field using particle position and mesh
143  template <class FT, class M, class C, class PT>
144  static
145  void gather(FT& /*pdata*/, const Field<FT,Dim,M,C>& f,
146  const Vektor<PT,Dim>& ppos, const M& mesh) {
147  CenteringTag<C> ctag;
148  Vektor<PT,Dim> gpos, dpos, delta;
149  NDIndex<Dim> ngp;
150  CompressedBrickIterator<FT,Dim> fiter;
151  int lgpoff[Dim];
152  unsigned int d;
153  // find nearest grid point for particle position, store in NDIndex obj
154  ngp = FindNGP(mesh, ppos, ctag);
155  // get position of ngp
156  FindPos(gpos, mesh, ngp, ctag);
157  // get mesh spacings
158  FindDelta(delta, mesh, ngp, ctag);
159  // Try to find ngp in local fields and get iterator
160  fiter = getFieldIter(f,ngp);
161  // Now we find the offset from ngp to next-lowest grip point (lgp)
162  for (d=0; d<Dim; ++d) {
163  if (gpos(d) > ppos(d)) {
164  lgpoff[d] = -1; // save the offset
165  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
166  }
167  else {
168  lgpoff[d] = 0; // save the offset
169  }
170  }
171  // adjust position of Field iterator to lgp position
172  fiter.moveBy(lgpoff);
173  // get distance from ppos to gpos
174  dpos = ppos - gpos;
175  // normalize dpos by mesh spacing
176  dpos /= delta;
177  // accumulate into particle attrib
178  ERRORMSG("IntCIC::gather: not implemented for Dim>3!!"<<endl);
179  return;
180  }
181 
182  // gather particle data from Field using particle position and mesh
183  // and cache mesh information for reuse
184  template <class FT, class M, class C, class PT>
185  static
186  void gather(FT& /*pdata*/, const Field<FT,Dim,M,C>& f,
187  const Vektor<PT,Dim>& ppos, const M& mesh,
188  NDIndex<Dim>& ngp, int lgpoff[Dim], Vektor<PT,Dim>& dpos) {
189  CenteringTag<C> ctag;
190  Vektor<PT,Dim> gpos, delta;
191  CompressedBrickIterator<FT,Dim> fiter;
192  unsigned int d;
193  // find nearest grid point for particle position, store in NDIndex obj
194  ngp = FindNGP(mesh, ppos, ctag);
195  // get position of ngp
196  FindPos(gpos, mesh, ngp, ctag);
197  // get mesh spacings
198  FindDelta(delta, mesh, ngp, ctag);
199  // Try to find ngp in local fields and get iterator
200  fiter = getFieldIter(f,ngp);
201  // Now we find the offset from ngp to next-lowest grip point (lgp)
202  for (d=0; d<Dim; ++d) {
203  if (gpos(d) > ppos(d)) {
204  lgpoff[d] = -1; // save the offset
205  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
206  }
207  else {
208  lgpoff[d] = 0; // save the offset
209  }
210  }
211  // adjust position of Field iterator to lgp position
212  fiter.moveBy(lgpoff);
213  // get distance from ppos to gpos
214  dpos = ppos - gpos;
215  // normalize dpos by mesh spacing
216  dpos /= delta;
217  // accumulate into particle attrib
218  ERRORMSG("IntCIC::gather: not implemented for Dim>3!!"<<endl);
219  return;
220  }
221 
222  // gather particle data from Field using cached mesh information
223  template <class FT, class M, class C, class PT>
224  static
225  void gather(FT& /*pdata*/, const Field<FT,Dim,M,C>& f,
226  const NDIndex<Dim>& ngp, const int lgpoff[Dim],
227  const Vektor<PT,Dim>&) {
228  CompressedBrickIterator<FT,Dim> fiter;
229  // Try to find ngp in local fields and get iterator
230  fiter = getFieldIter(f,ngp);
231  // adjust position of Field iterator to lgp position
232  fiter.moveBy(lgpoff);
233  // accumulate into particle attrib
234  ERRORMSG("IntCIC::gather: not implemented for Dim>3!!"<<endl);
235  return;
236  }
237 
238 };
239 
240 
241 template <>
242 class IntCICImpl<1U> : public Interpolator {
243 
244 public:
245  // constructor/destructor
246  IntCICImpl() {}
247  ~IntCICImpl() {}
248 
249  // gather/scatter functions
250 
251  // scatter particle data into Field using particle position and mesh
252  template <class FT, class M, class C, class PT>
253  static
254  void scatter(const FT& pdata, Field<FT,1U,M,C>& f,
255  const Vektor<PT,1U>& ppos, const M& mesh) {
256  CenteringTag<C> ctag;
257  Vektor<PT,1U> gpos, dpos, delta;
258  NDIndex<1U> ngp;
259  CompressedBrickIterator<FT,1U> fiter;
260  int lgpoff[1U];
261  unsigned int d;
262  // find nearest grid point for particle position, store in NDIndex obj
263  ngp = FindNGP(mesh, ppos, ctag);
264  // get position of ngp
265  FindPos(gpos, mesh, ngp, ctag);
266  // get mesh spacings
267  FindDelta(delta, mesh, ngp, ctag);
268  // Try to find ngp in local fields and get iterator
269  fiter = getFieldIter(f,ngp);
270  // Now we find the offset from ngp to next-lowest grip point (lgp)
271  for (d=0; d<1U; ++d) {
272  if (gpos(d) > ppos(d)) {
273  lgpoff[d] = -1; // save the offset
274  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
275  }
276  else {
277  lgpoff[d] = 0; // save the offset
278  }
279  }
280  // adjust position of Field iterator to lgp position
281  fiter.moveBy(lgpoff);
282  // get distance from ppos to gpos
283  dpos = ppos - gpos;
284  // normalize dpos by mesh spacing
285  dpos /= delta;
286  // accumulate into local elements
287  *fiter += (1 - dpos(0)) * pdata;
288  fiter.offset(1) += dpos(0) * pdata;
289  return;
290  }
291 
292  // scatter particle data into Field using particle position and mesh
293  // and cache mesh information for reuse
294  template <class FT, class M, class C, class PT>
295  static
296  void scatter(const FT& pdata, Field<FT,1U,M,C>& f,
297  const Vektor<PT,1U>& ppos, const M& mesh,
298  NDIndex<1U>& ngp, int lgpoff[1U], Vektor<PT,1U>& dpos) {
299  CenteringTag<C> ctag;
300  Vektor<PT,1U> gpos, delta;
301  CompressedBrickIterator<FT,1U> fiter;
302  unsigned int d;
303  // find nearest grid point for particle position, store in NDIndex obj
304  ngp = FindNGP(mesh, ppos, ctag);
305  // get position of ngp
306  FindPos(gpos, mesh, ngp, ctag);
307  // get mesh spacings
308  FindDelta(delta, mesh, ngp, ctag);
309  // Try to find ngp in local fields and get iterator
310  fiter = getFieldIter(f,ngp);
311  // Now we find the offset from ngp to next-lowest grip point (lgp)
312  for (d=0; d<1U; ++d) {
313  if (gpos(d) > ppos(d)) {
314  lgpoff[d] = -1; // save the offset
315  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
316  }
317  else {
318  lgpoff[d] = 0; // save the offset
319  }
320  }
321  // adjust position of Field iterator to lgp position
322  fiter.moveBy(lgpoff);
323  // get distance from ppos to gpos
324  dpos = ppos - gpos;
325  // normalize dpos by mesh spacing
326  dpos /= delta;
327  // accumulate into local elements
328  *fiter += (1 - dpos(0)) * pdata;
329  fiter.offset(1) += dpos(0) * pdata;
330  return;
331  }
332 
333  // scatter particle data into Field using cached mesh information
334  template <class FT, class M, class C, class PT>
335  static
336  void scatter(const FT& pdata, Field<FT,1U,M,C>& f,
337  const NDIndex<1U>& ngp, const int lgpoff[1U],
338  const Vektor<PT,1U>& dpos) {
339  CompressedBrickIterator<FT,1U> fiter;
340  // Try to find ngp in local fields and get iterator
341  fiter = getFieldIter(f,ngp);
342  // adjust position of Field iterator to lgp position
343  fiter.moveBy(lgpoff);
344  // accumulate into local elements
345  *fiter += (1 - dpos(0)) * pdata;
346  fiter.offset(1) += dpos(0) * pdata;
347  return;
348  }
349 
350  // gather particle data from Field using particle position and mesh
351  template <class FT, class M, class C, class PT>
352  static
353  void gather(FT& pdata, const Field<FT,1U,M,C>& f,
354  const Vektor<PT,1U>& ppos, const M& mesh) {
355  CenteringTag<C> ctag;
356  Vektor<PT,1U> gpos, dpos, delta;
357  NDIndex<1U> ngp;
358  CompressedBrickIterator<FT,1U> fiter;
359  int lgpoff[1U];
360  unsigned int d;
361  // find nearest grid point for particle position, store in NDIndex obj
362  ngp = FindNGP(mesh, ppos, ctag);
363  // get position of ngp
364  FindPos(gpos, mesh, ngp, ctag);
365  // get mesh spacings
366  FindDelta(delta, mesh, ngp, ctag);
367  // Try to find ngp in local fields and get iterator
368  fiter = getFieldIter(f,ngp);
369  // Now we find the offset from ngp to next-lowest grip point (lgp)
370  for (d=0; d<1U; ++d) {
371  if (gpos(d) > ppos(d)) {
372  lgpoff[d] = -1; // save the offset
373  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
374  }
375  else {
376  lgpoff[d] = 0; // save the offset
377  }
378  }
379  // adjust position of Field iterator to lgp position
380  fiter.moveBy(lgpoff);
381  // get distance from ppos to gpos
382  dpos = ppos - gpos;
383  // normalize dpos by mesh spacing
384  dpos /= delta;
385  // accumulate into particle attrib
386  pdata = (1 - dpos(0)) * (*fiter) +
387  dpos(0) * fiter.offset(1);
388  return;
389  }
390 
391  // gather particle data from Field using particle position and mesh
392  // and cache mesh information for reuse
393  template <class FT, class M, class C, class PT>
394  static
395  void gather(FT& pdata, const Field<FT,1U,M,C>& f,
396  const Vektor<PT,1U>& ppos, const M& mesh,
397  NDIndex<1U>& ngp, int lgpoff[1U], Vektor<PT,1U>& dpos) {
398  CenteringTag<C> ctag;
399  Vektor<PT,1U> gpos, delta;
400  CompressedBrickIterator<FT,1U> fiter;
401  unsigned int d;
402  // find nearest grid point for particle position, store in NDIndex obj
403  ngp = FindNGP(mesh, ppos, ctag);
404  // get position of ngp
405  FindPos(gpos, mesh, ngp, ctag);
406  // get mesh spacings
407  FindDelta(delta, mesh, ngp, ctag);
408  // Try to find ngp in local fields and get iterator
409  fiter = getFieldIter(f,ngp);
410  // Now we find the offset from ngp to next-lowest grip point (lgp)
411  for (d=0; d<1U; ++d) {
412  if (gpos(d) > ppos(d)) {
413  lgpoff[d] = -1; // save the offset
414  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
415  }
416  else {
417  lgpoff[d] = 0; // save the offset
418  }
419  }
420  // adjust position of Field iterator to lgp position
421  fiter.moveBy(lgpoff);
422  // get distance from ppos to gpos
423  dpos = ppos - gpos;
424  // normalize dpos by mesh spacing
425  dpos /= delta;
426  // accumulate into particle attrib
427  pdata = (1 - dpos(0)) * (*fiter) +
428  dpos(0) * fiter.offset(1);
429  return;
430  }
431 
432  // gather particle data from Field using cached mesh information
433  template <class FT, class M, class C, class PT>
434  static
435  void gather(FT& pdata, const Field<FT,1U,M,C>& f,
436  const NDIndex<1U>& ngp, const int lgpoff[1U],
437  const Vektor<PT,1U>& dpos) {
438  CompressedBrickIterator<FT,1U> fiter;
439  // Try to find ngp in local fields and get iterator
440  fiter = getFieldIter(f,ngp);
441  // adjust position of Field iterator to lgp position
442  fiter.moveBy(lgpoff);
443  // accumulate into particle attrib
444  pdata = (1 - dpos(0)) * (*fiter) +
445  dpos(0) * fiter.offset(1);
446  return;
447  }
448 
449 };
450 
451 
452 template <>
453 class IntCICImpl<2U> : public Interpolator {
454 
455 public:
456  // constructor/destructor
457  IntCICImpl() {}
458  ~IntCICImpl() {}
459 
460  // gather/scatter functions
461 
462  // scatter particle data into Field using particle position and mesh
463  template <class FT, class M, class C, class PT>
464  static
465  void scatter(const FT& pdata, Field<FT,2U,M,C>& f,
466  const Vektor<PT,2U>& ppos, const M& mesh) {
467  CenteringTag<C> ctag;
468  Vektor<PT,2U> gpos, dpos, delta;
469  NDIndex<2U> ngp;
470  CompressedBrickIterator<FT,2U> fiter;
471  int lgpoff[2U];
472  unsigned int d;
473  // find nearest grid point for particle position, store in NDIndex obj
474  ngp = FindNGP(mesh, ppos, ctag);
475  // get position of ngp
476  FindPos(gpos, mesh, ngp, ctag);
477  // get mesh spacings
478  FindDelta(delta, mesh, ngp, ctag);
479  // Try to find ngp in local fields and get iterator
480  fiter = getFieldIter(f,ngp);
481  // Now we find the offset from ngp to next-lowest grip point (lgp)
482  for (d=0; d<2U; ++d) {
483  if (gpos(d) > ppos(d)) {
484  lgpoff[d] = -1; // save the offset
485  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
486  }
487  else {
488  lgpoff[d] = 0; // save the offset
489  }
490  }
491  // adjust position of Field iterator to lgp position
492  fiter.moveBy(lgpoff);
493  // get distance from ppos to gpos
494  dpos = ppos - gpos;
495  // normalize dpos by mesh spacing
496  dpos /= delta;
497  // accumulate into local elements
498  *fiter += (1 - dpos(0)) * (1 - dpos(1)) * pdata;
499  fiter.offset(1,0) += dpos(0) * (1 - dpos(1)) * pdata;
500  fiter.offset(0,1) += (1 - dpos(0)) * dpos(1) * pdata;
501  fiter.offset(1,1) += dpos(0) * dpos(1) * pdata;
502  return;
503  }
504 
505  // scatter particle data into Field using particle position and mesh
506  // and cache mesh information for reuse
507  template <class FT, class M, class C, class PT>
508  static
509  void scatter(const FT& pdata, Field<FT,2U,M,C>& f,
510  const Vektor<PT,2U>& ppos, const M& mesh,
511  NDIndex<2U>& ngp, int lgpoff[2U], Vektor<PT,2U>& dpos) {
512  CenteringTag<C> ctag;
513  Vektor<PT,2U> gpos, delta;
514  CompressedBrickIterator<FT,2U> fiter;
515  unsigned int d;
516  // find nearest grid point for particle position, store in NDIndex obj
517  ngp = FindNGP(mesh, ppos, ctag);
518  // get position of ngp
519  FindPos(gpos, mesh, ngp, ctag);
520  // get mesh spacings
521  FindDelta(delta, mesh, ngp, ctag);
522  // Try to find ngp in local fields and get iterator
523  fiter = getFieldIter(f,ngp);
524  // Now we find the offset from ngp to next-lowest grip point (lgp)
525  for (d=0; d<2U; ++d) {
526  if (gpos(d) > ppos(d)) {
527  lgpoff[d] = -1; // save the offset
528  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
529  }
530  else {
531  lgpoff[d] = 0; // save the offset
532  }
533  }
534  // adjust position of Field iterator to lgp position
535  fiter.moveBy(lgpoff);
536  // get distance from ppos to gpos
537  dpos = ppos - gpos;
538  // normalize dpos by mesh spacing
539  dpos /= delta;
540  // accumulate into local elements
541  *fiter += (1 - dpos(0)) * (1 - dpos(1)) * pdata;
542  fiter.offset(1,0) += dpos(0) * (1 - dpos(1)) * pdata;
543  fiter.offset(0,1) += (1 - dpos(0)) * dpos(1) * pdata;
544  fiter.offset(1,1) += dpos(0) * dpos(1) * pdata;
545  return;
546  }
547 
548  // scatter particle data into Field using cached mesh information
549  template <class FT, class M, class C, class PT>
550  static
551  void scatter(const FT& pdata, Field<FT,2U,M,C>& f,
552  const NDIndex<2U>& ngp, const int lgpoff[2U],
553  const Vektor<PT,2U>& dpos) {
554  CompressedBrickIterator<FT,2U> fiter;
555  // Try to find ngp in local fields and get iterator
556  fiter = getFieldIter(f,ngp);
557  // adjust position of Field iterator to lgp position
558  fiter.moveBy(lgpoff);
559  // accumulate into local elements
560  *fiter += (1 - dpos(0)) * (1 - dpos(1)) * pdata;
561  fiter.offset(1,0) += dpos(0) * (1 - dpos(1)) * pdata;
562  fiter.offset(0,1) += (1 - dpos(0)) * dpos(1) * pdata;
563  fiter.offset(1,1) += dpos(0) * dpos(1) * pdata;
564  return;
565  }
566 
567  // gather particle data from Field using particle position and mesh
568  template <class FT, class M, class C, class PT>
569  static
570  void gather(FT& pdata, const Field<FT,2U,M,C>& f,
571  const Vektor<PT,2U>& ppos, const M& mesh) {
572  CenteringTag<C> ctag;
573  Vektor<PT,2U> gpos, dpos, delta;
574  NDIndex<2U> ngp;
575  CompressedBrickIterator<FT,2U> fiter;
576  int lgpoff[2U];
577  unsigned int d;
578  // find nearest grid point for particle position, store in NDIndex obj
579  ngp = FindNGP(mesh, ppos, ctag);
580  // get position of ngp
581  FindPos(gpos, mesh, ngp, ctag);
582  // get mesh spacings
583  FindDelta(delta, mesh, ngp, ctag);
584  // Try to find ngp in local fields and get iterator
585  fiter = getFieldIter(f,ngp);
586  // Now we find the offset from ngp to next-lowest grip point (lgp)
587  for (d=0; d<2U; ++d) {
588  if (gpos(d) > ppos(d)) {
589  lgpoff[d] = -1; // save the offset
590  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
591  }
592  else {
593  lgpoff[d] = 0; // save the offset
594  }
595  }
596  // adjust position of Field iterator to lgp position
597  fiter.moveBy(lgpoff);
598  // get distance from ppos to gpos
599  dpos = ppos - gpos;
600  // normalize dpos by mesh spacing
601  dpos /= delta;
602  // accumulate into particle attrib
603  pdata = (1 - dpos(0)) * (1 - dpos(1)) * (*fiter) +
604  dpos(0) * (1 - dpos(1)) * fiter.offset(1,0) +
605  (1 - dpos(0)) * dpos(1) * fiter.offset(0,1) +
606  dpos(0) * dpos(1) * fiter.offset(1,1);
607  return;
608  }
609 
610  // gather particle data from Field using particle position and mesh
611  // and cache mesh information for reuse
612  template <class FT, class M, class C, class PT>
613  static
614  void gather(FT& pdata, const Field<FT,2U,M,C>& f,
615  const Vektor<PT,2U>& ppos, const M& mesh,
616  NDIndex<2U>& ngp, int lgpoff[2U], Vektor<PT,2U>& dpos) {
617  CenteringTag<C> ctag;
618  Vektor<PT,2U> gpos, delta;
619  CompressedBrickIterator<FT,2U> fiter;
620  unsigned int d;
621  // find nearest grid point for particle position, store in NDIndex obj
622  ngp = FindNGP(mesh, ppos, ctag);
623  // get position of ngp
624  FindPos(gpos, mesh, ngp, ctag);
625  // get mesh spacings
626  FindDelta(delta, mesh, ngp, ctag);
627  // Try to find ngp in local fields and get iterator
628  fiter = getFieldIter(f,ngp);
629  // Now we find the offset from ngp to next-lowest grip point (lgp)
630  for (d=0; d<2U; ++d) {
631  if (gpos(d) > ppos(d)) {
632  lgpoff[d] = -1; // save the offset
633  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
634  }
635  else {
636  lgpoff[d] = 0; // save the offset
637  }
638  }
639  // adjust position of Field iterator to lgp position
640  fiter.moveBy(lgpoff);
641  // get distance from ppos to gpos
642  dpos = ppos - gpos;
643  // normalize dpos by mesh spacing
644  dpos /= delta;
645  // accumulate into particle attrib
646  pdata = (1 - dpos(0)) * (1 - dpos(1)) * (*fiter) +
647  dpos(0) * (1 - dpos(1)) * fiter.offset(1,0) +
648  (1 - dpos(0)) * dpos(1) * fiter.offset(0,1) +
649  dpos(0) * dpos(1) * fiter.offset(1,1);
650  return;
651  }
652 
653  // gather particle data from Field using cached mesh information
654  template <class FT, class M, class C, class PT>
655  static
656  void gather(FT& pdata, const Field<FT,2U,M,C>& f,
657  const NDIndex<2U>& ngp, const int lgpoff[2U],
658  const Vektor<PT,2U>& dpos) {
659  CompressedBrickIterator<FT,2U> fiter;
660  // Try to find ngp in local fields and get iterator
661  fiter = getFieldIter(f,ngp);
662  // adjust position of Field iterator to lgp position
663  fiter.moveBy(lgpoff);
664  // accumulate into particle attrib
665  pdata = (1 - dpos(0)) * (1 - dpos(1)) * (*fiter) +
666  dpos(0) * (1 - dpos(1)) * fiter.offset(1,0) +
667  (1 - dpos(0)) * dpos(1) * fiter.offset(0,1) +
668  dpos(0) * dpos(1) * fiter.offset(1,1);
669  return;
670  }
671 
672  // scatter particle data into Field using particle position and mesh
673  template <class FT, class M, class PT>
674  static
675  void scatter(const Vektor<FT,2U>& pdata, Field<Vektor<FT,2U>,2U,M,Edge>& f,
676  const Vektor<PT,2U>& ppos, const M& mesh) {
677  CenteringTag<Cell> ctag;
678  Vektor<PT,2U> ppos2, dpos, gpos, delta;
679  NDIndex<2U> ngp;
680 
681  CompressedBrickIterator<Vektor<FT,2U>,2U> fiter;
682 
683  // find nearest grid point for particle position, store in NDIndex obj
684  ngp = FindNGP(mesh, ppos, ctag);
685  FindDelta(delta, mesh, ngp, ctag);
686 
687  for (unsigned int comp = 0; comp < 2U; ++comp) {
688  Vektor<PT,2U> cppos = ppos;
689  cppos(comp) += 0.5 * delta(comp);
690  ngp = FindNGP(mesh, cppos, ctag);
691  ngp[comp] = ngp[comp] - 1;
692 
693  FindPos(gpos, mesh, ngp, CenteringTag<Vert>());
694  gpos(comp) += 0.5 * delta(comp);
695  dpos = ppos - gpos;
696  dpos /= delta;
697 
698  // Try to find ngp in local fields and get iterator
699  fiter = Interpolator::getFieldIter(f,ngp);
700 
701  // accumulate into local elements
702  (*fiter)(comp) += (1 - dpos(0)) * (1 - dpos(1)) * pdata(comp);
703  (fiter.offset(1,0))(comp) += dpos(0) * (1 - dpos(1)) * pdata(comp);
704  (fiter.offset(0,1))(comp) += (1 - dpos(0)) * dpos(1) * pdata(comp);
705  (fiter.offset(1,1))(comp) += dpos(0) * dpos(1) * pdata(comp);
706  }
707  return;
708  }
709 
710  template <class FT, class M, class PT>
711  static
712  void scatter(const FT& /*pdata*/, Field<FT,2U,M,Edge>& /*f*/,
713  const Vektor<PT,2U>& /*ppos*/, const M& /*mesh*/) {
714  ERRORMSG("IntCIC::scatter on Edge based field: not implemented for non-vectors!!"<<endl);
715  return;
716  }
717  // scatter particle data into Field using particle position and mesh
718  // and cache mesh information for reuse
719  template <class FT, class M, class PT>
720  static
721  void scatter(const FT& /*pdata*/, Field<FT,2U,M,Edge>& /*f*/,
722  const Vektor<PT,2U>& /*ppos*/, const M& /*mesh*/,
723  NDIndex<2U>& /*ngp*/, int /*lgpoff*/ [2U], Vektor<PT,2U>& /*dpos*/) {
724  ERRORMSG("IntCIC::scatter on Edge based field with cache: not implemented!!"<<endl);
725  return;
726  }
727 
728  template <class FT, class M, class PT>
729  static
730  void scatter(const FT&, Field<FT,2U,M,Edge>&,
731  const NDIndex<2U>&, const int [2U],
732  const Vektor<PT,2U>&) {
733  ERRORMSG("IntCIC::scatter on Edge based field with cache: not implemented!!"<<endl);
734  return;
735  }
736 
737  // gather particle data from Field using particle position and mesh
738  template <class FT, class M, class PT>
739  static
740  void gather(Vektor<FT,2U>& pdata, const Field<Vektor<FT,2U>,2U,M,Edge>& f,
741  const Vektor<PT,2U>& ppos, const M& mesh) {
742  Vektor<PT,2U> dpos, gpos, delta;
743  NDIndex<2U> ngp = mesh.getNearestVertex(ppos);;
744  CompressedBrickIterator<Vektor<FT,2U>,2U> fiter;
745 
746  FindDelta(delta, mesh, ngp, CenteringTag<Vert>());
747 
748  for (unsigned int comp = 0; comp < 2U; ++ comp) {
749  Vektor<PT,2U> cppos = ppos;
750  cppos(comp) += 0.5 * delta(comp);
751  ngp = mesh.getCellContaining(cppos);
752  ngp[comp] = ngp[comp] - 1;
753  gpos = mesh.getVertexPosition(ngp);
754  gpos(comp) += 0.5 * delta(comp);
755 
756  // get distance from ppos to gpos
757  dpos = ppos - gpos;
758  // normalize dpos by mesh spacing
759  dpos /= delta;
760  // accumulate into local elements
761 
762  // Try to find ngp in local fields and get iterator
763  fiter = Interpolator::getFieldIter(f,ngp);
764 
765  // accumulate into particle attrib
766  pdata(comp) = ((1 - dpos(0)) * (1 - dpos(1)) * (*fiter)(comp) +
767  dpos(0) * (1 - dpos(1)) * (fiter.offset(1,0))(comp) +
768  (1 - dpos(0)) * dpos(1) * (fiter.offset(0,1))(comp) +
769  dpos(0) * dpos(1) * (fiter.offset(1,1))(comp));
770  }
771  return;
772  }
773 
774  template <class FT, class M, class PT>
775  static
776  void gather(FT& /*pdata*/, const Field<FT,2U,M,Edge>& /*f*/,
777  const Vektor<PT,2U>& /*ppos*/, const M& /*mesh*/) {
778  ERRORMSG("IntCIC::gather on Edge based field: not implemented for non-vectors!!"<<endl);
779  return;
780  }
781 
782  // gather particle data from Field using particle position and mesh
783  // and cache mesh information for reuse
784  template <class FT, class M, class PT>
785  static
786  void gather(FT& /*pdata*/, const Field<FT,2U,M,Edge>& /*f*/,
787  const Vektor<PT,2U>& /*ppos*/, const M& /*mesh*/,
788  NDIndex<2U>& /*ngp*/, int /*lgpoff*/[2U], Vektor<PT,2U>& /*dpos*/) {
789  ERRORMSG("IntCIC::gather on Edge based field with cache: not implemented!!"<<endl);
790  return;
791  }
792 
793  // gather particle data from Field using cached mesh information
794  template <class FT, class M, class PT>
795  static
796  void gather(FT& /*pdata*/, const Field<FT,2U,M,Edge>& /*f*/,
797  const NDIndex<2U>& /*ngp*/, const int /*lgpoff*/[2U],
798  const Vektor<PT,2U>& /*dpos*/) {
799  ERRORMSG("IntCIC::gather on Edge based field with cache: not implemented!!"<<endl);
800  return;
801  }
802 };
803 
804 
805 template <>
806 class IntCICImpl<3U> : public Interpolator {
807 
808 public:
809  // constructor/destructor
810  IntCICImpl() {}
811  ~IntCICImpl() {}
812 
813  // gather/scatter functions
814 
815  // scatter particle data into Field using particle position and mesh
816  template <class FT, class M, class C, class PT>
817  static
818  void scatter(const FT& pdata, Field<FT,3U,M,C>& f,
819  const Vektor<PT,3U>& ppos, const M& mesh) {
820  CenteringTag<C> ctag;
821  Vektor<PT,3U> gpos, dpos, delta;
822  NDIndex<3U> ngp;
823  CompressedBrickIterator<FT,3U> fiter;
824  int lgpoff[3U];
825  unsigned int d;
826  // find nearest grid point for particle position, store in NDIndex obj
827  ngp = FindNGP(mesh, ppos, ctag);
828  // get position of ngp
829  FindPos(gpos, mesh, ngp, ctag);
830  // get mesh spacings
831  FindDelta(delta, mesh, ngp, ctag);
832  // Try to find ngp in local fields and get iterator
833  fiter = getFieldIter(f,ngp);
834  // Now we find the offset from ngp to next-lowest grip point (lgp)
835  for (d=0; d<3U; ++d) {
836  if (gpos(d) > ppos(d)) {
837  lgpoff[d] = -1; // save the offset
838  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
839  }
840  else {
841  lgpoff[d] = 0; // save the offset
842  }
843  }
844  // adjust position of Field iterator to lgp position
845  fiter.moveBy(lgpoff);
846  // get distance from ppos to gpos
847  dpos = ppos - gpos;
848  // normalize dpos by mesh spacing
849  dpos /= delta;
850  // accumulate into local elements
851  *fiter += (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * pdata;
852  fiter.offset(1,0,0) += dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * pdata;
853  fiter.offset(0,1,0) += (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * pdata;
854  fiter.offset(1,1,0) += dpos(0) * dpos(1) * (1 - dpos(2)) * pdata;
855  fiter.offset(0,0,1) += (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * pdata;
856  fiter.offset(1,0,1) += dpos(0) * (1 - dpos(1)) * dpos(2) * pdata;
857  fiter.offset(0,1,1) += (1 - dpos(0)) * dpos(1) * dpos(2) * pdata;
858  fiter.offset(1,1,1) += dpos(0) * dpos(1) * dpos(2) * pdata;
859  return;
860  }
861 
862  // scatter particle data into Field using particle position and mesh
863  // and cache mesh information for reuse
864  template <class FT, class M, class C, class PT>
865  static
866  void scatter(const FT& pdata, Field<FT,3U,M,C>& f,
867  const Vektor<PT,3U>& ppos, const M& mesh,
868  NDIndex<3U>& ngp, int lgpoff[3U], Vektor<PT,3U>& dpos) {
869  CenteringTag<C> ctag;
870  Vektor<PT,3U> gpos, delta;
871  CompressedBrickIterator<FT,3U> fiter;
872  unsigned int d;
873  // find nearest grid point for particle position, store in NDIndex obj
874  ngp = FindNGP(mesh, ppos, ctag);
875  // get position of ngp
876  FindPos(gpos, mesh, ngp, ctag);
877  // get mesh spacings
878  FindDelta(delta, mesh, ngp, ctag);
879  // Try to find ngp in local fields and get iterator
880  fiter = getFieldIter(f,ngp);
881  // Now we find the offset from ngp to next-lowest grip point (lgp)
882  for (d=0; d<3U; ++d) {
883  if (gpos(d) > ppos(d)) {
884  lgpoff[d] = -1; // save the offset
885  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
886  }
887  else {
888  lgpoff[d] = 0; // save the offset
889  }
890  }
891  // adjust position of Field iterator to lgp position
892  fiter.moveBy(lgpoff);
893  // get distance from ppos to gpos
894  dpos = ppos - gpos;
895  // normalize dpos by mesh spacing
896  dpos /= delta;
897  // accumulate into local elements
898  *fiter += (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * pdata;
899  fiter.offset(1,0,0) += dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * pdata;
900  fiter.offset(0,1,0) += (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * pdata;
901  fiter.offset(1,1,0) += dpos(0) * dpos(1) * (1 - dpos(2)) * pdata;
902  fiter.offset(0,0,1) += (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * pdata;
903  fiter.offset(1,0,1) += dpos(0) * (1 - dpos(1)) * dpos(2) * pdata;
904  fiter.offset(0,1,1) += (1 - dpos(0)) * dpos(1) * dpos(2) * pdata;
905  fiter.offset(1,1,1) += dpos(0) * dpos(1) * dpos(2) * pdata;
906  return;
907  }
908 
909  // scatter particle data into Field using cached mesh information
910  template <class FT, class M, class C, class PT>
911  static
912  void scatter(const FT& pdata, Field<FT,3U,M,C>& f,
913  const NDIndex<3U>& ngp, const int lgpoff[3U],
914  const Vektor<PT,3U>& dpos) {
915  CompressedBrickIterator<FT,3U> fiter;
916  // Try to find ngp in local fields and get iterator
917  fiter = getFieldIter(f,ngp);
918  // adjust position of Field iterator to lgp position
919  fiter.moveBy(lgpoff);
920  // accumulate into local elements
921  *fiter += (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * pdata;
922  fiter.offset(1,0,0) += dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * pdata;
923  fiter.offset(0,1,0) += (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * pdata;
924  fiter.offset(1,1,0) += dpos(0) * dpos(1) * (1 - dpos(2)) * pdata;
925  fiter.offset(0,0,1) += (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * pdata;
926  fiter.offset(1,0,1) += dpos(0) * (1 - dpos(1)) * dpos(2) * pdata;
927  fiter.offset(0,1,1) += (1 - dpos(0)) * dpos(1) * dpos(2) * pdata;
928  fiter.offset(1,1,1) += dpos(0) * dpos(1) * dpos(2) * pdata;
929  return;
930  }
931 
932  // gather particle data from Field using particle position and mesh
933  template <class FT, class M, class C, class PT>
934  static
935  void gather(FT& pdata, const Field<FT,3U,M,C>& f,
936  const Vektor<PT,3U>& ppos, const M& mesh) {
937  CenteringTag<C> ctag;
938  Vektor<PT,3U> gpos, dpos, delta;
939  NDIndex<3U> ngp;
940  CompressedBrickIterator<FT,3U> fiter;
941  int lgpoff[3U];
942  unsigned int d;
943  // find nearest grid point for particle position, store in NDIndex obj
944  ngp = FindNGP(mesh, ppos, ctag);
945  // get position of ngp
946  FindPos(gpos, mesh, ngp, ctag);
947  // get mesh spacings
948  FindDelta(delta, mesh, ngp, ctag);
949  // Try to find ngp in local fields and get iterator
950  fiter = getFieldIter(f,ngp);
951  // Now we find the offset from ngp to next-lowest grip point (lgp)
952  for (d=0; d<3U; ++d) {
953  if (gpos(d) > ppos(d)) {
954  lgpoff[d] = -1; // save the offset
955  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
956  }
957  else {
958  lgpoff[d] = 0; // save the offset
959  }
960  }
961  // adjust position of Field iterator to lgp position
962  fiter.moveBy(lgpoff);
963  // get distance from ppos to gpos
964  dpos = ppos - gpos;
965  // normalize dpos by mesh spacing
966  dpos /= delta;
967  // accumulate into particle attrib
968  pdata = (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * (*fiter)
969  + dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * fiter.offset(1,0,0)
970  + (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * fiter.offset(0,1,0)
971  + dpos(0) * dpos(1) * (1 - dpos(2)) * fiter.offset(1,1,0)
972  + (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * fiter.offset(0,0,1)
973  + dpos(0) * (1 - dpos(1)) * dpos(2) * fiter.offset(1,0,1)
974  + (1 - dpos(0)) * dpos(1) * dpos(2) * fiter.offset(0,1,1)
975  + dpos(0) * dpos(1) * dpos(2) * fiter.offset(1,1,1);
976  return;
977  }
978 
979  // gather particle data from Field using particle position and mesh
980  // and cache mesh information for reuse
981  template <class FT, class M, class C, class PT>
982  static
983  void gather(FT& pdata, const Field<FT,3U,M,C>& f,
984  const Vektor<PT,3U>& ppos, const M& mesh,
985  NDIndex<3U>& ngp, int lgpoff[3U], Vektor<PT,3U>& dpos) {
986  CenteringTag<C> ctag;
987  Vektor<PT,3U> gpos, delta;
988  CompressedBrickIterator<FT,3U> fiter;
989  unsigned int d;
990  // find nearest grid point for particle position, store in NDIndex obj
991  ngp = FindNGP(mesh, ppos, ctag);
992  // get position of ngp
993  FindPos(gpos, mesh, ngp, ctag);
994  // get mesh spacings
995  FindDelta(delta, mesh, ngp, ctag);
996  // Try to find ngp in local fields and get iterator
997  fiter = getFieldIter(f,ngp);
998  // Now we find the offset from ngp to next-lowest grip point (lgp)
999  for (d=0; d<3U; ++d) {
1000  if (gpos(d) > ppos(d)) {
1001  lgpoff[d] = -1; // save the offset
1002  gpos(d) = gpos(d) - delta(d); // adjust gpos to lgp position
1003  }
1004  else {
1005  lgpoff[d] = 0; // save the offset
1006  }
1007  }
1008  // adjust position of Field iterator to lgp position
1009  fiter.moveBy(lgpoff);
1010  // get distance from ppos to gpos
1011  dpos = ppos - gpos;
1012  // normalize dpos by mesh spacing
1013  dpos /= delta;
1014  // accumulate into particle attrib
1015  pdata = (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * (*fiter)
1016  + dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * fiter.offset(1,0,0)
1017  + (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * fiter.offset(0,1,0)
1018  + dpos(0) * dpos(1) * (1 - dpos(2)) * fiter.offset(1,1,0)
1019  + (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * fiter.offset(0,0,1)
1020  + dpos(0) * (1 - dpos(1)) * dpos(2) * fiter.offset(1,0,1)
1021  + (1 - dpos(0)) * dpos(1) * dpos(2) * fiter.offset(0,1,1)
1022  + dpos(0) * dpos(1) * dpos(2) * fiter.offset(1,1,1);
1023  return;
1024  }
1025 
1026  // gather particle data from Field using cached mesh information
1027  template <class FT, class M, class C, class PT>
1028  static
1029  void gather(FT& pdata, const Field<FT,3U,M,C>& f,
1030  const NDIndex<3U>& ngp, const int lgpoff[3U],
1031  const Vektor<PT,3U>& dpos) {
1032  CompressedBrickIterator<FT,3U> fiter;
1033  // Try to find ngp in local fields and get iterator
1034  fiter = getFieldIter(f,ngp);
1035  // adjust position of Field iterator to lgp position
1036  fiter.moveBy(lgpoff);
1037  // accumulate into particle attrib
1038  pdata = (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * (*fiter)
1039  + dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * fiter.offset(1,0,0)
1040  + (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * fiter.offset(0,1,0)
1041  + dpos(0) * dpos(1) * (1 - dpos(2)) * fiter.offset(1,1,0)
1042  + (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * fiter.offset(0,0,1)
1043  + dpos(0) * (1 - dpos(1)) * dpos(2) * fiter.offset(1,0,1)
1044  + (1 - dpos(0)) * dpos(1) * dpos(2) * fiter.offset(0,1,1)
1045  + dpos(0) * dpos(1) * dpos(2) * fiter.offset(1,1,1);
1046  return;
1047  }
1048 
1049  // scatter particle data into Field using particle position and mesh
1050  template <class FT, class M, class PT>
1051  static
1052  void scatter(const Vektor<FT,3U>& pdata, Field<Vektor<FT,3U>,3U,M,Edge>& f,
1053  const Vektor<PT,3U>& ppos, const M& mesh) {
1054  CenteringTag<Cell> ctag;
1055  Vektor<PT,3U> ppos2, dpos, gpos, delta;
1056  NDIndex<3U> ngp;
1057 
1058  CompressedBrickIterator<Vektor<FT,3U>,3U> fiter;
1059 
1060  // find nearest grid point for particle position, store in NDIndex obj
1061  ngp = FindNGP(mesh, ppos, ctag);
1062  FindDelta(delta, mesh, ngp, ctag);
1063 
1064  for (unsigned int comp = 0; comp < 3U; ++comp) {
1065  Vektor<PT,3U> cppos = ppos;
1066  cppos(comp) += 0.5 * delta(comp);
1067  ngp = FindNGP(mesh, cppos, ctag);
1068  ngp[comp] = ngp[comp] - 1;
1069 
1070  FindPos(gpos, mesh, ngp, CenteringTag<Vert>());
1071  gpos(comp) += 0.5 * delta(comp);
1072  dpos = ppos - gpos;
1073  dpos /= delta;
1074 
1075  // Try to find ngp in local fields and get iterator
1076  fiter = Interpolator::getFieldIter(f,ngp);
1077 
1078  // accumulate into local elements
1079  (*fiter)(comp) += (1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * pdata(comp);
1080  (fiter.offset(1,0,0))(comp) += dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * pdata(comp);
1081  (fiter.offset(0,1,0))(comp) += (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * pdata(comp);
1082  (fiter.offset(1,1,0))(comp) += dpos(0) * dpos(1) * (1 - dpos(2)) * pdata(comp);
1083  (fiter.offset(0,0,1))(comp) += (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * pdata(comp);
1084  (fiter.offset(1,0,1))(comp) += dpos(0) * (1 - dpos(1)) * dpos(2) * pdata(comp);
1085  (fiter.offset(0,1,1))(comp) += (1 - dpos(0)) * dpos(1) * dpos(2) * pdata(comp);
1086  (fiter.offset(1,1,1))(comp) += dpos(0) * dpos(1) * dpos(2) * pdata(comp);
1087  }
1088  return;
1089  }
1090 
1091  template <class FT, class M, class PT>
1092  static
1093  void scatter(const FT& /*pdata*/, Field<FT,3U,M,Edge>& /*f*/,
1094  const Vektor<PT,3U>& /*ppos*/, const M& /*mesh*/) {
1095  ERRORMSG("IntCIC::scatter on Edge based field: not implemented for non-vectors!!"<<endl);
1096  return;
1097  }
1098  // scatter particle data into Field using particle position and mesh
1099  // and cache mesh information for reuse
1100  template <class FT, class M, class PT>
1101  static
1102  void scatter(const FT& /*pdata*/, Field<FT,3U,M,Edge>& /*f*/,
1103  const Vektor<PT,3U>& /*ppos*/, const M& /*mesh*/,
1104  NDIndex<3U>& /*ngp*/, int /*lgpoff*/[3U], Vektor<PT,3U>& /*dpos*/) {
1105  ERRORMSG("IntCIC::scatter on Edge based field with cache: not implemented!!"<<endl);
1106  return;
1107  }
1108 
1109  template <class FT, class M, class PT>
1110  static
1111  void scatter(const FT& /*pdata*/, Field<FT,3U,M,Edge>& /*f*/,
1112  const NDIndex<3U>& /*ngp*/, const int /*lgpoff*/[3U],
1113  const Vektor<PT,3U>& /*dpos*/) {
1114  ERRORMSG("IntCIC::scatter on Edge based field with cache: not implemented!!"<<endl);
1115  return;
1116  }
1117 
1118  // gather particle data from Field using particle position and mesh
1119  template <class FT, class M, class PT>
1120  static
1121  void gather(Vektor<FT,3U>& pdata, const Field<Vektor<FT,3U>,3U,M,Edge>& f,
1122  const Vektor<PT,3U>& ppos, const M& mesh) {
1123  Vektor<PT,3U> dpos, gpos, delta;
1124  NDIndex<3U> ngp = mesh.getNearestVertex(ppos);;
1125  CompressedBrickIterator<Vektor<FT,3U>,3U> fiter;
1126 
1127  FindDelta(delta, mesh, ngp, CenteringTag<Vert>());
1128 
1129  for (unsigned int comp = 0; comp < 3U; ++ comp) {
1130  Vektor<PT,3U> cppos = ppos;
1131  cppos(comp) += 0.5 * delta(comp);
1132  ngp = mesh.getCellContaining(cppos);
1133  ngp[comp] = ngp[comp] - 1;
1134  gpos = mesh.getVertexPosition(ngp);
1135  gpos(comp) += 0.5 * delta(comp);
1136 
1137  // get distance from ppos to gpos
1138  dpos = ppos - gpos;
1139  // normalize dpos by mesh spacing
1140  dpos /= delta;
1141  // accumulate into local elements
1142 
1143  // Try to find ngp in local fields and get iterator
1144  fiter = Interpolator::getFieldIter(f,ngp);
1145 
1146  // accumulate into particle attrib
1147  pdata(comp) = ((1 - dpos(0)) * (1 - dpos(1)) * (1 - dpos(2)) * (*fiter)(comp) +
1148  dpos(0) * (1 - dpos(1)) * (1 - dpos(2)) * (fiter.offset(1,0,0))(comp) +
1149  (1 - dpos(0)) * dpos(1) * (1 - dpos(2)) * (fiter.offset(0,1,0))(comp) +
1150  dpos(0) * dpos(1) * (1 - dpos(2)) * (fiter.offset(1,1,0))(comp) +
1151  (1 - dpos(0)) * (1 - dpos(1)) * dpos(2) * (fiter.offset(0,0,1))(comp) +
1152  dpos(0) * (1 - dpos(1)) * dpos(2) * (fiter.offset(1,0,1))(comp) +
1153  (1 - dpos(0)) * dpos(1) * dpos(2) * (fiter.offset(0,1,1))(comp) +
1154  dpos(0) * dpos(1) * dpos(2) * (fiter.offset(1,1,1))(comp));
1155  }
1156  return;
1157  }
1158 
1159  template <class FT, class M, class PT>
1160  static
1161  void gather(FT& /*pdata*/, const Field<FT,3U,M,Edge>& /*f*/,
1162  const Vektor<PT,3U>& /*ppos*/, const M& /*mesh*/) {
1163  ERRORMSG("IntCIC::gather on Edge based field: not implemented for non-vectors!!"<<endl);
1164  return;
1165  }
1166 
1167  // gather particle data from Field using particle position and mesh
1168  // and cache mesh information for reuse
1169  template <class FT, class M, class PT>
1170  static
1171  void gather(FT& /*pdata*/, const Field<FT,3U,M,Edge>& /*f*/,
1172  const Vektor<PT,3U>& /*ppos*/, const M& /*mesh*/,
1173  NDIndex<3U>& /*ngp*/, int /*lgpoff*/[3U], Vektor<PT,3U>& /*dpos*/) {
1174  ERRORMSG("IntCIC::gather on Edge based field with cache: not implemented!!"<<endl);
1175  return;
1176  }
1177 
1178  // gather particle data from Field using cached mesh information
1179  template <class FT, class M, class PT>
1180  static
1181  void gather(FT& /*pdata*/, const Field<FT,3U,M,Edge>& /*f*/,
1182  const NDIndex<3U>& /*ngp*/, const int /*lgpoff*/[3U],
1183  const Vektor<PT,3U>& /*dpos*/) {
1184  ERRORMSG("IntCIC::gather on Edge based field with cache: not implemented!!"<<endl);
1185  return;
1186  }
1187 };
1188 
1189 
1190 // IntCIC class -- what the user sees
1191 class IntCIC {
1192 
1193 public:
1194  // constructor/destructor
1195  IntCIC() {}
1196  ~IntCIC() {}
1197 
1198  // gather/scatter functions
1199 
1200  // scatter particle data into Field using particle position and mesh
1201  template <class FT, unsigned Dim, class M, class C, class PT>
1202  static
1203  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
1204  const Vektor<PT,Dim>& ppos, const M& mesh) {
1205  IntCICImpl<Dim>::scatter(pdata,f,ppos,mesh);
1206  }
1207 
1208  // scatter particle data into Field using particle position and mesh
1209  // and cache mesh information for reuse
1210  template <class FT, unsigned Dim, class M, class C, class PT>
1211  static
1212  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
1213  const Vektor<PT,Dim>& ppos, const M& mesh,
1214  CacheDataCIC<PT,Dim>& cache) {
1215  IntCICImpl<Dim>::scatter(pdata,f,ppos,mesh,cache.Index_m,
1216  cache.Offset_m,cache.Delta_m);
1217  }
1218 
1219  // scatter particle data into Field using cached mesh information
1220  template <class FT, unsigned Dim, class M, class C, class PT>
1221  static
1222  void scatter(const FT& pdata, Field<FT,Dim,M,C>& f,
1223  const CacheDataCIC<PT,Dim>& cache) {
1224  IntCICImpl<Dim>::scatter(pdata,f,cache.Index_m,
1225  cache.Offset_m,cache.Delta_m);
1226  }
1227 
1228  // gather particle data from Field using particle position and mesh
1229  template <class FT, unsigned Dim, class M, class C, class PT>
1230  static
1231  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
1232  const Vektor<PT,Dim>& ppos, const M& mesh) {
1233  IntCICImpl<Dim>::gather(pdata,f,ppos,mesh);
1234  }
1235 
1236  // gather particle data from Field using particle position and mesh
1237  // and cache mesh information for reuse
1238  template <class FT, unsigned Dim, class M, class C, class PT>
1239  static
1240  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
1241  const Vektor<PT,Dim>& ppos, const M& mesh,
1242  CacheDataCIC<PT,Dim>& cache) {
1243  IntCICImpl<Dim>::gather(pdata,f,ppos,mesh,cache.Index_m,
1244  cache.Offset_m,cache.Delta_m);
1245  }
1246 
1247  // gather particle data from Field using cached mesh information
1248  template <class FT, unsigned Dim, class M, class C, class PT>
1249  static
1250  void gather(FT& pdata, const Field<FT,Dim,M,C>& f,
1251  const CacheDataCIC<PT,Dim>& cache) {
1252  IntCICImpl<Dim>::gather(pdata,f,cache.Index_m,
1253  cache.Offset_m,cache.Delta_m);
1254  }
1255 
1256 };
1257 
1258 #endif // INT_CIC_H
1259 
Dim
const unsigned Dim
Definition: P3MPoissonSolver.h:25
FindNGP
NDIndex< Dim > FindNGP(const M &mesh, const Vektor< PT, Dim > &ppos, CenteringTag< Cell >)
Definition: Interpolator.h:37
FindDelta
void FindDelta(Vektor< PT, Dim > &delta, const M &mesh, const NDIndex< Dim > &gp, CenteringTag< Cell >)
Definition: Interpolator.h:95
FindPos
void FindPos(Vektor< PT, Dim > &pos, const M &mesh, const NDIndex< Dim > &indices, CenteringTag< Cell >)
Definition: Interpolator.h:62
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
ERRORMSG
#define ERRORMSG(msg)
Definition: IpplInfo.h:350
Vektor
Definition: Vektor.h:32
Field
Definition: Field.h:33
BrickIterator::offset
T & offset(int i) const
Definition: BrickIterator.h:107
BrickIterator::moveBy
void moveBy(int i)
Definition: BrickIterator.h:135
Edge
Definition: Centering.h:50
InterpolatorTraits< T, Dim, IntCIC >::Cache_t
CacheDataCIC< T, Dim > Cache_t
Definition: IntCIC.h:30
IntCICImpl::~IntCICImpl
~IntCICImpl()
Definition: IntCIC.h:42
IntCICImpl::IntCICImpl
IntCICImpl()
Definition: IntCIC.h:41
IntCICImpl::gather
static void gather(FT &, const Field< FT, Dim, M, C > &f, const NDIndex< Dim > &ngp, const int lgpoff[Dim], const Vektor< PT, Dim > &)
Definition: IntCIC.h:225
IntCICImpl::gather
static void gather(FT &, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntCIC.h:145
IntCICImpl::gather
static void gather(FT &, 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: IntCIC.h:186
IntCICImpl::scatter
static void scatter(const FT &, 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: IntCIC.h:90
IntCICImpl::scatter
static void scatter(const FT &, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntCIC.h:49
IntCICImpl::scatter
static void scatter(const FT &, Field< FT, Dim, M, C > &f, const NDIndex< Dim > &ngp, const int lgpoff[Dim], const Vektor< PT, Dim > &)
Definition: IntCIC.h:129
IntCICImpl< 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: IntCIC.h:336
IntCICImpl< 1U >::gather
static void gather(FT &pdata, const Field< FT, 1U, M, C > &f, const Vektor< PT, 1U > &ppos, const M &mesh)
Definition: IntCIC.h:353
IntCICImpl< 1U >::scatter
static void scatter(const FT &pdata, Field< FT, 1U, M, C > &f, const Vektor< PT, 1U > &ppos, const M &mesh)
Definition: IntCIC.h:254
IntCICImpl< 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: IntCIC.h:395
IntCICImpl< 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: IntCIC.h:435
IntCICImpl< 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: IntCIC.h:296
IntCICImpl< 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: IntCIC.h:614
IntCICImpl< 2U >::gather
static void gather(FT &, const Field< FT, 2U, M, Edge > &, const Vektor< PT, 2U > &, const M &, NDIndex< 2U > &, int[2U], Vektor< PT, 2U > &)
Definition: IntCIC.h:786
IntCICImpl< 2U >::scatter
static void scatter(const FT &, Field< FT, 2U, M, Edge > &, const Vektor< PT, 2U > &, const M &)
Definition: IntCIC.h:712
IntCICImpl< 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: IntCIC.h:656
IntCICImpl< 2U >::scatter
static void scatter(const Vektor< FT, 2U > &pdata, Field< Vektor< FT, 2U >, 2U, M, Edge > &f, const Vektor< PT, 2U > &ppos, const M &mesh)
Definition: IntCIC.h:675
IntCICImpl< 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: IntCIC.h:551
IntCICImpl< 2U >::gather
static void gather(FT &, const Field< FT, 2U, M, Edge > &, const NDIndex< 2U > &, const int[2U], const Vektor< PT, 2U > &)
Definition: IntCIC.h:796
IntCICImpl< 2U >::gather
static void gather(FT &, const Field< FT, 2U, M, Edge > &, const Vektor< PT, 2U > &, const M &)
Definition: IntCIC.h:776
IntCICImpl< 2U >::gather
static void gather(Vektor< FT, 2U > &pdata, const Field< Vektor< FT, 2U >, 2U, M, Edge > &f, const Vektor< PT, 2U > &ppos, const M &mesh)
Definition: IntCIC.h:740
IntCICImpl< 2U >::scatter
static void scatter(const FT &pdata, Field< FT, 2U, M, C > &f, const Vektor< PT, 2U > &ppos, const M &mesh)
Definition: IntCIC.h:465
IntCICImpl< 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: IntCIC.h:509
IntCICImpl< 2U >::scatter
static void scatter(const FT &, Field< FT, 2U, M, Edge > &, const Vektor< PT, 2U > &, const M &, NDIndex< 2U > &, int[2U], Vektor< PT, 2U > &)
Definition: IntCIC.h:721
IntCICImpl< 2U >::gather
static void gather(FT &pdata, const Field< FT, 2U, M, C > &f, const Vektor< PT, 2U > &ppos, const M &mesh)
Definition: IntCIC.h:570
IntCICImpl< 2U >::scatter
static void scatter(const FT &, Field< FT, 2U, M, Edge > &, const NDIndex< 2U > &, const int[2U], const Vektor< PT, 2U > &)
Definition: IntCIC.h:730
IntCICImpl< 3U >::scatter
static void scatter(const FT &, Field< FT, 3U, M, Edge > &, const Vektor< PT, 3U > &, const M &)
Definition: IntCIC.h:1093
IntCICImpl< 3U >::gather
static void gather(FT &pdata, const Field< FT, 3U, M, C > &f, const Vektor< PT, 3U > &ppos, const M &mesh)
Definition: IntCIC.h:935
IntCICImpl< 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: IntCIC.h:983
IntCICImpl< 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: IntCIC.h:1029
IntCICImpl< 3U >::gather
static void gather(FT &, const Field< FT, 3U, M, Edge > &, const NDIndex< 3U > &, const int[3U], const Vektor< PT, 3U > &)
Definition: IntCIC.h:1181
IntCICImpl< 3U >::scatter
static void scatter(const FT &pdata, Field< FT, 3U, M, C > &f, const Vektor< PT, 3U > &ppos, const M &mesh)
Definition: IntCIC.h:818
IntCICImpl< 3U >::gather
static void gather(Vektor< FT, 3U > &pdata, const Field< Vektor< FT, 3U >, 3U, M, Edge > &f, const Vektor< PT, 3U > &ppos, const M &mesh)
Definition: IntCIC.h:1121
IntCICImpl< 3U >::scatter
static void scatter(const Vektor< FT, 3U > &pdata, Field< Vektor< FT, 3U >, 3U, M, Edge > &f, const Vektor< PT, 3U > &ppos, const M &mesh)
Definition: IntCIC.h:1052
IntCICImpl< 3U >::scatter
static void scatter(const FT &, Field< FT, 3U, M, Edge > &, const Vektor< PT, 3U > &, const M &, NDIndex< 3U > &, int[3U], Vektor< PT, 3U > &)
Definition: IntCIC.h:1102
IntCICImpl< 3U >::gather
static void gather(FT &, const Field< FT, 3U, M, Edge > &, const Vektor< PT, 3U > &, const M &)
Definition: IntCIC.h:1161
IntCICImpl< 3U >::gather
static void gather(FT &, const Field< FT, 3U, M, Edge > &, const Vektor< PT, 3U > &, const M &, NDIndex< 3U > &, int[3U], Vektor< PT, 3U > &)
Definition: IntCIC.h:1171
IntCICImpl< 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: IntCIC.h:912
IntCICImpl< 3U >::scatter
static void scatter(const FT &, Field< FT, 3U, M, Edge > &, const NDIndex< 3U > &, const int[3U], const Vektor< PT, 3U > &)
Definition: IntCIC.h:1111
IntCICImpl< 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: IntCIC.h:866
IntCIC::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntCIC.h:1203
IntCIC::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh)
Definition: IntCIC.h:1231
IntCIC::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const CacheDataCIC< PT, Dim > &cache)
Definition: IntCIC.h:1222
IntCIC::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh, CacheDataCIC< PT, Dim > &cache)
Definition: IntCIC.h:1240
IntCIC::gather
static void gather(FT &pdata, const Field< FT, Dim, M, C > &f, const CacheDataCIC< PT, Dim > &cache)
Definition: IntCIC.h:1250
IntCIC::scatter
static void scatter(const FT &pdata, Field< FT, Dim, M, C > &f, const Vektor< PT, Dim > &ppos, const M &mesh, CacheDataCIC< PT, Dim > &cache)
Definition: IntCIC.h:1212
IntCIC::IntCIC
IntCIC()
Definition: IntCIC.h:1195
IntCIC::~IntCIC
~IntCIC()
Definition: IntCIC.h:1196
CacheDataCIC::Index_m
NDIndex< Dim > Index_m
Definition: Interpolator.h:148
CacheDataCIC::Offset_m
int Offset_m[Dim]
Definition: Interpolator.h:149
CacheDataCIC::Delta_m
Vektor< T, Dim > Delta_m
Definition: Interpolator.h:150
Interpolator::getFieldIter
static CompressedBrickIterator< T, Dim > getFieldIter(const BareField< T, Dim > &f, const NDIndex< Dim > &pt)
Definition: Interpolator.h:183
Generated on Wed Aug 25 2021 16:40:59 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.1