OPAL (Object Oriented Parallel Accelerator Library) 2022.1
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
25class IntCIC;
26
27// specialization of InterpolatorTraits
28template <class T, unsigned Dim>
29struct InterpolatorTraits<T,Dim,IntCIC> {
30 typedef CacheDataCIC<T,Dim> Cache_t;
31};
32
33
34
35// IntCICImpl class definition
36template <unsigned Dim>
37class IntCICImpl : public Interpolator {
38
39public:
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
241template <>
242class IntCICImpl<1U> : public Interpolator {
243
244public:
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
452template <>
453class IntCICImpl<2U> : public Interpolator {
454
455public:
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
805template <>
806class IntCICImpl<3U> : public Interpolator {
807
808public:
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
1191class IntCIC {
1192
1193public:
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:26
FindDelta
void FindDelta(Vektor< PT, Dim > &delta, const M &mesh, const NDIndex< Dim > &gp, CenteringTag< Cell >)
Definition: Interpolator.h:95
FindNGP
NDIndex< Dim > FindNGP(const M &mesh, const Vektor< PT, Dim > &ppos, CenteringTag< Cell >)
Definition: Interpolator.h:37
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 Thu Oct 20 2022 17:40:08 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.3