OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
UniformCartesian.hpp
Go to the documentation of this file.
1 // -*- C++ -*-
2 /***************************************************************************
3  *
4  * The IPPL Framework
5  *
6  ***************************************************************************/
7 
8 // UniformCartesian.cpp
9 // Implementations for UniformCartesian mesh class (uniform spacings)
10 
11 // include files
12 #include "Utility/PAssert.h"
13 #include "Utility/IpplInfo.h"
14 #include "Field/BareField.h"
15 #include "Field/BrickExpression.h"
16 #include "Field/LField.h"
17 #include "Field/Assign.h"
18 #include "Field/AssignDefs.h"
19 
20 //-----------------------------------------------------------------------------
21 // Setup chores common to all constructors:
22 //-----------------------------------------------------------------------------
23 template <unsigned Dim, class MFLOAT>
24 void
25 UniformCartesian<Dim,MFLOAT>::
26 setup()
27 {
28  hasSpacingFields = false;
29  FlCell = 0;
30  FlVert = 0;
31  VertSpacings = 0;
32  CellSpacings = 0;
33  volume = 0.0;
34 }
35 
36 //-----------------------------------------------------------------------------
37 // Constructors from NDIndex object:
38 //-----------------------------------------------------------------------------
39 template <unsigned Dim, class MFLOAT>
40 UniformCartesian<Dim,MFLOAT>::
41 UniformCartesian(const NDIndex<Dim>& ndi)
42 {
43  setup();
44  for (unsigned int d=0; d<Dim; d++) {
45  gridSizes[d] = ndi[d].length(); // Number of vertices along this dimension.
46  meshSpacing[d] = ndi[d].stride(); // Default mesh spacing from stride()
47  origin(d) = ndi[d].first(); // Default origin at ndi[d].first
48  }
49  volume = 1.0; // Default mesh has unit cell volume.
50  set_Dvc(); // Set derivative coefficients from spacings.
51 }
52 // Also specify mesh spacings:
53 template <unsigned Dim, class MFLOAT>
54 UniformCartesian<Dim,MFLOAT>::
55 UniformCartesian(const NDIndex<Dim>& ndi, MFLOAT* const delX)
56 {
57  setup();
58  for (unsigned int d=0; d<Dim; d++) {
59  gridSizes[d] = ndi[d].length(); // Number of vertices along this dimension.
60  origin(d) = ndi[d].first(); // Default origin at ndi[d].first
61  }
62  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
63 }
64 // Also specify mesh spacings and origin:
65 template <unsigned Dim, class MFLOAT>
66 UniformCartesian<Dim,MFLOAT>::
67 UniformCartesian(const NDIndex<Dim>& ndi, MFLOAT* const delX,
68  const Vektor<MFLOAT,Dim>& orig)
69 {
70  setup();
71  for (unsigned int d=0; d<Dim; d++) {
72  gridSizes[d] = ndi[d].length(); // Number of vertices along this dimension.
73  }
74  set_origin(orig); // Set origin.
75  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
76 }
77 //-----------------------------------------------------------------------------
78 // Constructors from Index objects:
79 //-----------------------------------------------------------------------------
80 
81 //===========1D============
82 template <unsigned Dim, class MFLOAT>
83 UniformCartesian<Dim,MFLOAT>::
84 UniformCartesian(const Index& I)
85 {
86  PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
87  setup();
88  gridSizes[0] = I.length(); // Number of vertices along this dimension.
89  meshSpacing[0] = I.stride(); // Default mesh spacing from stride()
90  origin(0) = I.first(); // Default origin at I.first()
91 
92  volume = 1.0; // Default mesh has unit cell volume.
93  set_Dvc(); // Set derivative coefficients from spacings.
94 }
95 // Also specify mesh spacings:
96 template <unsigned Dim, class MFLOAT>
97 UniformCartesian<Dim,MFLOAT>::
98 UniformCartesian(const Index& I, MFLOAT* const delX)
99 {
100  PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
101  setup();
102  gridSizes[0] = I.length(); // Number of vertices along this dimension.
103  origin(0) = I.first(); // Default origin at I.first()
104 
105  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
106 }
107 // Also specify mesh spacings and origin:
108 template <unsigned Dim, class MFLOAT>
109 UniformCartesian<Dim,MFLOAT>::
110 UniformCartesian(const Index& I, MFLOAT* const delX,
111  const Vektor<MFLOAT,Dim>& orig)
112 {
113  PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
114  setup();
115  gridSizes[0] = I.length(); // Number of vertices along this dimension.
116  set_origin(orig); // Set origin.
117  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
118 }
119 
120 //===========2D============
121 template <unsigned Dim, class MFLOAT>
122 UniformCartesian<Dim,MFLOAT>::
123 UniformCartesian(const Index& I, const Index& J)
124 {
125  PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
126  setup();
127  gridSizes[0] = I.length(); // Number of vertices along this dimension.
128  gridSizes[1] = J.length(); // Number of vertices along this dimension.
129  meshSpacing[0] = I.stride(); // Default mesh spacing from stride()
130  meshSpacing[1] = J.stride();
131  origin(0) = I.first(); // Default origin at (I.first(),J.first())
132  origin(1) = J.first();
133 
134  volume = 1.0; // Default mesh has unit cell volume.
135  set_Dvc(); // Set derivative coefficients from spacings.
136 }
137 // Also specify mesh spacings:
138 template <unsigned Dim, class MFLOAT>
139 UniformCartesian<Dim,MFLOAT>::
140 UniformCartesian(const Index& I, const Index& J, MFLOAT* const delX)
141 {
142  PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
143  setup();
144  gridSizes[0] = I.length(); // Number of vertices along this dimension.
145  gridSizes[1] = J.length(); // Number of vertices along this dimension.
146  origin(0) = I.first(); // Default origin at (I.first(),J.first())
147  origin(1) = J.first();
148  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
149 }
150 // Also specify mesh spacings and origin:
151 template <unsigned Dim, class MFLOAT>
152 UniformCartesian<Dim,MFLOAT>::
153 UniformCartesian(const Index& I, const Index& J, MFLOAT* const delX,
154  const Vektor<MFLOAT,Dim>& orig)
155 {
156  PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
157  setup();
158  gridSizes[0] = I.length(); // Number of vertices along this dimension.
159  gridSizes[1] = J.length(); // Number of vertices along this dimension.
160  set_origin(orig); // Set origin.
161  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
162 }
163 
164 //===========3D============
165 template <unsigned Dim, class MFLOAT>
166 UniformCartesian<Dim,MFLOAT>::
167 UniformCartesian(const Index& I, const Index& J, const Index& K)
168 {
169  PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
170  setup();
171  gridSizes[0] = I.length(); // Number of vertices along this dimension.
172  gridSizes[1] = J.length(); // Number of vertices along this dimension.
173  gridSizes[2] = K.length(); // Number of vertices along this dimension.
174  meshSpacing[0] = I.stride(); // Default mesh spacing from stride()
175  meshSpacing[1] = J.stride();
176  meshSpacing[2] = K.stride();
177  origin(0) = I.first(); // Default origin at (I.first(),J.first(),K.first())
178  origin(1) = J.first();
179  origin(2) = K.first();
180 
181  volume = 1.0; // Default mesh has unit cell volume.
182  set_Dvc(); // Set derivative coefficients from spacings.
183 }
184 // Also specify mesh spacings:
185 template <unsigned Dim, class MFLOAT>
186 UniformCartesian<Dim,MFLOAT>::
187 UniformCartesian(const Index& I, const Index& J, const Index& K,
188  MFLOAT* const delX)
189 {
190  PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
191  setup();
192  gridSizes[0] = I.length(); // Number of vertices along this dimension.
193  gridSizes[1] = J.length(); // Number of vertices along this dimension.
194  gridSizes[2] = K.length(); // Number of vertices along this dimension.
195  origin(0) = I.first(); // Default origin at (I.first(),J.first(),K.first())
196  origin(1) = J.first();
197  origin(2) = K.first();
198  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
199 }
200 // Also specify mesh spacings and origin:
201 template <unsigned Dim, class MFLOAT>
202 UniformCartesian<Dim,MFLOAT>::
203 UniformCartesian(const Index& I, const Index& J, const Index& K,
204  MFLOAT* const delX, const Vektor<MFLOAT,Dim>& orig)
205 {
206  PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
207  setup();
208  gridSizes[0] = I.length(); // Number of vertices along this dimension.
209  gridSizes[1] = J.length(); // Number of vertices along this dimension.
210  gridSizes[2] = K.length(); // Number of vertices along this dimension.
211  set_origin(orig); // Set origin.
212  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
213 }
214 
215 //-----------------------------------------------------------------------------
216 // initialize with NDIndex object:
217 //-----------------------------------------------------------------------------
218 template <unsigned Dim, class MFLOAT>
219 void
220 UniformCartesian<Dim,MFLOAT>::
221 initialize(const NDIndex<Dim>& ndi)
222 {
223  setup();
224  for (unsigned int d=0; d<Dim; d++) {
225  gridSizes[d] = ndi[d].length(); // Number of vertices along this dimension.
226  meshSpacing[d] = ndi[d].stride(); // Default mesh spacing from stride()
227  origin(d) = ndi[d].first(); // Default origin at ndi[d].first
228  }
229  volume = 1.0; // Default mesh has unit cell volume.
230  set_Dvc(); // Set derivative coefficients from spacings.
231 }
232 // Also specify mesh spacings:
233 template <unsigned Dim, class MFLOAT>
234 void
235 UniformCartesian<Dim,MFLOAT>::
236 initialize(const NDIndex<Dim>& ndi, MFLOAT* const delX)
237 {
238  setup();
239  for (unsigned int d=0; d<Dim; d++) {
240  gridSizes[d] = ndi[d].length(); // Number of vertices along this dimension.
241  origin(d) = ndi[d].first(); // Default origin at ndi[d].first
242  }
243  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
244 }
245 // Also specify mesh spacings and origin:
246 template <unsigned Dim, class MFLOAT>
247 void
248 UniformCartesian<Dim,MFLOAT>::
249 initialize(const NDIndex<Dim>& ndi, MFLOAT* const delX,
250  const Vektor<MFLOAT,Dim>& orig)
251 {
252  setup();
253  for (unsigned int d=0; d<Dim; d++) {
254  gridSizes[d] = ndi[d].length(); // Number of vertices along this dimension.
255  }
256  set_origin(orig); // Set origin.
257  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
258 }
259 //-----------------------------------------------------------------------------
260 // initialize from Index objects:
261 //-----------------------------------------------------------------------------
262 
263 //===========1D============
264 template <unsigned Dim, class MFLOAT>
265 void
266 UniformCartesian<Dim,MFLOAT>::
267 initialize(const Index& I)
268 {
269  PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
270  setup();
271  gridSizes[0] = I.length(); // Number of vertices along this dimension.
272  meshSpacing[0] = I.stride(); // Default mesh spacing from stride()
273  origin(0) = I.first(); // Default origin at I.first()
274 
275  volume = 1.0; // Default mesh has unit cell volume.
276  set_Dvc(); // Set derivative coefficients from spacings.
277 }
278 // Also specify mesh spacings:
279 template <unsigned Dim, class MFLOAT>
280 void
281 UniformCartesian<Dim,MFLOAT>::
282 initialize(const Index& I, MFLOAT* const delX)
283 {
284  PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
285  setup();
286  gridSizes[0] = I.length(); // Number of vertices along this dimension.
287  origin(0) = I.first(); // Default origin at I.first()
288 
289  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
290 }
291 // Also specify mesh spacings and origin:
292 template <unsigned Dim, class MFLOAT>
293 void
294 UniformCartesian<Dim,MFLOAT>::
295 initialize(const Index& I, MFLOAT* const delX,
296  const Vektor<MFLOAT,Dim>& orig)
297 {
298  PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
299  setup();
300  gridSizes[0] = I.length(); // Number of vertices along this dimension.
301  set_origin(orig); // Set origin.
302  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
303 }
304 
305 //===========2D============
306 template <unsigned Dim, class MFLOAT>
307 void
308 UniformCartesian<Dim,MFLOAT>::
309 initialize(const Index& I, const Index& J)
310 {
311  PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
312  setup();
313  gridSizes[0] = I.length(); // Number of vertices along this dimension.
314  gridSizes[1] = J.length(); // Number of vertices along this dimension.
315  meshSpacing[0] = I.stride(); // Default mesh spacing from stride()
316  meshSpacing[1] = J.stride();
317  origin(0) = I.first(); // Default origin at (I.first(),J.first())
318  origin(1) = J.first();
319 
320  volume = 1.0; // Default mesh has unit cell volume.
321  set_Dvc(); // Set derivative coefficients from spacings.
322 }
323 // Also specify mesh spacings:
324 template <unsigned Dim, class MFLOAT>
325 void
326 UniformCartesian<Dim,MFLOAT>::
327 initialize(const Index& I, const Index& J, MFLOAT* const delX)
328 {
329  PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
330  setup();
331  gridSizes[0] = I.length(); // Number of vertices along this dimension.
332  gridSizes[1] = J.length(); // Number of vertices along this dimension.
333  origin(0) = I.first(); // Default origin at (I.first(),J.first())
334  origin(1) = J.first();
335  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
336 }
337 // Also specify mesh spacings and origin:
338 template <unsigned Dim, class MFLOAT>
339 void
340 UniformCartesian<Dim,MFLOAT>::
341 initialize(const Index& I, const Index& J, MFLOAT* const delX,
342  const Vektor<MFLOAT,Dim>& orig)
343 {
344  PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
345  setup();
346  gridSizes[0] = I.length(); // Number of vertices along this dimension.
347  gridSizes[1] = J.length(); // Number of vertices along this dimension.
348  set_origin(orig); // Set origin.
349  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
350 }
351 
352 //===========3D============
353 template <unsigned Dim, class MFLOAT>
354 void
355 UniformCartesian<Dim,MFLOAT>::
356 initialize(const Index& I, const Index& J, const Index& K)
357 {
358  PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
359  setup();
360  gridSizes[0] = I.length(); // Number of vertices along this dimension.
361  gridSizes[1] = J.length(); // Number of vertices along this dimension.
362  gridSizes[2] = K.length(); // Number of vertices along this dimension.
363  meshSpacing[0] = I.stride(); // Default mesh spacing from stride()
364  meshSpacing[1] = J.stride();
365  meshSpacing[2] = K.stride();
366  origin(0) = I.first(); // Default origin at (I.first(),J.first(),K.first())
367  origin(1) = J.first();
368  origin(2) = K.first();
369 
370  volume = 1.0; // Default mesh has unit cell volume.
371  set_Dvc(); // Set derivative coefficients from spacings.
372 }
373 // Also specify mesh spacings:
374 template <unsigned Dim, class MFLOAT>
375 void
376 UniformCartesian<Dim,MFLOAT>::
377 initialize(const Index& I, const Index& J, const Index& K,
378  MFLOAT* const delX)
379 {
380  PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
381  setup();
382  gridSizes[0] = I.length(); // Number of vertices along this dimension.
383  gridSizes[1] = J.length(); // Number of vertices along this dimension.
384  gridSizes[2] = K.length(); // Number of vertices along this dimension.
385  origin(0) = I.first(); // Default origin at (I.first(),J.first(),K.first())
386  origin(1) = J.first();
387  origin(2) = K.first();
388  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
389 }
390 // Also specify mesh spacings and origin:
391 template <unsigned Dim, class MFLOAT>
392 void
393 UniformCartesian<Dim,MFLOAT>::
394 initialize(const Index& I, const Index& J, const Index& K,
395  MFLOAT* const delX, const Vektor<MFLOAT,Dim>& orig)
396 {
397  PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
398  setup();
399  gridSizes[0] = I.length(); // Number of vertices along this dimension.
400  gridSizes[1] = J.length(); // Number of vertices along this dimension.
401  gridSizes[2] = K.length(); // Number of vertices along this dimension.
402  set_origin(orig); // Set origin.
403  set_meshSpacing(delX); // Set mesh spacings and compute cell volume
404 }
405 
406 //-----------------------------------------------------------------------------
407 // Set/accessor functions for member data:
408 //-----------------------------------------------------------------------------
409 // Set the origin of mesh vertex positions:
410 template<unsigned Dim, class MFLOAT>
411 void UniformCartesian<Dim,MFLOAT>::
412 set_origin(const Vektor<MFLOAT,Dim>& o)
413 {
414  origin = o;
415  this->notifyOfChange();
416 }
417 // Get the origin of mesh vertex positions:
418 template<unsigned Dim, class MFLOAT>
419 Vektor<MFLOAT,Dim> UniformCartesian<Dim,MFLOAT>::
420 get_origin() const
421 {
422  return origin;
423 }
424 
425 // Set the spacings of mesh vertex positions (recompute Dvc, cell volume):
426 template<unsigned Dim, class MFLOAT>
427 void UniformCartesian<Dim,MFLOAT>::
428 set_meshSpacing(MFLOAT* const del)
429 {
430  unsigned d;
431  volume = 1.0;
432  for (d=0; d<Dim; ++d) {
433  meshSpacing[d] = del[d];
434  volume *= del[d];
435  }
436  set_Dvc();
437  // if spacing fields exist, we must recompute values
438  if (hasSpacingFields) storeSpacingFields();
439  this->notifyOfChange();
440 }
441 // Set only the derivative constants, using pre-set spacings:
442 template<unsigned Dim, class MFLOAT>
443 void UniformCartesian<Dim,MFLOAT>::
444 set_Dvc()
445 {
446  unsigned d;
447  MFLOAT coef = 1.0;
448  for (d=1;d<Dim;++d) coef *= 0.5;
449 
450  for (d=0;d<Dim;++d) {
451  MFLOAT dvc = coef/meshSpacing[d];
452  for (unsigned b=0; b<(1u<<Dim); ++b) {
453  int s = ( b&(1<<d) ) ? 1 : -1;
454  Dvc[b][d] = s*dvc;
455  }
456  }
457 }
458 // Get the spacings of mesh vertex positions along specified direction:
459 template<unsigned Dim, class MFLOAT>
460 MFLOAT UniformCartesian<Dim,MFLOAT>::
461 get_meshSpacing(unsigned d) const
462 {
463  PAssert_LT(d, Dim);
464  MFLOAT ms = meshSpacing[d];
465  return ms;
466 }
467 
468 // Get the cell volume:
469 template<unsigned Dim, class MFLOAT>
470 MFLOAT UniformCartesian<Dim,MFLOAT>::
471 get_volume() const
472 {
473  return volume;
474 }
475 
477 
478 // Applicative templates for Mesh BC PETE_apply() functions, used
479 // by BrickExpression in storeSpacingFields()
480 
481 // Reflective/None (all are the same for uniform, incl. Periodic).
482 template<class T>
483 struct OpUMeshExtrapolate
484 {
485  OpUMeshExtrapolate(T& o, T& s) : Offset(o), Slope(s) {}
486  T Offset, Slope;
487 };
488 
489 template<class T>
490 inline void PETE_apply(OpUMeshExtrapolate<T> e, T& a, T b)
491 {
492  a = b*e.Slope+e.Offset;
493 }
494 
496 
497 // Create BareField's of vertex and cell spacings
498 // Special prototypes taking no args or FieldLayout ctor args:
499 // No-arg case:
500 template<unsigned Dim, class MFLOAT>
501 void UniformCartesian<Dim,MFLOAT>::
502 storeSpacingFields()
503 {
504  // Set up default FieldLayout parameters:
505  e_dim_tag et[Dim];
506  for (unsigned int d=0; d<Dim; d++) et[d] = PARALLEL;
507  storeSpacingFields(et, -1);
508 }
509 // 1D
510 template<unsigned Dim, class MFLOAT>
511 void UniformCartesian<Dim,MFLOAT>::
512 storeSpacingFields(e_dim_tag p1, int vnodes)
513 {
514  e_dim_tag et[1];
515  et[0] = p1;
516  storeSpacingFields(et, vnodes);
517 }
518 // 2D
519 template<unsigned Dim, class MFLOAT>
520 void UniformCartesian<Dim,MFLOAT>::
521 storeSpacingFields(e_dim_tag p1, e_dim_tag p2, int vnodes)
522 {
523  e_dim_tag et[2];
524  et[0] = p1;
525  et[1] = p2;
526  storeSpacingFields(et, vnodes);
527 }
528 // 3D
529 template<unsigned Dim, class MFLOAT>
530 void UniformCartesian<Dim,MFLOAT>::
531 storeSpacingFields(e_dim_tag p1, e_dim_tag p2, e_dim_tag p3, int vnodes)
532 {
533  e_dim_tag et[3];
534  et[0] = p1;
535  et[1] = p2;
536  et[2] = p3;
537  storeSpacingFields(et, vnodes);
538 }
539 // The general storeSpacingfields() function; others invoke this internally:
540 template<unsigned Dim, class MFLOAT>
541 void UniformCartesian<Dim,MFLOAT>::
542 storeSpacingFields(e_dim_tag* et, int vnodes)
543 {
544  // VERTEX-VERTEX SPACINGS (same as CELL-CELL SPACINGS for uniform):
545  NDIndex<Dim> cells, verts;
546  unsigned int d;
547  for (d=0; d<Dim; d++) {
548  cells[d] = Index(gridSizes[d]-1);
549  verts[d] = Index(gridSizes[d]);
550  }
551  if (!hasSpacingFields) {
552  // allocate layout and spacing field
553  FlCell = new FieldLayout<Dim>(cells, et, vnodes);
554  // Note: enough guard cells only for existing Div(), etc. implementations:
555  // (not really used by Div() etc for UniformCartesian); someday should make
556  // this user-settable.
557  VertSpacings =
558  new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlCell,GuardCellSizes<Dim>(1));
559  // Added 12/8/98 --TJW:
560  FlVert =
561  new FieldLayout<Dim>(verts, et, vnodes);
562  // Note: enough guard cells only for existing Div(), etc. implementations:
563  CellSpacings =
564  new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlVert,GuardCellSizes<Dim>(1));
565  }
566  BareField<Vektor<MFLOAT,Dim>,Dim>& vertSpacings = *VertSpacings;
567  Vektor<MFLOAT,Dim> vertexSpacing;
568  for (d=0; d<Dim; d++)
569  vertexSpacing[d] = meshSpacing[d];
570  vertSpacings = vertexSpacing;
571  //-------------------------------------------------
572  // Now the hard part, filling in the guard cells:
573  //-------------------------------------------------
574  // The easy part of the hard part is filling so that all the internal
575  // guard layers are right:
576  vertSpacings.fillGuardCells();
577  // The hard part of the hard part is filling the external guard layers,
578  // using the mesh BC to figure out how:
579  // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
580  // Temporaries used in loop over faces
581  Vektor<MFLOAT,Dim> v0,v1; v0 = 0.0; v1 = 1.0; // Used for Reflective mesh BC
582  unsigned int face;
583  typedef Vektor<MFLOAT,Dim> T; // Used multipple places in loop below
584  typename BareField<T,Dim>::iterator_if vfill_i; // Iterator used below
585  int voffset; // Pointer offsets used with LField::iterator below
586  // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
587  for (face=0; face < 2*Dim; face++) {
588  // NDIndex's spanning elements and guard elements:
589  NDIndex<Dim> vSlab = AddGuardCells(cells,vertSpacings.getGuardCellSizes());
590  // Shrink it down to be the guards along the active face:
591  d = face/2;
592  // The following bitwise AND logical test returns true if face is odd
593  // (meaning the "high" or "right" face in the numbering convention) and
594  // returns false if face is even (meaning the "low" or "left" face in
595  // the numbering convention):
596  if ( face & 1u ) {
597  vSlab[d] = Index(cells[d].max() + 1,
598  cells[d].max() + vertSpacings.rightGuard(d));
599  } else {
600  vSlab[d] = Index(cells[d].min() - vertSpacings.leftGuard(d),
601  cells[d].min() - 1);
602  }
603  // Compute pointer offsets used with LField::iterator below:
604  // Treat all as Reflective BC (see Cartesian for comparison); for
605  // uniform cartesian mesh, all mesh BC's equivalent for this purpose:
606  if ( face & 1u ) {
607  voffset = 2*cells[d].max() + 1 - 1;
608  } else {
609  voffset = 2*cells[d].min() - 1 + 1;
610  }
611 
612  // +++++++++++++++vertSpacings++++++++++++++
613  for (vfill_i=vertSpacings.begin_if();
614  vfill_i!=vertSpacings.end_if(); ++vfill_i)
615  {
616  // Cache some things we will use often below.
617  // Pointer to the data for the current LField (right????):
618  LField<T,Dim> &fill = *(*vfill_i).second;
619  // NDIndex spanning all elements in the LField, including the guards:
620  const NDIndex<Dim> &fill_alloc = fill.getAllocated();
621  // If the previously-created boundary guard-layer NDIndex "cSlab"
622  // contains any of the elements in this LField (they will be guard
623  // elements if it does), assign the values into them here by applying
624  // the boundary condition:
625  if ( vSlab.touches( fill_alloc ) )
626  {
627  // Find what it touches in this LField.
628  NDIndex<Dim> dest = vSlab.intersect( fill_alloc );
629 
630  // For exrapolation boundary conditions, the boundary guard-layer
631  // elements are typically copied from interior values; the "src"
632  // NDIndex specifies the interior elements to be copied into the
633  // "dest" boundary guard-layer elements (possibly after some
634  // mathematical operations like multipplying by minus 1 later):
635  NDIndex<Dim> src = dest; // Create dest equal to src
636  // Now calculate the interior elements; the voffset variable
637  // computed above makes this right for "low" or "high" face cases:
638  src[d] = voffset - src[d];
639 
640  // TJW: Why is there another loop over LField's here??????????
641  // Loop over the ones that src touches.
642  typename BareField<T,Dim>::iterator_if from_i;
643  for (from_i=vertSpacings.begin_if();
644  from_i!=vertSpacings.end_if(); ++from_i)
645  {
646  // Cache a few things.
647  LField<T,Dim> &from = *(*from_i).second;
648  const NDIndex<Dim> &from_owned = from.getOwned();
649  const NDIndex<Dim> &from_alloc = from.getAllocated();
650  // If src touches this LField...
651  if ( src.touches( from_owned ) )
652  {
653  NDIndex<Dim> from_it = src.intersect( from_alloc );
654  NDIndex<Dim> vfill_it = dest.plugBase( from_it );
655  // Build iterators for the copy...
656  typedef typename LField<T,Dim>::iterator LFI;
657  LFI lhs = fill.begin(vfill_it);
658  LFI rhs = from.begin(from_it);
659  // And do the assignment (reflective BC hardwired):
660  BrickExpression<Dim,LFI,LFI,OpUMeshExtrapolate<T> >
661  (lhs,rhs,OpUMeshExtrapolate<T>(v0,v1)).apply();
662  }
663  }
664  }
665  }
666 
667  }
668 
669  // For uniform cartesian mesh, cell-cell spacings are identical to
670  // vert-vert spacings:
671  //12/8/98 CellSpacings = VertSpacings;
672  // Added 12/8/98 --TJW:
673  BareField<Vektor<MFLOAT,Dim>,Dim>& cellSpacings = *CellSpacings;
674  cellSpacings = vertexSpacing;
675 
676  hasSpacingFields = true; // Flag this as having been done to this object.
677 }
678 
679 // These specify both the total number of vnodes and the numbers of vnodes
680 // along each dimension for the partitioning of the index space. Obviously
681 // this restricts the number of vnodes to be a product of the numbers along
682 // each dimension (the constructor implementation checks this): Special
683 // cases for 1-3 dimensions, ala FieldLayout ctors (see FieldLayout.h for
684 // more relevant comments, including definition of recurse):
685 
686 // 1D
687 template<unsigned Dim, class MFLOAT>
688 void UniformCartesian<Dim,MFLOAT>::
689 storeSpacingFields(e_dim_tag p1,
690  unsigned vnodes1,
691  bool recurse,
692  int vnodes) {
693  e_dim_tag et[1];
694  et[0] = p1;
695  unsigned vnodesPerDirection[Dim];
696  vnodesPerDirection[0] = vnodes1;
697  storeSpacingFields(et, vnodesPerDirection, recurse, vnodes);
698 }
699 // 2D
700 template<unsigned Dim, class MFLOAT>
701 void UniformCartesian<Dim,MFLOAT>::
702 storeSpacingFields(e_dim_tag p1, e_dim_tag p2,
703  unsigned vnodes1, unsigned vnodes2,
704  bool recurse,int vnodes) {
705  e_dim_tag et[2];
706  et[0] = p1;
707  et[1] = p2;
708  unsigned vnodesPerDirection[Dim];
709  vnodesPerDirection[0] = vnodes1;
710  vnodesPerDirection[1] = vnodes2;
711  storeSpacingFields(et, vnodesPerDirection, recurse, vnodes);
712 }
713 // 3D
714 template<unsigned Dim, class MFLOAT>
715 void UniformCartesian<Dim,MFLOAT>::
716 storeSpacingFields(e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
717  unsigned vnodes1, unsigned vnodes2, unsigned vnodes3,
718  bool recurse, int vnodes) {
719  e_dim_tag et[3];
720  et[0] = p1;
721  et[1] = p2;
722  et[2] = p3;
723  unsigned vnodesPerDirection[Dim];
724  vnodesPerDirection[0] = vnodes1;
725  vnodesPerDirection[1] = vnodes2;
726  vnodesPerDirection[2] = vnodes3;
727  storeSpacingFields(et, vnodesPerDirection, recurse, vnodes);
728 }
729 
730 // TJW: Note: should clean up here eventually, and put redundant code from
731 // this and the other general storeSpacingFields() implementation into one
732 // function. Need to check this in quickly for Blanca right now --12/8/98
733 // The general storeSpacingfields() function; others invoke this internally:
734 template<unsigned Dim, class MFLOAT>
735 void UniformCartesian<Dim,MFLOAT>::
736 storeSpacingFields(e_dim_tag *p,
737  unsigned* vnodesPerDirection,
738  bool recurse, int vnodes)
739 {
740  // VERTEX-VERTEX SPACINGS (same as CELL-CELL SPACINGS for uniform):
741  NDIndex<Dim> cells, verts;
742  unsigned int d;
743  for (d=0; d<Dim; d++) {
744  cells[d] = Index(gridSizes[d]-1);
745  verts[d] = Index(gridSizes[d]);
746  }
747  if (!hasSpacingFields) {
748  // allocate layout and spacing field
749  FlCell =
750  new FieldLayout<Dim>(cells, p, vnodesPerDirection, recurse, vnodes);
751  // Note: enough guard cells only for existing Div(), etc. implementations:
752  // (not really used by Div() etc for UniformCartesian); someday should make
753  // this user-settable.
754  VertSpacings =
755  new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlCell,GuardCellSizes<Dim>(1));
756  // Added 12/8/98 --TJW:
757  FlVert =
758  new FieldLayout<Dim>(verts, p, vnodesPerDirection, recurse, vnodes);
759  // Note: enough guard cells only for existing Div(), etc. implementations:
760  CellSpacings =
761  new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlVert,GuardCellSizes<Dim>(1));
762  }
763  BareField<Vektor<MFLOAT,Dim>,Dim>& vertSpacings = *VertSpacings;
764  Vektor<MFLOAT,Dim> vertexSpacing;
765  for (d=0; d<Dim; d++)
766  vertexSpacing[d] = meshSpacing[d];
767  vertSpacings = vertexSpacing;
768  //-------------------------------------------------
769  // Now the hard part, filling in the guard cells:
770  //-------------------------------------------------
771  // The easy part of the hard part is filling so that all the internal
772  // guard layers are right:
773  vertSpacings.fillGuardCells();
774  // The hard part of the hard part is filling the external guard layers,
775  // using the mesh BC to figure out how:
776  // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
777  // Temporaries used in loop over faces
778  Vektor<MFLOAT,Dim> v0,v1; v0 = 0.0; v1 = 1.0; // Used for Reflective mesh BC
779  unsigned int face;
780  typedef Vektor<MFLOAT,Dim> T; // Used multipple places in loop below
781  typename BareField<T,Dim>::iterator_if vfill_i; // Iterator used below
782  int voffset; // Pointer offsets used with LField::iterator below
783  // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
784  for (face=0; face < 2*Dim; face++) {
785  // NDIndex's spanning elements and guard elements:
786  NDIndex<Dim> vSlab = AddGuardCells(cells,vertSpacings.getGuardCellSizes());
787  // Shrink it down to be the guards along the active face:
788  d = face/2;
789  // The following bitwise AND logical test returns true if face is odd
790  // (meaning the "high" or "right" face in the numbering convention) and
791  // returns false if face is even (meaning the "low" or "left" face in
792  // the numbering convention):
793  if ( face & 1u ) {
794  vSlab[d] = Index(cells[d].max() + 1,
795  cells[d].max() + vertSpacings.rightGuard(d));
796  } else {
797  vSlab[d] = Index(cells[d].min() - vertSpacings.leftGuard(d),
798  cells[d].min() - 1);
799  }
800  // Compute pointer offsets used with LField::iterator below:
801  // Treat all as Reflective BC (see Cartesian for comparison); for
802  // uniform cartesian mesh, all mesh BC's equivalent for this purpose:
803  if ( face & 1u ) {
804  voffset = 2*cells[d].max() + 1 - 1;
805  } else {
806  voffset = 2*cells[d].min() - 1 + 1;
807  }
808 
809  // +++++++++++++++vertSpacings++++++++++++++
810  for (vfill_i=vertSpacings.begin_if();
811  vfill_i!=vertSpacings.end_if(); ++vfill_i)
812  {
813  // Cache some things we will use often below.
814  // Pointer to the data for the current LField (right????):
815  LField<T,Dim> &fill = *(*vfill_i).second;
816  // NDIndex spanning all elements in the LField, including the guards:
817  const NDIndex<Dim> &fill_alloc = fill.getAllocated();
818  // If the previously-created boundary guard-layer NDIndex "cSlab"
819  // contains any of the elements in this LField (they will be guard
820  // elements if it does), assign the values into them here by applying
821  // the boundary condition:
822  if ( vSlab.touches( fill_alloc ) )
823  {
824  // Find what it touches in this LField.
825  NDIndex<Dim> dest = vSlab.intersect( fill_alloc );
826 
827  // For exrapolation boundary conditions, the boundary guard-layer
828  // elements are typically copied from interior values; the "src"
829  // NDIndex specifies the interior elements to be copied into the
830  // "dest" boundary guard-layer elements (possibly after some
831  // mathematical operations like multipplying by minus 1 later):
832  NDIndex<Dim> src = dest; // Create dest equal to src
833  // Now calculate the interior elements; the voffset variable
834  // computed above makes this right for "low" or "high" face cases:
835  src[d] = voffset - src[d];
836 
837  // TJW: Why is there another loop over LField's here??????????
838  // Loop over the ones that src touches.
839  typename BareField<T,Dim>::iterator_if from_i;
840  for (from_i=vertSpacings.begin_if();
841  from_i!=vertSpacings.end_if(); ++from_i)
842  {
843  // Cache a few things.
844  LField<T,Dim> &from = *(*from_i).second;
845  const NDIndex<Dim> &from_owned = from.getOwned();
846  const NDIndex<Dim> &from_alloc = from.getAllocated();
847  // If src touches this LField...
848  if ( src.touches( from_owned ) )
849  {
850  NDIndex<Dim> from_it = src.intersect( from_alloc );
851  NDIndex<Dim> vfill_it = dest.plugBase( from_it );
852  // Build iterators for the copy...
853  typedef typename LField<T,Dim>::iterator LFI;
854  LFI lhs = fill.begin(vfill_it);
855  LFI rhs = from.begin(from_it);
856  // And do the assignment (reflective BC hardwired):
857  BrickExpression<Dim,LFI,LFI,OpUMeshExtrapolate<T> >
858  (lhs,rhs,OpUMeshExtrapolate<T>(v0,v1)).apply();
859  }
860  }
861  }
862  }
863 
864  }
865 
866  // For uniform cartesian mesh, cell-cell spacings are identical to
867  // vert-vert spacings:
868  //12/8/98 CellSpacings = VertSpacings;
869  // Added 12/8/98 --TJW:
870  BareField<Vektor<MFLOAT,Dim>,Dim>& cellSpacings = *CellSpacings;
871  cellSpacings = vertexSpacing;
872 
873  hasSpacingFields = true; // Flag this as having been done to this object.
874 }
875 
876 //-----------------------------------------------------------------------------
877 // I/O:
878 //-----------------------------------------------------------------------------
879 // Formatted output of UniformCartesian object:
880 template< unsigned Dim, class MFLOAT >
881 void
882 UniformCartesian<Dim,MFLOAT>::
883 print(std::ostream& out)
884 {
885  Inform info("", out);
886  print(info);
887 }
888 
889 template< unsigned Dim, class MFLOAT >
890 void
891 UniformCartesian<Dim,MFLOAT>::
892 print(Inform& out)
893 {
894  out << "======UniformCartesian<" << Dim << ",MFLOAT>==begin======\n";
895  unsigned int d;
896  for (d=0; d < Dim; d++)
897  out << "gridSizes[" << d << "] = " << gridSizes[d] << "\n";
898  out << "origin = " << origin << "\n";
899  for (d=0; d < Dim; d++)
900  out << "meshSpacing[" << d << "] = " << meshSpacing[d] << "\n";
901  for (d=0; d < (1u<<Dim); d++)
902  out << "Dvc[" << d << "] = " << Dvc[d] << "\n";
903  out << "cell volume = " << volume << "\n";
904  out << "======UniformCartesian<" << Dim << ",MFLOAT>==end========\n";
905 }
906 
907 //--------------------------------------------------------------------------
908 // Various (UniformCartesian) mesh mechanisms:
909 //--------------------------------------------------------------------------
910 
911 // Volume of single cell indexed by NDIndex:
912 template <unsigned Dim, class MFLOAT>
913 MFLOAT
914 UniformCartesian<Dim,MFLOAT>::
915 getCellVolume(const NDIndex<Dim>& ndi) const
916 {
917  unsigned int d;
918  for (d=0; d<Dim; d++)
919  if (ndi[d].length() != 1)
920  ERRORMSG("UniformCartesian::getCellVolume() error: arg is not a NDIndex"
921  << "specifying a single element" << endl);
922  return volume;
923 }
924 // Field of volumes of all cells:
925 template <unsigned Dim, class MFLOAT>
926 Field<MFLOAT,Dim,UniformCartesian<Dim,MFLOAT>,Cell>&
927 UniformCartesian<Dim,MFLOAT>::
928 getCellVolumeField(Field<MFLOAT,Dim,UniformCartesian<Dim,MFLOAT>,Cell>&
929  volumes) const
930 {
931  volumes = volume;
932  return volumes;
933 }
934 // Volume of range of cells bounded by verticies specified by input NDIndex;
935 template <unsigned Dim, class MFLOAT>
936 MFLOAT
937 UniformCartesian<Dim,MFLOAT>::
938 getVertRangeVolume(const NDIndex<Dim>& ndi) const
939 {
940  // Get vertex positions of extremal cells:
941  Vektor<MFLOAT,Dim> v0, v1;
942  NDIndex<Dim> ndi0, ndi1;
943  unsigned int d;
944  int i0, i1;
945  for (d=0; d<Dim; d++) {
946  i0 = ndi[d].first();
947  i1 = ndi[d].last();
948  ndi0[d] = Index(i0,i0,1); // Bounding vertex (from below)
949  ndi1[d] = Index(i1,i1,1); // Bounding vertex (from above)
950  }
951  v0 = getVertexPosition(ndi0);
952  v1 = getVertexPosition(ndi1);
953  // Compute volume of rectangular solid beweeen these extremal vertices:
954  MFLOAT volume = 1.0;
955  for (d=0; d<Dim; d++) volume *= abs(v1(d) - v0(d));
956  return volume;
957 }
958 // Volume of range of cells spanned by input NDIndex (index of cells):
959 template <unsigned Dim, class MFLOAT>
960 MFLOAT
961 UniformCartesian<Dim,MFLOAT>::
962 getCellRangeVolume(const NDIndex<Dim>& ndi) const
963 {
964  // Get vertex positions bounding extremal cells:
965  Vektor<MFLOAT,Dim> v0, v1;
966  NDIndex<Dim> ndi0, ndi1;
967  unsigned int d;
968  int i0, i1;
969  for (d=0; d<Dim; d++) {
970  i0 = ndi[d].first();
971  i1 = ndi[d].last() + 1;
972  ndi0[d] = Index(i0,i0,1); // Bounding vertex (from below)
973  ndi1[d] = Index(i1,i1,1); // Bounding vertex (from above)
974  }
975  v0 = getVertexPosition(ndi0);
976  v1 = getVertexPosition(ndi1);
977  // Compute volume of rectangular solid beweeen these extremal vertices:
978  MFLOAT volume = 1.0;
979  for (d=0; d<Dim; d++) volume *= abs(v1(d) - v0(d));
980  return volume;
981 }
982 
983 // Nearest vertex index to (x,y,z):
984 template <unsigned Dim, class MFLOAT>
985 NDIndex<Dim>
986 UniformCartesian<Dim,MFLOAT>::
987 getNearestVertex(const Vektor<MFLOAT,Dim>& x) const
988 {
989  // Find coordinate vectors of the vertices just above and just below the
990  // input point (extremal vertices on cell containing point):
991  NDIndex<Dim> ndi;
992  int i;
993  for (unsigned int d=0; d<Dim; d++) {
994  i = (int)((x(d) - origin(d))/meshSpacing[d] + 0.5);
995  if (x(d) >= origin(d))
996  ndi[d] = Index(i,i);
997  else
998  ndi[d] = Index(i-1,i-1);
999  }
1000  return ndi;
1001 }
1002 // Nearest vertex index with all vertex coordinates below (x,y,z):
1003 template <unsigned Dim, class MFLOAT>
1004 NDIndex<Dim>
1005 UniformCartesian<Dim,MFLOAT>::
1006 getVertexBelow(const Vektor<MFLOAT,Dim>& x) const
1007 {
1008  // Find coordinate vectors of the vertices just above and just below the
1009  // input point (extremal vertices on cell containing point):
1010  NDIndex<Dim> ndi;
1011  int i;
1012  for (unsigned int d=0; d<Dim; d++) {
1013  i = (int)((x(d) - origin(d))/meshSpacing[d]);
1014  if (x(d) >= origin(d))
1015  ndi[d] = Index(i,i);
1016  else
1017  ndi[d] = Index(i-1,i-1);
1018  }
1019  return ndi;
1020 }
1021 // (x,y,z) coordinates of indexed vertex:
1022 template <unsigned Dim, class MFLOAT>
1023 Vektor<MFLOAT,Dim>
1024 UniformCartesian<Dim,MFLOAT>::
1025 getVertexPosition(const NDIndex<Dim>& ndi) const
1026 {
1027  unsigned int d;
1028  for (d=0; d<Dim; d++) {
1029  if (ndi[d].length() != 1)
1030  ERRORMSG("UniformCartesian::getVertexPosition() error: arg is not a"
1031  << " NDIndex specifying a single element" << endl);
1032  }
1033  // N.B.: following may need modification to be right for periodic Mesh BC:
1034  Vektor<MFLOAT,Dim> vertexPosition;
1035  for (d=0; d<Dim; d++)
1036  vertexPosition(d) = ndi[d].first()*meshSpacing[d] + origin(d);
1037  return vertexPosition;
1038 }
1039 // Field of (x,y,z) coordinates of all vertices:
1040 template <unsigned Dim, class MFLOAT>
1041 Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Vert>&
1042 UniformCartesian<Dim,MFLOAT>::
1043 getVertexPositionField(Field<Vektor<MFLOAT,Dim>,Dim,
1044  UniformCartesian<Dim,MFLOAT>,Vert>& vertexPositions) const
1045 {
1046  unsigned int d;
1047  int currentLocation[Dim];
1048  Vektor<MFLOAT,Dim> vertexPosition;
1049  vertexPositions.Uncompress(); // uncompress field before entering values!
1050  typename Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Vert>::iterator fi,
1051  fi_end = vertexPositions.end();
1052  for (fi = vertexPositions.begin(); fi != fi_end; ++fi) {
1053  fi.GetCurrentLocation(currentLocation);
1054  for (d=0; d<Dim; d++)
1055  vertexPosition(d) = origin(d) + currentLocation[d]*meshSpacing[d];
1056  *fi = vertexPosition;
1057  }
1058  return vertexPositions;
1059 }
1060 
1061 // (x,y,z) coordinates of indexed cell:
1062 template <unsigned Dim, class MFLOAT>
1063 Vektor<MFLOAT,Dim>
1064 UniformCartesian<Dim,MFLOAT>::
1065 getCellPosition(const NDIndex<Dim>& ndi) const
1066 {
1067  unsigned int d;
1068  for (d=0; d<Dim; d++) {
1069  if (ndi[d].length() != 1)
1070  ERRORMSG("UniformCartesian::getCellPosition() error: arg is not a"
1071  << " NDIndex specifying a single element" << endl);
1072  }
1073  // N.B.: following may need modification to be right for periodic Mesh BC:
1074  Vektor<MFLOAT,Dim> cellPosition;
1075  for (d=0; d<Dim; d++)
1076  cellPosition(d) = (ndi[d].first()+0.5)*meshSpacing[d] + origin(d);
1077  return cellPosition;
1078 }
1079 // Field of (x,y,z) coordinates of all cells:
1080 template <unsigned Dim, class MFLOAT>
1081 Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Cell>&
1082 UniformCartesian<Dim,MFLOAT>::
1083 getCellPositionField(Field<Vektor<MFLOAT,Dim>,Dim,
1084  UniformCartesian<Dim,MFLOAT>,Cell>& cellPositions) const
1085 {
1086  unsigned int d;
1087  int currentLocation[Dim];
1088  Vektor<MFLOAT,Dim> cellPosition;
1089  cellPositions.Uncompress(); // uncompress field before entering values!
1090  typename Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Cell>::iterator fi,
1091  fi_end = cellPositions.end();
1092  for (fi = cellPositions.begin(); fi != fi_end; ++fi) {
1093  fi.GetCurrentLocation(currentLocation);
1094  for (d=0; d<Dim; d++)
1095  cellPosition(d) = origin(d) + (currentLocation[d]+0.5)*meshSpacing[d];
1096  *fi = cellPosition;
1097  }
1098  return cellPositions;
1099 }
1100 
1101 // Vertex-vertex grid spacing of indexed cell:
1102 template <unsigned Dim, class MFLOAT>
1103 Vektor<MFLOAT,Dim>
1104 UniformCartesian<Dim,MFLOAT>::
1105 getDeltaVertex(const NDIndex<Dim>& ndi) const
1106 {
1107  // bfh: updated to compute interval for a whole index range
1108  Vektor<MFLOAT,Dim> vertexVertexSpacing;
1109  for (unsigned int d=0; d<Dim; d++)
1110  vertexVertexSpacing[d] = meshSpacing[d] * ndi[d].length();
1111  return vertexVertexSpacing;
1112 }
1113 
1114 // Field of vertex-vertex grid spacings of all cells:
1115 template <unsigned Dim, class MFLOAT>
1116 Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Cell>&
1117 UniformCartesian<Dim,MFLOAT>::
1118 getDeltaVertexField(Field<Vektor<MFLOAT,Dim>,Dim,
1119  UniformCartesian<Dim,MFLOAT>,Cell>& vertexSpacings) const
1120 {
1121  Vektor<MFLOAT,Dim> vertexVertexSpacing;
1122  unsigned int d;
1123  for (d=0; d<Dim; d++) vertexVertexSpacing(d) = meshSpacing[d];
1124  vertexSpacings = vertexVertexSpacing;
1125  return vertexSpacings;
1126 }
1127 
1128 // Cell-cell grid spacing of indexed vertex:
1129 template <unsigned Dim, class MFLOAT>
1130 Vektor<MFLOAT,Dim>
1131 UniformCartesian<Dim,MFLOAT>::
1132 getDeltaCell(const NDIndex<Dim>& ndi) const
1133 {
1134  // bfh: updated to compute interval for a whole index range
1135  Vektor<MFLOAT,Dim> cellCellSpacing;
1136  for (unsigned int d=0; d<Dim; d++)
1137  cellCellSpacing[d] = meshSpacing[d] * ndi[d].length();
1138  return cellCellSpacing;
1139 }
1140 
1141 // Field of cell-cell grid spacings of all vertices:
1142 template <unsigned Dim, class MFLOAT>
1143 Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Vert>&
1144 UniformCartesian<Dim,MFLOAT>::
1145 getDeltaCellField(Field<Vektor<MFLOAT,Dim>,Dim,
1146  UniformCartesian<Dim,MFLOAT>,Vert>& cellSpacings) const
1147 {
1148  Vektor<MFLOAT,Dim> cellCellSpacing;
1149  unsigned int d;
1150  for (d=0; d<Dim; d++) cellCellSpacing(d) = meshSpacing[d];
1151  cellSpacings = cellCellSpacing;
1152  return cellSpacings;
1153 }
1154 // Array of surface normals to cells adjoining indexed cell:
1155 template <unsigned Dim, class MFLOAT>
1156 Vektor<MFLOAT,Dim>*
1157 UniformCartesian<Dim,MFLOAT>::
1158 getSurfaceNormals(const NDIndex<Dim>& /*ndi*/) const
1159 {
1160  Vektor<MFLOAT,Dim>* surfaceNormals = new Vektor<MFLOAT,Dim>[2*Dim];
1161  unsigned int d, i;
1162  for (d=0; d<Dim; d++) {
1163  for (i=0; i<Dim; i++) {
1164  surfaceNormals[2*d](i) = 0.0;
1165  surfaceNormals[2*d+1](i) = 0.0;
1166  }
1167  surfaceNormals[2*d](d) = -1.0;
1168  surfaceNormals[2*d+1](d) = 1.0;
1169  }
1170  return surfaceNormals;
1171 }
1172 // Array of (pointers to) Fields of surface normals to all cells:
1173 template <unsigned Dim, class MFLOAT>
1174 void
1175 UniformCartesian<Dim,MFLOAT>::
1176 getSurfaceNormalFields(Field<Vektor<MFLOAT,Dim>, Dim,
1177  UniformCartesian<Dim,MFLOAT>,Cell>**
1178  surfaceNormalsFields ) const
1179 {
1180  Vektor<MFLOAT,Dim>* surfaceNormals = new Vektor<MFLOAT,Dim>[2*Dim];
1181  unsigned int d, i;
1182  for (d=0; d<Dim; d++) {
1183  for (i=0; i<Dim; i++) {
1184  surfaceNormals[2*d](i) = 0.0;
1185  surfaceNormals[2*d+1](i) = 0.0;
1186  }
1187  surfaceNormals[2*d](d) = -1.0;
1188  surfaceNormals[2*d+1](d) = 1.0;
1189  }
1190  for (d=0; d<2*Dim; d++) assign((*(surfaceNormalsFields[d])),
1191  surfaceNormals[d]);
1192  // return surfaceNormalsFields;
1193 }
1194 // Similar functions, but specify the surface normal to a single face, using
1195 // the following numbering convention: 0 means low face of 1st dim, 1 means
1196 // high face of 1st dim, 2 means low face of 2nd dim, 3 means high face of
1197 // 2nd dim, and so on:
1198 // Surface normal to face on indexed cell:
1199 template <unsigned Dim, class MFLOAT>
1200 Vektor<MFLOAT,Dim>
1201 UniformCartesian<Dim,MFLOAT>::
1202 getSurfaceNormal(const NDIndex<Dim>& /*ndi*/, unsigned face) const
1203 {
1204  Vektor<MFLOAT,Dim> surfaceNormal;
1205  unsigned int d;
1206  // The following bitwise AND logical test returns true if face is odd
1207  // (meaning the "high" or "right" face in the numbering convention) and
1208  // returns false if face is even (meaning the "low" or "left" face in the
1209  // numbering convention):
1210  if ( face & 1 ) {
1211  for (d=0; d<Dim; d++) {
1212  if ((face/2) == d) {
1213  surfaceNormal(face) = -1.0;
1214  } else {
1215  surfaceNormal(face) = 0.0;
1216  }
1217  }
1218  } else {
1219  for (d=0; d<Dim; d++) {
1220  if ((face/2) == d) {
1221  surfaceNormal(face) = 1.0;
1222  } else {
1223  surfaceNormal(face) = 0.0;
1224  }
1225  }
1226  }
1227  return surfaceNormal;
1228 }
1229 // Field of surface normals to face on all cells:
1230 template <unsigned Dim, class MFLOAT>
1231 Field<Vektor<MFLOAT,Dim>,Dim,UniformCartesian<Dim,MFLOAT>,Cell>&
1232 UniformCartesian<Dim,MFLOAT>::
1233 getSurfaceNormalField(Field<Vektor<MFLOAT,Dim>, Dim,
1234  UniformCartesian<Dim,MFLOAT>,Cell>& surfaceNormalField,
1235  unsigned face) const
1236 {
1237  Vektor<MFLOAT,Dim> surfaceNormal;
1238  unsigned int d;
1239  // The following bitwise AND logical test returns true if face is odd
1240  // (meaning the "high" or "right" face in the numbering convention) and
1241  // returns false if face is even (meaning the "low" or "left" face in the
1242  // numbering convention):
1243  if ( face & 1 ) {
1244  for (d=0; d<Dim; d++) {
1245  if ((face/2) == d) {
1246  surfaceNormal(face) = -1.0;
1247  } else {
1248  surfaceNormal(face) = 0.0;
1249  }
1250  }
1251  } else {
1252  for (d=0; d<Dim; d++) {
1253  if ((face/2) == d) {
1254  surfaceNormal(face) = 1.0;
1255  } else {
1256  surfaceNormal(face) = 0.0;
1257  }
1258  }
1259  }
1260  surfaceNormalField = surfaceNormal;
1261  return surfaceNormalField;
1262 }
1263 
1264 
1265 //--------------------------------------------------------------------------
1266 // Global functions
1267 //--------------------------------------------------------------------------
1268 
1269 
1270 //*****************************************************************************
1271 // Stuff taken from old Cartesian.h, which now applies to UniformCartesian:
1272 //*****************************************************************************
1273 
1274 //----------------------------------------------------------------------
1275 // Divergence Vektor/Vert -> Scalar/Cell
1276 //----------------------------------------------------------------------
1277 template < class T, class MFLOAT >
1278 Field<T,1U,UniformCartesian<1U,MFLOAT>,Cell>&
1279 Div(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
1280  Field<T,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
1281 {
1282  const NDIndex<1U>& domain = r.getDomain();
1283  Index I = domain[0];
1284 
1285  assign( r[I] ,
1286  dot(x[I ], x.get_mesh().Dvc[0]) +
1287  dot(x[I+1], x.get_mesh().Dvc[1]));
1288  return r;
1289 }
1290 //----------------------------------------------------------------------
1291 template < class T, class MFLOAT >
1292 Field<T,2U,UniformCartesian<2U,MFLOAT>,Cell>&
1293 Div(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
1294  Field<T,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
1295 {
1296  const NDIndex<2U>& domain = r.getDomain();
1297  Index I = domain[0];
1298  Index J = domain[1];
1299 
1300  assign( r[I][J] ,
1301  dot(x[I ][J ], x.get_mesh().Dvc[0]) +
1302  dot(x[I+1][J ], x.get_mesh().Dvc[1]) +
1303  dot(x[I ][J+1], x.get_mesh().Dvc[2]) +
1304  dot(x[I+1][J+1], x.get_mesh().Dvc[3]));
1305  return r;
1306 }
1307 //----------------------------------------------------------------------
1308 template < class T, class MFLOAT >
1309 Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>&
1310 Div(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
1311  Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
1312 {
1313  const NDIndex<3U>& domain = r.getDomain();
1314  Index I = domain[0];
1315  Index J = domain[1];
1316  Index K = domain[2];
1317 
1318  assign( r[I][J][K] ,
1319  dot(x[I ][J ][K ], x.get_mesh().Dvc[0]) +
1320  dot(x[I+1][J ][K ], x.get_mesh().Dvc[1]) +
1321  dot(x[I ][J+1][K ], x.get_mesh().Dvc[2]) +
1322  dot(x[I+1][J+1][K ], x.get_mesh().Dvc[3]) +
1323  dot(x[I ][J ][K+1], x.get_mesh().Dvc[4]) +
1324  dot(x[I+1][J ][K+1], x.get_mesh().Dvc[5]) +
1325  dot(x[I ][J+1][K+1], x.get_mesh().Dvc[6]) +
1326  dot(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]));
1327  return r;
1328 }
1329 //----------------------------------------------------------------------
1330 // Divergence Vektor/Cell -> Scalar/Vert
1331 //----------------------------------------------------------------------
1332 template < class T, class MFLOAT >
1333 Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>&
1334 Div(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
1335  Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
1336 {
1337  const NDIndex<1U>& domain = r.getDomain();
1338  Index I = domain[0];
1339 
1340  assign( r[I] ,
1341  dot(x[I-1], x.get_mesh().Dvc[0]) +
1342  dot(x[I ], x.get_mesh().Dvc[1]));
1343  return r;
1344 }
1345 //----------------------------------------------------------------------
1346 // (tjw: the one in cv.cpp)
1347 template < class T, class MFLOAT >
1348 Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>&
1349 Div(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
1350  Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
1351 {
1352  const NDIndex<2U>& domain = r.getDomain();
1353  Index I = domain[0];
1354  Index J = domain[1];
1355 
1356  assign( r[I][J] ,
1357  dot(x[I-1][J-1], x.get_mesh().Dvc[0]) +
1358  dot(x[I ][J-1], x.get_mesh().Dvc[1]) +
1359  dot(x[I-1][J ], x.get_mesh().Dvc[2]) +
1360  dot(x[I ][J ], x.get_mesh().Dvc[3]));
1361  return r;
1362 }
1363 //----------------------------------------------------------------------
1364 template < class T, class MFLOAT >
1365 Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>&
1366 Div(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
1367  Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
1368 {
1369  const NDIndex<3U>& domain = r.getDomain();
1370  Index I = domain[0];
1371  Index J = domain[1];
1372  Index K = domain[2];
1373 
1374  assign( r[I][J][K] ,
1375  dot(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]) +
1376  dot(x[I ][J-1][K-1], x.get_mesh().Dvc[1]) +
1377  dot(x[I-1][J ][K-1], x.get_mesh().Dvc[2]) +
1378  dot(x[I ][J ][K-1], x.get_mesh().Dvc[3]) +
1379  dot(x[I-1][J-1][K ], x.get_mesh().Dvc[4]) +
1380  dot(x[I ][J-1][K ], x.get_mesh().Dvc[5]) +
1381  dot(x[I-1][J ][K ], x.get_mesh().Dvc[6]) +
1382  dot(x[I ][J ][K ], x.get_mesh().Dvc[7]));
1383  return r;
1384 }
1385 //----------------------------------------------------------------------
1386 // Divergence Vektor/Vert -> Scalar/Vert
1387 //----------------------------------------------------------------------
1388 template < class T, class MFLOAT >
1389 Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>&
1390 Div(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
1391  Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
1392 {
1393  const NDIndex<1U>& domain = r.getDomain();
1394  Index I = domain[0];
1395  Vektor<T,1U> idx;
1396  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
1397 
1398  assign( r[I] , dot( idx , (x[I+1] - x[I-1]) ) );
1399  return r;
1400 }
1401 //----------------------------------------------------------------------
1402 template < class T, class MFLOAT >
1403 Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>&
1404 Div(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
1405  Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
1406 {
1407  const NDIndex<2U>& domain = r.getDomain();
1408  Index I = domain[0];
1409  Index J = domain[1];
1410  Vektor<T,2U> idx,idy;
1411  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
1412  idx[1] = 0.0;
1413  idy[0] = 0.0;
1414  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
1415 
1416  assign( r[I][J] ,
1417  dot( idx , (x[I+1][J ] - x[I-1][J ])) +
1418  dot( idy , (x[I ][J+1] - x[I ][J-1])) );
1419  return r;
1420 }
1421 //----------------------------------------------------------------------
1422 template < class T, class MFLOAT >
1423 Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>&
1424 Div(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
1425  Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
1426 {
1427  const NDIndex<3U>& domain = r.getDomain();
1428  Index I = domain[0];
1429  Index J = domain[1];
1430  Index K = domain[2];
1431  Vektor<T,3U> idx,idy,idz;
1432  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
1433  idx[1] = 0.0;
1434  idx[2] = 0.0;
1435  idy[0] = 0.0;
1436  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
1437  idy[2] = 0.0;
1438  idz[0] = 0.0;
1439  idz[1] = 0.0;
1440  idz[2] = 0.5/x.get_mesh().get_meshSpacing(2);
1441 
1442  assign( r[I][J][K] ,
1443  dot(idx , (x[I+1][J ][K ] - x[I-1][J ][K ] )) +
1444  dot(idy , (x[I ][J+1][K ] - x[I ][J-1][K ] )) +
1445  dot(idz , (x[I ][J ][K+1] - x[I ][J ][K-1] )) );
1446  return r;
1447 }
1448 //----------------------------------------------------------------------
1449 // Divergence Vektor/Edge -> Scalar/Vert
1450 //----------------------------------------------------------------------
1451 template < class T, class MFLOAT >
1452 Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>&
1453 Div(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Edge>& x,
1454  Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
1455 {
1456  const NDIndex<1U>& domain = r.getDomain();
1457  Index I = domain[0];
1458  Vektor<T,1U> idx;
1459  idx[0] = 1.0/x.get_mesh().get_meshSpacing(0);
1460 
1461  assign( r[I] , dot( idx , (x[I] - x[I-1]) ) );
1462  return r;
1463 }
1464 //----------------------------------------------------------------------
1465 template < class T, class MFLOAT >
1466 Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>&
1467 Div(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Edge>& x,
1468  Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
1469 {
1470  const NDIndex<2U>& domain = r.getDomain();
1471  Index I = domain[0];
1472  Index J = domain[1];
1473  Vektor<T,2U> idx,idy;
1474  idx[0] = 1.0/x.get_mesh().get_meshSpacing(0);
1475  idx[1] = 0.0;
1476  idy[0] = 0.0;
1477  idy[1] = 1.0/x.get_mesh().get_meshSpacing(1);
1478 
1479  assign( r[I][J] ,
1480  dot( idx , (x[I][J] - x[I-1][J ])) +
1481  dot( idy , (x[I][J] - x[I ][J-1])) );
1482  return r;
1483 }
1484 //----------------------------------------------------------------------
1485 template < class T, class MFLOAT >
1486 Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>&
1487 Div(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Edge>& x,
1488  Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
1489 {
1490  const NDIndex<3U>& domain = r.getDomain();
1491  Index I = domain[0];
1492  Index J = domain[1];
1493  Index K = domain[2];
1494  Vektor<T,3U> idx,idy,idz;
1495  idx[0] = 1.0/x.get_mesh().get_meshSpacing(0);
1496  idx[1] = 0.0;
1497  idx[2] = 0.0;
1498  idy[0] = 0.0;
1499  idy[1] = 1.0/x.get_mesh().get_meshSpacing(1);
1500  idy[2] = 0.0;
1501  idz[0] = 0.0;
1502  idz[1] = 0.0;
1503  idz[2] = 1.0/x.get_mesh().get_meshSpacing(2);
1504 
1505  assign( r[I][J][K] ,
1506  dot(idx , (x[I][J][K] - x[I-1][J ][K ] )) +
1507  dot(idy , (x[I][J][K] - x[I ][J-1][K ] )) +
1508  dot(idz , (x[I][J][K] - x[I ][J ][K-1] )) );
1509  return r;
1510 }
1511 //----------------------------------------------------------------------
1512 // Divergence Vektor/Cell -> Scalar/Cell
1513 //----------------------------------------------------------------------
1514 template < class T, class MFLOAT >
1515 Field<T,1U,UniformCartesian<1U,MFLOAT>,Cell>&
1516 Div(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
1517  Field<T,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
1518 {
1519  const NDIndex<1U>& domain = r.getDomain();
1520  Index I = domain[0];
1521  Vektor<T,1U> idx;
1522  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
1523 
1524  assign( r[I] , dot( idx , (x[I+1] - x[I-1]) ) );
1525  return r;
1526 }
1527 //----------------------------------------------------------------------
1528 template < class T, class MFLOAT >
1529 Field<T,2U,UniformCartesian<2U,MFLOAT>,Cell>&
1530 Div(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
1531  Field<T,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
1532 {
1533  const NDIndex<2U>& domain = r.getDomain();
1534  Index I = domain[0];
1535  Index J = domain[1];
1536  Vektor<T,2U> idx,idy;
1537  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
1538  idx[1] = 0.0;
1539  idy[0] = 0.0;
1540  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
1541 
1542  assign( r[I][J] ,
1543  dot( idx , (x[I+1][J ] - x[I-1][J ])) +
1544  dot( idy , (x[I ][J+1] - x[I ][J-1])) );
1545  return r;
1546 }
1547 //----------------------------------------------------------------------
1548 template < class T, class MFLOAT >
1549 Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>&
1550 Div(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
1551  Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
1552 {
1553  const NDIndex<3U>& domain = r.getDomain();
1554  Index I = domain[0];
1555  Index J = domain[1];
1556  Index K = domain[2];
1557  Vektor<T,3U> idx,idy,idz;
1558  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
1559  idx[1] = 0.0;
1560  idx[2] = 0.0;
1561  idy[0] = 0.0;
1562  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
1563  idy[2] = 0.0;
1564  idz[0] = 0.0;
1565  idz[1] = 0.0;
1566  idz[2] = 0.5/x.get_mesh().get_meshSpacing(2);
1567 
1568  assign( r[I][J][K] ,
1569  dot(idx , (x[I+1][J ][K ] - x[I-1][J ][K ] )) +
1570  dot(idy , (x[I ][J+1][K ] - x[I ][J-1][K ] )) +
1571  dot(idz , (x[I ][J ][K+1] - x[I ][J ][K-1] )) );
1572  return r;
1573 }
1574 //----------------------------------------------------------------------
1575 // Divergence Tenzor/Vert -> Vektor/Cell
1576 //----------------------------------------------------------------------
1577 template < class T, class MFLOAT >
1578 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>&
1579 Div(Field<Tenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
1580  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
1581 {
1582  const NDIndex<1U>& domain = r.getDomain();
1583  Index I = domain[0];
1584 
1585  assign( r[I] ,
1586  dot(x[I ], x.get_mesh().Dvc[0]) +
1587  dot(x[I+1], x.get_mesh().Dvc[1]));
1588  return r;
1589 }
1590 //----------------------------------------------------------------------
1591 template < class T, class MFLOAT >
1592 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>&
1593 Div(Field<Tenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
1594  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
1595 {
1596  const NDIndex<2U>& domain = r.getDomain();
1597  Index I = domain[0];
1598  Index J = domain[1];
1599 
1600  assign( r[I][J] ,
1601  dot(x[I ][J ], x.get_mesh().Dvc[0]) +
1602  dot(x[I+1][J ], x.get_mesh().Dvc[1]) +
1603  dot(x[I ][J+1], x.get_mesh().Dvc[2]) +
1604  dot(x[I+1][J+1], x.get_mesh().Dvc[3]));
1605 
1606  return r;
1607 }
1608 //----------------------------------------------------------------------
1609 template < class T, class MFLOAT >
1610 inline Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>&
1611 Div(Field<Tenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
1612  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
1613 {
1614  const NDIndex<3U>& domain = r.getDomain();
1615  Index I = domain[0];
1616  Index J = domain[1];
1617  Index K = domain[2];
1618 
1619  assign( r[I][J][K] ,
1620  dot(x[I ][J ][K ], x.get_mesh().Dvc[0]) +
1621  dot(x[I+1][J ][K ], x.get_mesh().Dvc[1]) +
1622  dot(x[I ][J+1][K ], x.get_mesh().Dvc[2]) +
1623  dot(x[I+1][J+1][K ], x.get_mesh().Dvc[3]) +
1624  dot(x[I ][J ][K+1], x.get_mesh().Dvc[4]) +
1625  dot(x[I+1][J ][K+1], x.get_mesh().Dvc[5]) +
1626  dot(x[I ][J+1][K+1], x.get_mesh().Dvc[6]) +
1627  dot(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]));
1628 
1629  return r;
1630 }
1631 //----------------------------------------------------------------------
1632 // Divergence SymTenzor/Vert -> Vektor/Cell
1633 //----------------------------------------------------------------------
1634 template < class T, class MFLOAT >
1635 inline Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>&
1636 Div(Field<SymTenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
1637  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
1638 {
1639  const NDIndex<1U>& domain = r.getDomain();
1640  Index I = domain[0];
1641 
1642  assign( r[I] ,
1643  dot(x[I ], x.get_mesh().Dvc[0]) +
1644  dot(x[I+1], x.get_mesh().Dvc[1]));
1645  return r;
1646 }
1647 //----------------------------------------------------------------------
1648 template < class T, class MFLOAT >
1649 inline Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>&
1650 Div(Field<SymTenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
1651  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
1652 {
1653  const NDIndex<2U>& domain = r.getDomain();
1654  Index I = domain[0];
1655  Index J = domain[1];
1656 
1657  assign( r[I][J] ,
1658  dot(x[I ][J ], x.get_mesh().Dvc[0]) +
1659  dot(x[I+1][J ], x.get_mesh().Dvc[1]) +
1660  dot(x[I ][J+1], x.get_mesh().Dvc[2]) +
1661  dot(x[I+1][J+1], x.get_mesh().Dvc[3]));
1662  return r;
1663 }
1664 //----------------------------------------------------------------------
1665 template < class T, class MFLOAT >
1666 inline Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>&
1667 Div(Field<SymTenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
1668  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
1669 {
1670  const NDIndex<3U>& domain = r.getDomain();
1671  Index I = domain[0];
1672  Index J = domain[1];
1673  Index K = domain[2];
1674 
1675  assign( r[I][J][K] ,
1676  dot(x[I ][J ][K ], x.get_mesh().Dvc[0]) +
1677  dot(x[I+1][J ][K ], x.get_mesh().Dvc[1]) +
1678  dot(x[I ][J+1][K ], x.get_mesh().Dvc[2]) +
1679  dot(x[I+1][J+1][K ], x.get_mesh().Dvc[3]) +
1680  dot(x[I ][J ][K+1], x.get_mesh().Dvc[4]) +
1681  dot(x[I+1][J ][K+1], x.get_mesh().Dvc[5]) +
1682  dot(x[I ][J+1][K+1], x.get_mesh().Dvc[6]) +
1683  dot(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]));
1684  return r;
1685 }
1686 
1687 //----------------------------------------------------------------------
1688 // Divergence Tenzor/Cell -> Vektor/Vert
1689 //----------------------------------------------------------------------
1690 template < class T, class MFLOAT >
1691 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>&
1692 Div(Field<Tenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
1693  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
1694 {
1695  const NDIndex<1U>& domain = r.getDomain();
1696  Index I = domain[0];
1697 
1698  assign( r[I] ,
1699  dot(x[I-1], x.get_mesh().Dvc[0]) +
1700  dot(x[I ], x.get_mesh().Dvc[1]));
1701  return r;
1702 }
1703 //----------------------------------------------------------------------
1704 template < class T, class MFLOAT >
1705 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>&
1706 Div(Field<Tenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
1707  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
1708 {
1709  const NDIndex<2U>& domain = r.getDomain();
1710  Index I = domain[0];
1711  Index J = domain[1];
1712 
1713  assign( r[I][J] ,
1714  dot(x[I-1][J-1], x.get_mesh().Dvc[0]) +
1715  dot(x[I ][J-1], x.get_mesh().Dvc[1]) +
1716  dot(x[I-1][J ], x.get_mesh().Dvc[2]) +
1717  dot(x[I ][J ], x.get_mesh().Dvc[3]));
1718  return r;
1719 }
1720 //----------------------------------------------------------------------
1721 template < class T, class MFLOAT >
1722 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>&
1723 Div(Field<Tenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
1724  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
1725 {
1726  const NDIndex<3U>& domain = r.getDomain();
1727  Index I = domain[0];
1728  Index J = domain[1];
1729  Index K = domain[2];
1730 
1731  assign( r[I][J][K] ,
1732  dot(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]) +
1733  dot(x[I ][J-1][K-1], x.get_mesh().Dvc[1]) +
1734  dot(x[I-1][J ][K-1], x.get_mesh().Dvc[2]) +
1735  dot(x[I ][J ][K-1], x.get_mesh().Dvc[3]) +
1736  dot(x[I-1][J-1][K ], x.get_mesh().Dvc[4]) +
1737  dot(x[I ][J-1][K ], x.get_mesh().Dvc[5]) +
1738  dot(x[I-1][J ][K ], x.get_mesh().Dvc[6]) +
1739  dot(x[I ][J ][K ], x.get_mesh().Dvc[7]));
1740  return r;
1741 }
1742 
1743 //----------------------------------------------------------------------
1744 // Divergence SymTenzor/Cell -> Vektor/Vert
1745 //----------------------------------------------------------------------
1746 template < class T, class MFLOAT >
1747 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>&
1748 Div(Field<SymTenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
1749  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
1750 {
1751  const NDIndex<1U>& domain = r.getDomain();
1752  Index I = domain[0];
1753 
1754  assign( r[I] ,
1755  dot(x[I-1], x.get_mesh().Dvc[0]) +
1756  dot(x[I ], x.get_mesh().Dvc[1]));
1757  return r;
1758 }
1759 //----------------------------------------------------------------------
1760 template < class T, class MFLOAT >
1761 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>&
1762 Div(Field<SymTenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
1763  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
1764 {
1765  const NDIndex<2U>& domain = r.getDomain();
1766  Index I = domain[0];
1767  Index J = domain[1];
1768 
1769  assign( r[I][J] ,
1770  dot(x[I-1][J-1], x.get_mesh().Dvc[0]) +
1771  dot(x[I ][J-1], x.get_mesh().Dvc[1]) +
1772  dot(x[I-1][J ], x.get_mesh().Dvc[2]) +
1773  dot(x[I ][J ], x.get_mesh().Dvc[3]));
1774  return r;
1775 }
1776 //----------------------------------------------------------------------
1777 template < class T, class MFLOAT >
1778 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>&
1779 Div(Field<SymTenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
1780  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
1781 {
1782  const NDIndex<3U>& domain = r.getDomain();
1783  Index I = domain[0];
1784  Index J = domain[1];
1785  Index K = domain[2];
1786 
1787  assign( r[I][J][K] ,
1788  dot(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]) +
1789  dot(x[I ][J-1][K-1], x.get_mesh().Dvc[1]) +
1790  dot(x[I-1][J ][K-1], x.get_mesh().Dvc[2]) +
1791  dot(x[I ][J ][K-1], x.get_mesh().Dvc[3]) +
1792  dot(x[I-1][J-1][K ], x.get_mesh().Dvc[4]) +
1793  dot(x[I ][J-1][K ], x.get_mesh().Dvc[5]) +
1794  dot(x[I-1][J ][K ], x.get_mesh().Dvc[6]) +
1795  dot(x[I ][J ][K ], x.get_mesh().Dvc[7]));
1796  return r;
1797 }
1798 
1799 //----------------------------------------------------------------------
1800 // Grad Scalar/Vert -> Vektor/Cell
1801 //----------------------------------------------------------------------
1802 template < class T, class MFLOAT >
1803 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>&
1804 Grad(Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
1805  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
1806 {
1807  const NDIndex<1U>& domain = r.getDomain();
1808  Index I = domain[0];
1809 
1810  assign( r[I] ,
1811  x[I ] * x.get_mesh().Dvc[0] +
1812  x[I+1] * x.get_mesh().Dvc[1]);
1813  return r;
1814 }
1815 //----------------------------------------------------------------------
1816 template < class T, class MFLOAT >
1817 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>&
1818 Grad(Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
1819  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
1820 {
1821  const NDIndex<2U>& domain = r.getDomain();
1822  Index I = domain[0];
1823  Index J = domain[1];
1824 
1825  assign( r[I][J] ,
1826  x[I ][J ] * x.get_mesh().Dvc[0] +
1827  x[I+1][J ] * x.get_mesh().Dvc[1] +
1828  x[I ][J+1] * x.get_mesh().Dvc[2] +
1829  x[I+1][J+1] * x.get_mesh().Dvc[3]);
1830  return r;
1831 }
1832 //----------------------------------------------------------------------
1833 template < class T, class MFLOAT >
1834 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>&
1835 Grad(Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
1836  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
1837 {
1838  const NDIndex<3U>& domain = r.getDomain();
1839  Index I = domain[0];
1840  Index J = domain[1];
1841  Index K = domain[2];
1842 
1843  assign( r[I][J][K] ,
1844  x[I ][J ][K ] * x.get_mesh().Dvc[0] +
1845  x[I+1][J ][K ] * x.get_mesh().Dvc[1] +
1846  x[I ][J+1][K ] * x.get_mesh().Dvc[2] +
1847  x[I+1][J+1][K ] * x.get_mesh().Dvc[3] +
1848  x[I ][J ][K+1] * x.get_mesh().Dvc[4] +
1849  x[I+1][J ][K+1] * x.get_mesh().Dvc[5] +
1850  x[I ][J+1][K+1] * x.get_mesh().Dvc[6] +
1851  x[I+1][J+1][K+1] * x.get_mesh().Dvc[7]);
1852 
1853  return r;
1854 }
1855 //----------------------------------------------------------------------
1856 // Grad Scalar/Vert -> Vektor/Edge
1857 //----------------------------------------------------------------------
1858 template < class T, class MFLOAT >
1859 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Edge>&
1860 Grad(Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
1861  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Edge>& r)
1862 {
1863  const NDIndex<1U>& domain = r.getDomain();
1864  Index I = Index(domain[0].first(), domain[0].last()-1);
1865 
1866  assign( r[I] ,
1867  x[I ] * x.get_mesh().Dvc[0] +
1868  x[I+1] * x.get_mesh().Dvc[1]);
1869  return r;
1870 }
1871 //----------------------------------------------------------------------
1872 template < class T, class MFLOAT >
1873 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Edge>&
1874 Grad(Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
1875  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Edge>& r)
1876 {
1877  const NDIndex<2U>& domain = r.getDomain();
1878  Index I = Index(domain[0].first(), domain[0].last()-1);
1879  Index J = Index(domain[1].first(), domain[1].last()-1);
1880 
1881  Vektor<T,2U> idx,idy;
1882  idx[0] = 1.0/x.get_mesh().get_meshSpacing(0);
1883  idx[1] = 0.0;
1884  idy[0] = 0.0;
1885  idy[1] = 1.0/x.get_mesh().get_meshSpacing(1);
1886 
1887  assign( r[I][J] ,
1888  (x[I+1][J ] - x[I ][J ]) * idx +
1889  (x[I ][J+1] - x[I ][J ]) * idy);
1890  I = Index(domain[0].last(), domain[0].last());
1891  assign( r[I][J](1),
1892  (x[I][J+1] - x[I][J]));
1893  I = Index(domain[0].first(), domain[0].last()-1);
1894  J = Index(domain[1].last(), domain[1].last());
1895  assign( r[I][J](0),
1896  (x[I+1][J] - x[I][J]));
1897  J = Index(domain[1].first(), domain[1].last()-1);
1898  return r;
1899 }
1900 //----------------------------------------------------------------------
1901 template < class T, class MFLOAT >
1902 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Edge>&
1903 Grad(Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
1904  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Edge>& r)
1905 {
1906  const NDIndex<3U>& domain = r.getDomain();
1907  Index I = Index(domain[0].first(), domain[0].last()-1);
1908  Index J = Index(domain[1].first(), domain[1].last()-1);
1909  Index K = Index(domain[2].first(), domain[2].last()-1);
1910 
1911  Vektor<T,3U> idx,idy,idz;
1912  idx[0] = 1.0/x.get_mesh().get_meshSpacing(0);
1913  idx[1] = 0.0;
1914  idx[2] = 0.0;
1915  idy[0] = 0.0;
1916  idy[1] = 1.0/x.get_mesh().get_meshSpacing(1);
1917  idy[2] = 0.0;
1918  idz[0] = 0.0;
1919  idz[1] = 0.0;
1920  idz[2] = 1.0/x.get_mesh().get_meshSpacing(2);
1921 
1922  assign( r[I][J][K] ,
1923  (x[I+1][J ][K ] - x[I ][J ][K ]) * idx +
1924  (x[I ][J+1][K ] - x[I ][J ][K ]) * idy +
1925  (x[I ][J ][K+1] - x[I ][J ][K ]) * idz);
1926  I = Index(domain[0].last(), domain[0].last());
1927  assign( r[I][J][K](1),
1928  (x[I ][J+1][K ] - x[I ][J ][K ]));
1929  assign( r[I][J][K](2),
1930  (x[I ][J ][K+1] - x[I ][J ][K ]));
1931  I = Index(domain[0].first(), domain[0].last()-1);
1932  J = Index(domain[1].last(), domain[1].last());
1933  assign( r[I][J][K](0),
1934  (x[I+1][J ][K ] - x[I ][J ][K ]));
1935  assign( r[I][J][K](2),
1936  (x[I ][J ][K+1] - x[I ][J ][K ]));
1937  J = Index(domain[1].first(), domain[1].last()-1);
1938  K = Index(domain[2].last(), domain[2].last());
1939  assign( r[I][J][K](0),
1940  (x[I+1][J ][K ] - x[I ][J ][K ]));
1941  assign( r[I][J][K](1),
1942  (x[I ][J+1][K ] - x[I ][J ][K ]));
1943  K = Index(domain[2].first(), domain[2].last()-1);
1944 
1945  return r;
1946 }
1947 //----------------------------------------------------------------------
1948 // Grad Scalar/Cell -> Vektor/Vert
1949 //----------------------------------------------------------------------
1950 template < class T, class MFLOAT >
1951 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>&
1952 Grad(Field<T,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
1953  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
1954 {
1955  const NDIndex<1U>& domain = r.getDomain();
1956  Index I = domain[0];
1957 
1958  assign( r[I] ,
1959  x[I-1] * x.get_mesh().Dvc[0] +
1960  x[I ] * x.get_mesh().Dvc[1]);
1961  return r;
1962 }
1963 //----------------------------------------------------------------------
1964 template < class T, class MFLOAT >
1965 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>&
1966 Grad(Field<T,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
1967  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
1968 {
1969  const NDIndex<2U>& domain = r.getDomain();
1970  Index I = domain[0];
1971  Index J = domain[1];
1972 
1973  assign( r[I][J] ,
1974  x[I-1][J-1] * x.get_mesh().Dvc[0] +
1975  x[I ][J-1] * x.get_mesh().Dvc[1] +
1976  x[I-1][J ] * x.get_mesh().Dvc[2] +
1977  x[I ][J ] * x.get_mesh().Dvc[3]);
1978  return r;
1979 }
1980 //----------------------------------------------------------------------
1981 template < class T, class MFLOAT >
1982 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>&
1983 Grad(Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
1984  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
1985 {
1986  const NDIndex<3U>& domain = r.getDomain();
1987  Index I = domain[0];
1988  Index J = domain[1];
1989  Index K = domain[2];
1990 
1991  assign( r[I][J][K] ,
1992  x[I-1][J-1][K-1] * x.get_mesh().Dvc[0] +
1993  x[I ][J-1][K-1] * x.get_mesh().Dvc[1] +
1994  x[I-1][J ][K-1] * x.get_mesh().Dvc[2] +
1995  x[I ][J ][K-1] * x.get_mesh().Dvc[3] +
1996  x[I-1][J-1][K ] * x.get_mesh().Dvc[4] +
1997  x[I ][J-1][K ] * x.get_mesh().Dvc[5] +
1998  x[I-1][J ][K ] * x.get_mesh().Dvc[6] +
1999  x[I ][J ][K ] * x.get_mesh().Dvc[7]);
2000  return r;
2001 }
2002 //----------------------------------------------------------------------
2003 // Grad Scalar/Vert -> Vektor/Vert
2004 //----------------------------------------------------------------------
2005 template < class T, class MFLOAT >
2006 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>&
2007 Grad(Field<T,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
2008  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
2009 {
2010  const NDIndex<1U>& domain = r.getDomain();
2011  Index I = domain[0];
2012 
2013  Vektor<T,1U> idx;
2014  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2015 
2016  assign( r[I] , idx * (x[I+1] - x[I-1] ) );
2017  return r;
2018 }
2019 //----------------------------------------------------------------------
2020 template < class T, class MFLOAT >
2021 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>&
2022 Grad(Field<T,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
2023  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
2024 {
2025  const NDIndex<2U>& domain = r.getDomain();
2026  Index I = domain[0];
2027  Index J = domain[1];
2028 
2029  Vektor<T,2U> idx,idy;
2030  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2031  idx[1] = 0.0;
2032  idy[0] = 0.0;
2033  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
2034 
2035  assign( r[I][J] ,
2036  idx * (x[I+1][J ] - x[I-1][J ]) +
2037  idy * (x[I ][J+1] - x[I ][J-1]) );
2038  return r;
2039 }
2040 //----------------------------------------------------------------------
2041 template < class T, class MFLOAT >
2042 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>&
2043 Grad(Field<T,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
2044  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
2045 {
2046  const NDIndex<3U>& domain = r.getDomain();
2047  Index I = domain[0];
2048  Index J = domain[1];
2049  Index K = domain[2];
2050 
2051  Vektor<T,3U> idx,idy,idz;
2052  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2053  idx[1] = 0.0;
2054  idx[2] = 0.0;
2055  idy[0] = 0.0;
2056  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
2057  idy[2] = 0.0;
2058  idz[0] = 0.0;
2059  idz[1] = 0.0;
2060  idz[2] = 0.5/x.get_mesh().get_meshSpacing(2);
2061 
2062  assign(r[I][J][K] ,
2063  idx * (x[I+1][J ][K ] - x[I-1][J ][K ]) +
2064  idy * (x[I ][J+1][K ] - x[I ][J-1][K ]) +
2065  idz * (x[I ][J ][K+1] - x[I ][J ][K-1]));
2066  return r;
2067 }
2068 //----------------------------------------------------------------------
2069 // Grad Scalar/Cell -> Vektor/Cell
2070 //----------------------------------------------------------------------
2071 template < class T, class MFLOAT >
2072 Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>&
2073 Grad(Field<T,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
2074  Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
2075 {
2076  const NDIndex<1U>& domain = r.getDomain();
2077  Index I = domain[0];
2078 
2079  Vektor<T,1U> idx;
2080  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2081 
2082  assign( r[I] , idx * (x[I+1] - x[I-1] ) );
2083  return r;
2084 }
2085 //----------------------------------------------------------------------
2086 template < class T, class MFLOAT >
2087 Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>&
2088 Grad(Field<T,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
2089  Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
2090 {
2091  const NDIndex<2U>& domain = r.getDomain();
2092  Index I = domain[0];
2093  Index J = domain[1];
2094 
2095  Vektor<T,2U> idx,idy;
2096  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2097  idx[1] = 0.0;
2098  idy[0] = 0.0;
2099  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
2100 
2101  assign( r[I][J] ,
2102  idx * (x[I+1][J ] - x[I-1][J ]) +
2103  idy * (x[I ][J+1] - x[I ][J-1]) );
2104 
2105  return r;
2106 }
2107 //----------------------------------------------------------------------
2108 template < class T, class MFLOAT >
2109 Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>&
2110 Grad(Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
2111  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
2112 {
2113  const NDIndex<3U>& domain = r.getDomain();
2114  Index I = domain[0];
2115  Index J = domain[1];
2116  Index K = domain[2];
2117 
2118  Vektor<T,3U> idx,idy,idz;
2119  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2120  idx[1] = 0.0;
2121  idx[2] = 0.0;
2122  idy[0] = 0.0;
2123  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
2124  idy[2] = 0.0;
2125  idz[0] = 0.0;
2126  idz[1] = 0.0;
2127  idz[2] = 0.5/x.get_mesh().get_meshSpacing(2);
2128 
2129  assign(r[I][J][K] ,
2130  idx * (x[I+1][J ][K ] - x[I-1][J ][K ]) +
2131  idy * (x[I ][J+1][K ] - x[I ][J-1][K ]) +
2132  idz * (x[I ][J ][K+1] - x[I ][J ][K-1]));
2133 
2134  //approximate at the border up to a term of order O(h^2)
2135  Vektor<T,3U> xx = Vektor<T,3U>(1.0,0.0,0.0);
2136  Vektor<T,3U> yy = Vektor<T,3U>(0.0,1.0,0.0);
2137  Vektor<T,3U> zz = Vektor<T,3U>(0.0,0.0,1.0);
2138 
2139  Index lo(0, 0);
2140 
2141  int maxI = I.length() - 1;
2142  int maxJ = J.length() - 1;
2143  int maxK = K.length() - 1;
2144 
2145  Index Iup(maxI, maxI);
2146  Index Jup(maxJ, maxJ);
2147  Index Kup(maxK, maxK);
2148 
2149  r[lo][J][K] = (idx * (- 1.0*x[lo+2][J][K]
2150  + 4.0*x[lo+1][J][K]
2151  - 3.0*x[lo ][J][K])
2152  + yy*r[lo][J][K]
2153  + zz*r[lo][J][K]);
2154  r[Iup][J][K] = (idx * ( 1.0*x[Iup-2][J][K]
2155  - 4.0*x[Iup-1][J][K]
2156  + 3.0*x[Iup ][J][K])
2157  + yy*r[Iup][J][K]
2158  + zz*r[Iup][J][K]);
2159 
2160  r[I][lo][K] = (idy * (- 1.0*x[I][lo+2][K]
2161  + 4.0*x[I][lo+1][K]
2162  - 3.0*x[I][lo ][K])
2163  + xx*r[I][lo][K]
2164  + zz*r[I][lo][K]);
2165  r[I][Jup][K] = (idy * ( 1.0*x[I][Jup-2][K]
2166  - 4.0*x[I][Jup-1][K]
2167  + 3.0*x[I][Jup ][K])
2168  + xx*r[I][Jup][K]
2169  + zz*r[I][Jup][K]);
2170 
2171  r[I][J][lo] = (idz * (- 1.0*x[I][J][lo+2]
2172  + 4.0*x[I][J][lo+1]
2173  - 3.0*x[I][J][lo ])
2174  + xx*r[I][J][lo]
2175  + yy*r[I][J][lo]);
2176  r[I][J][Kup] = (idz * ( 1.0*x[I][J][Kup-2]
2177  - 4.0*x[I][J][Kup-1]
2178  + 3.0*x[I][J][Kup ])
2179  + xx*r[I][J][Kup]
2180  + yy*r[I][J][Kup]);
2181 
2182  return r;
2183 }
2184 
2185 
2186 template < class T, class MFLOAT >
2187 inline Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>&
2188 Grad1Ord(Field<T,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
2189  Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
2190 {
2191  const NDIndex<3U>& domain = r.getDomain();
2192  Index I = domain[0];
2193  Index J = domain[1];
2194  Index K = domain[2];
2195 
2196  Vektor<T,3U> idx,idy,idz;
2197  idx[0] = 0.5/x.get_mesh().get_meshSpacing(0);
2198  idx[1] = 0.0;
2199  idx[2] = 0.0;
2200  idy[0] = 0.0;
2201  idy[1] = 0.5/x.get_mesh().get_meshSpacing(1);
2202  idy[2] = 0.0;
2203  idz[0] = 0.0;
2204  idz[1] = 0.0;
2205  idz[2] = 0.5/x.get_mesh().get_meshSpacing(2);
2206 
2207  assign(r[I][J][K] ,
2208  idx * (x[I+1][J ][K ] - x[I-1][J ][K ]) +
2209  idy * (x[I ][J+1][K ] - x[I ][J-1][K ]) +
2210  idz * (x[I ][J ][K+1] - x[I ][J ][K-1]));
2211  return r;
2212 }
2213 //----------------------------------------------------------------------
2214 // Grad Vektor/Vert -> Tenzor/Cell
2215 //----------------------------------------------------------------------
2216 template < class T, class MFLOAT >
2217 Field<Tenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>&
2218 Grad(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
2219  Field<Tenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
2220 {
2221  const NDIndex<1U>& domain = r.getDomain();
2222  Index I = domain[0];
2223 
2224  assign( r[I] ,
2225  outerProduct( x[I ] , x.get_mesh().Dvc[0] ) +
2226  outerProduct( x[I+1] , x.get_mesh().Dvc[1])) ;
2227  return r;
2228 }
2229 //----------------------------------------------------------------------
2230 template < class T, class MFLOAT >
2231 Field<Tenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>&
2232 Grad(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
2233  Field<Tenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
2234 {
2235  const NDIndex<2U>& domain = r.getDomain();
2236  Index I = domain[0];
2237  Index J = domain[1];
2238 
2239  assign( r[I][J] ,
2240  outerProduct( x[I ][J ] , x.get_mesh().Dvc[0] ) +
2241  outerProduct( x[I+1][J ] , x.get_mesh().Dvc[1] ) +
2242  outerProduct( x[I ][J+1] , x.get_mesh().Dvc[2] ) +
2243  outerProduct( x[I+1][J+1] , x.get_mesh().Dvc[3] )) ;
2244  return r;
2245 }
2246 //----------------------------------------------------------------------
2247 template < class T, class MFLOAT >
2248 Field<Tenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>&
2249 Grad(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
2250  Field<Tenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
2251 {
2252  const NDIndex<3U>& domain = r.getDomain();
2253  Index I = domain[0];
2254  Index J = domain[1];
2255  Index K = domain[2];
2256 
2257  assign( r[I][J][K] ,
2258  outerProduct( x[I ][J ][K ] , x.get_mesh().Dvc[0] ) +
2259  outerProduct( x[I+1][J ][K ] , x.get_mesh().Dvc[1] ) +
2260  outerProduct( x[I ][J+1][K ] , x.get_mesh().Dvc[2] ) +
2261  outerProduct( x[I+1][J+1][K ] , x.get_mesh().Dvc[3] ) +
2262  outerProduct( x[I ][J ][K+1] , x.get_mesh().Dvc[4] ) +
2263  outerProduct( x[I+1][J ][K+1] , x.get_mesh().Dvc[5] ) +
2264  outerProduct( x[I ][J+1][K+1] , x.get_mesh().Dvc[6] ) +
2265  outerProduct( x[I+1][J+1][K+1] , x.get_mesh().Dvc[7] ));
2266 
2267  return r;
2268 }
2269 //----------------------------------------------------------------------
2270 // Grad Vektor/Cell -> Tenzor/Vert
2271 //----------------------------------------------------------------------
2272 template < class T, class MFLOAT >
2273 Field<Tenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>&
2274 Grad(Field<Vektor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
2275  Field<Tenzor<T,1U>,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
2276 {
2277  const NDIndex<1U>& domain = r.getDomain();
2278  Index I = domain[0];
2279 
2280  assign( r[I] ,
2281  outerProduct( x[I-1] , x.get_mesh().Dvc[0] ) +
2282  outerProduct( x[I ] , x.get_mesh().Dvc[1] ));
2283  return r;
2284 }
2285 //----------------------------------------------------------------------
2286 template < class T, class MFLOAT >
2287 Field<Tenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>&
2288 Grad(Field<Vektor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
2289  Field<Tenzor<T,2U>,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
2290 {
2291  const NDIndex<2U>& domain = r.getDomain();
2292  Index I = domain[0];
2293  Index J = domain[1];
2294 
2295  assign( r[I][J] ,
2296  outerProduct( x[I-1][J-1] , x.get_mesh().Dvc[0] ) +
2297  outerProduct( x[I ][J-1] , x.get_mesh().Dvc[1] ) +
2298  outerProduct( x[I-1][J ] , x.get_mesh().Dvc[2] ) +
2299  outerProduct( x[I ][J ] , x.get_mesh().Dvc[3] ));
2300  return r;
2301 }
2302 //----------------------------------------------------------------------
2303 template < class T, class MFLOAT >
2304 Field<Tenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>&
2305 Grad(Field<Vektor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
2306  Field<Tenzor<T,3U>,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
2307 {
2308  const NDIndex<3U>& domain = r.getDomain();
2309  Index I = domain[0];
2310  Index J = domain[1];
2311  Index K = domain[2];
2312 
2313  assign( r[I][J][K] ,
2314  outerProduct( x[I-1][J-1][K-1] , x.get_mesh().Dvc[0] ) +
2315  outerProduct( x[I ][J-1][K-1] , x.get_mesh().Dvc[1] ) +
2316  outerProduct( x[I-1][J ][K-1] , x.get_mesh().Dvc[2] ) +
2317  outerProduct( x[I ][J ][K-1] , x.get_mesh().Dvc[3] ) +
2318  outerProduct( x[I-1][J-1][K ] , x.get_mesh().Dvc[4] ) +
2319  outerProduct( x[I ][J-1][K ] , x.get_mesh().Dvc[5] ) +
2320  outerProduct( x[I-1][J ][K ] , x.get_mesh().Dvc[6] ) +
2321  outerProduct( x[I ][J ][K ] , x.get_mesh().Dvc[7] ));
2322  return r;
2323 }
2324 
2325 namespace IPPL {
2326 
2327 //----------------------------------------------------------------------
2328 // Weighted average Cell to Vert
2329 //----------------------------------------------------------------------
2330 template < class T1, class T2, class MFLOAT >
2331 Field<T1,1U,UniformCartesian<1U,MFLOAT>,Vert>&
2332 Average(Field<T1,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
2333  Field<T2,1U,UniformCartesian<1U,MFLOAT>,Cell>& w,
2334  Field<T1,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
2335 {
2336  const NDIndex<1U>& domain = r.getDomain();
2337  Index I = domain[0];
2338  assign( r[I] ,
2339  ( x[I-1] * w[I-1] + x[I ] * w[I ] )/
2340  ( w[I-1] + w[I ] ) );
2341  return r;
2342 }
2343 //----------------------------------------------------------------------
2344 template < class T1, class T2, class MFLOAT >
2345 Field<T1,2U,UniformCartesian<2U,MFLOAT>,Vert>&
2346 Average(Field<T1,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
2347  Field<T2,2U,UniformCartesian<2U,MFLOAT>,Cell>& w,
2348  Field<T1,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
2349 {
2350  const NDIndex<2U>& domain = r.getDomain();
2351  Index I = domain[0];
2352  Index J = domain[1];
2353 
2354  assign( r[I][J] ,
2355  ( x[I-1][J-1] * w[I-1][J-1] +
2356  x[I ][J-1] * w[I ][J-1] +
2357  x[I-1][J ] * w[I-1][J ] +
2358  x[I ][J ] * w[I ][J ] )/
2359  ( w[I-1][J-1] +
2360  w[I ][J-1] +
2361  w[I-1][J ] +
2362  w[I ][J ] ) );
2363  return r;
2364 }
2365 //----------------------------------------------------------------------
2366 template < class T1, class T2, class MFLOAT >
2367 Field<T1,3U,UniformCartesian<3U,MFLOAT>,Vert>&
2368 Average(Field<T1,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
2369  Field<T2,3U,UniformCartesian<3U,MFLOAT>,Cell>& w,
2370  Field<T1,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
2371 {
2372  const NDIndex<3U>& domain = r.getDomain();
2373  Index I = domain[0];
2374  Index J = domain[1];
2375  Index K = domain[2];
2376 
2377  assign( r[I][J][K] ,
2378  ( x[I-1][J-1][K-1] * w[I-1][J-1][K-1] +
2379  x[I ][J-1][K-1] * w[I ][J-1][K-1] +
2380  x[I-1][J ][K-1] * w[I-1][J ][K-1] +
2381  x[I ][J ][K-1] * w[I ][J ][K-1] +
2382  x[I-1][J-1][K ] * w[I-1][J-1][K ] +
2383  x[I ][J-1][K ] * w[I ][J-1][K ] +
2384  x[I-1][J ][K ] * w[I-1][J ][K ] +
2385  x[I ][J ][K ] * w[I ][J ][K ] )/
2386  ( w[I-1][J-1][K-1] +
2387  w[I ][J-1][K-1] +
2388  w[I-1][J ][K-1] +
2389  w[I ][J ][K-1] +
2390  w[I-1][J-1][K ] +
2391  w[I ][J-1][K ] +
2392  w[I-1][J ][K ] +
2393  w[I ][J ][K ] ) );
2394  return r;
2395 }
2396 //----------------------------------------------------------------------
2397 // Weighted average Vert to Cell
2398 // N.B.: won't work except for unit-stride (& zero-base?) Field's.
2399 //----------------------------------------------------------------------
2400 template < class T1, class T2, class MFLOAT >
2401 Field<T1,1U,UniformCartesian<1U,MFLOAT>,Cell>&
2402 Average(Field<T1,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
2403  Field<T2,1U,UniformCartesian<1U,MFLOAT>,Vert>& w,
2404  Field<T1,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
2405 {
2406  const NDIndex<1U>& domain = r.getDomain();
2407  Index I = domain[0];
2408  assign( r[I] ,
2409  ( x[I+1] * w[I+1] + x[I ] * w[I ] )/
2410  ( w[I+1] + w[I ] ) );
2411  return r;
2412 }
2413 //----------------------------------------------------------------------
2414 template < class T1, class T2, class MFLOAT >
2415 Field<T1,2U,UniformCartesian<2U,MFLOAT>,Cell>&
2416 Average(Field<T1,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
2417  Field<T2,2U,UniformCartesian<2U,MFLOAT>,Vert>& w,
2418  Field<T1,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
2419 {
2420  const NDIndex<2U>& domain = r.getDomain();
2421  Index I = domain[0];
2422  Index J = domain[1];
2423 
2424  assign( r[I][J] ,
2425  ( x[I+1][J+1] * w[I+1][J+1] +
2426  x[I ][J+1] * w[I ][J+1] +
2427  x[I+1][J ] * w[I+1][J ] +
2428  x[I ][J ] * w[I ][J ] )/
2429  ( w[I+1][J+1] +
2430  w[I ][J+1] +
2431  w[I+1][J ] +
2432  w[I ][J ] ) );
2433  return r;
2434 }
2435 //----------------------------------------------------------------------
2436 template < class T1, class T2, class MFLOAT >
2437 Field<T1,3U,UniformCartesian<3U,MFLOAT>,Cell>&
2438 Average(Field<T1,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
2439  Field<T2,3U,UniformCartesian<3U,MFLOAT>,Vert>& w,
2440  Field<T1,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
2441 {
2442  const NDIndex<3U>& domain = r.getDomain();
2443  Index I = domain[0];
2444  Index J = domain[1];
2445  Index K = domain[2];
2446 
2447  assign( r[I][J][K] ,
2448  ( x[I+1][J+1][K+1] * w[I+1][J+1][K+1] +
2449  x[I ][J+1][K+1] * w[I ][J+1][K+1] +
2450  x[I+1][J ][K+1] * w[I+1][J ][K+1] +
2451  x[I ][J ][K+1] * w[I ][J ][K+1] +
2452  x[I+1][J+1][K ] * w[I+1][J+1][K ] +
2453  x[I ][J+1][K ] * w[I ][J+1][K ] +
2454  x[I+1][J ][K ] * w[I+1][J ][K ] +
2455  x[I ][J ][K ] * w[I ][J ][K ] )/
2456  ( w[I+1][J+1][K+1] +
2457  w[I ][J+1][K+1] +
2458  w[I+1][J ][K+1] +
2459  w[I ][J ][K+1] +
2460  w[I+1][J+1][K ] +
2461  w[I ][J+1][K ] +
2462  w[I+1][J ][K ] +
2463  w[I ][J ][K ] ) );
2464  return r;
2465 }
2466 
2467 //----------------------------------------------------------------------
2468 // Unweighted average Cell to Vert
2469 //----------------------------------------------------------------------
2470 template < class T1, class MFLOAT >
2471 Field<T1,1U,UniformCartesian<1U,MFLOAT>,Vert>&
2472 Average(Field<T1,1U,UniformCartesian<1U,MFLOAT>,Cell>& x,
2473  Field<T1,1U,UniformCartesian<1U,MFLOAT>,Vert>& r)
2474 {
2475  const NDIndex<1U>& domain = r.getDomain();
2476  Index I = domain[0];
2477  r[I] = 0.5*(x[I-1] + x[I ]);
2478  return r;
2479 }
2480 //----------------------------------------------------------------------
2481 template < class T1, class MFLOAT >
2482 Field<T1,2U,UniformCartesian<2U,MFLOAT>,Vert>&
2483 Average(Field<T1,2U,UniformCartesian<2U,MFLOAT>,Cell>& x,
2484  Field<T1,2U,UniformCartesian<2U,MFLOAT>,Vert>& r)
2485 {
2486  const NDIndex<2U>& domain = r.getDomain();
2487  Index I = domain[0];
2488  Index J = domain[1];
2489  r[I][J] = 0.25*(x[I-1][J-1] + x[I ][J-1] + x[I-1][J ] + x[I ][J ]);
2490  return r;
2491 }
2492 //----------------------------------------------------------------------
2493 template < class T1, class MFLOAT >
2494 Field<T1,3U,UniformCartesian<3U,MFLOAT>,Vert>&
2495 Average(Field<T1,3U,UniformCartesian<3U,MFLOAT>,Cell>& x,
2496  Field<T1,3U,UniformCartesian<3U,MFLOAT>,Vert>& r)
2497 {
2498  const NDIndex<3U>& domain = r.getDomain();
2499  Index I = domain[0];
2500  Index J = domain[1];
2501  Index K = domain[2];
2502  r[I][J][K] = 0.125*(x[I-1][J-1][K-1] + x[I ][J-1][K-1] + x[I-1][J ][K-1] +
2503  x[I ][J ][K-1] + x[I-1][J-1][K ] + x[I ][J-1][K ] +
2504  x[I-1][J ][K ] + x[I ][J ][K ]);
2505  return r;
2506 }
2507 //----------------------------------------------------------------------
2508 // Unweighted average Vert to Cell
2509 // N.B.: won't work except for unit-stride (& zero-base?) Field's.
2510 //----------------------------------------------------------------------
2511 template < class T1, class MFLOAT >
2512 Field<T1,1U,UniformCartesian<1U,MFLOAT>,Cell>&
2513 Average(Field<T1,1U,UniformCartesian<1U,MFLOAT>,Vert>& x,
2514 
2515  Field<T1,1U,UniformCartesian<1U,MFLOAT>,Cell>& r)
2516 {
2517  const NDIndex<1U>& domain = r.getDomain();
2518  Index I = domain[0];
2519  r[I] = 0.5*(x[I+1] + x[I ]);
2520  return r;
2521 }
2522 //----------------------------------------------------------------------
2523 template < class T1, class MFLOAT >
2524 Field<T1,2U,UniformCartesian<2U,MFLOAT>,Cell>&
2525 Average(Field<T1,2U,UniformCartesian<2U,MFLOAT>,Vert>& x,
2526  Field<T1,2U,UniformCartesian<2U,MFLOAT>,Cell>& r)
2527 {
2528  const NDIndex<2U>& domain = r.getDomain();
2529  Index I = domain[0];
2530  Index J = domain[1];
2531  r[I][J] = 0.25*(x[I+1][J+1] + x[I ][J+1] + x[I+1][J ] + x[I ][J ]);
2532  return r;
2533 }
2534 //----------------------------------------------------------------------
2535 template < class T1, class MFLOAT >
2536 Field<T1,3U,UniformCartesian<3U,MFLOAT>,Cell>&
2537 Average(Field<T1,3U,UniformCartesian<3U,MFLOAT>,Vert>& x,
2538  Field<T1,3U,UniformCartesian<3U,MFLOAT>,Cell>& r)
2539 {
2540  const NDIndex<3U>& domain = r.getDomain();
2541  Index I = domain[0];
2542  Index J = domain[1];
2543  Index K = domain[2];
2544  r[I][J][K] = 0.125*(x[I+1][J+1][K+1] + x[I ][J+1][K+1] + x[I+1][J ][K+1] +
2545  x[I ][J ][K+1] + x[I+1][J+1][K ] + x[I ][J+1][K ] +
2546  x[I+1][J ][K ] + x[I ][J ][K ]);
2547  return r;
2548 }
2549 
2550 }
dot
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
Definition: Vector3D.cpp:118
Dim
const unsigned Dim
Definition: P3MPoissonSolver.h:25
outerProduct
Tenzor< typename PETEBinaryReturn< T1, T2, OpMultipply >::type, D > outerProduct(const Vektor< T1, D > &v1, const Vektor< T2, D > &v2)
Definition: Tenzor.h:443
assign
void assign(const BareField< T, Dim > &a, RHS b, OP op, ExprTag< true >)
AddGuardCells
NDIndex< Dim > AddGuardCells(const NDIndex< Dim > &idx, const GuardCellSizes< Dim > &g)
Definition: GuardCellSizes.h:69
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
min
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
e_dim_tag
e_dim_tag
Definition: FieldLayout.h:55
PARALLEL
@ PARALLEL
Definition: FieldLayout.h:55
Div
Field< T, 1U, UniformCartesian< 1U, MFLOAT >, Cell > & Div(Field< Vektor< T, 1U >, 1U, UniformCartesian< 1U, MFLOAT >, Vert > &x, Field< T, 1U, UniformCartesian< 1U, MFLOAT >, Cell > &r)
Definition: UniformCartesian.hpp:1279
Grad
Field< Vektor< T, 1U >, 1U, UniformCartesian< 1U, MFLOAT >, Cell > & Grad(Field< T, 1U, UniformCartesian< 1U, MFLOAT >, Vert > &x, Field< Vektor< T, 1U >, 1U, UniformCartesian< 1U, MFLOAT >, Cell > &r)
Definition: UniformCartesian.hpp:1804
Grad1Ord
Field< Vektor< T, 3U >, 3U, UniformCartesian< 3U, MFLOAT >, Cell > & Grad1Ord(Field< T, 3U, UniformCartesian< 3U, MFLOAT >, Cell > &x, Field< Vektor< T, 3U >, 3U, UniformCartesian< 3U, MFLOAT >, Cell > &r)
Definition: UniformCartesian.hpp:2188
PETE_apply
void PETE_apply(OpUMeshExtrapolate< T > e, T &a, T b)
Definition: UniformCartesian.hpp:490
a
std::complex< double > a
Definition: IpplExpressions.h:56
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
ERRORMSG
#define ERRORMSG(msg)
Definition: IpplInfo.h:350
PAssert_LT
#define PAssert_LT(a, b)
Definition: PAssert.h:106
PInsist
#define PInsist(c, m)
Definition: PAssert.h:120
Physics::e
constexpr double e
The value of.
Definition: Physics.h:39
mslang::iterator
std::string::iterator iterator
Definition: MSLang.h:16
Options::info
bool info
Info flag.
Definition: Options.cpp:28
IPPL::Average
Field< T1, 1U, Cartesian< 1U, MFLOAT >, Vert > & Average(Field< T1, 1U, Cartesian< 1U, MFLOAT >, Cell > &x, Field< T2, 1U, Cartesian< 1U, MFLOAT >, Cell > &w, Field< T1, 1U, Cartesian< 1U, MFLOAT >, Vert > &r)
Definition: Cartesian.hpp:3113
Offset
Definition: Offset.h:66
Tenzor
Definition: Tenzor.h:35
Vektor
Definition: Vektor.h:32
Field
Definition: Field.h:33
NDIndex::plugBase
NDIndex< Dim > plugBase(const NDIndex< D > &i) const
Definition: NDIndex.h:114
NDIndex::touches
bool touches(const NDIndex< Dim > &) const
Definition: NDIndexInlines.h:80
NDIndex::intersect
NDIndex< Dim > intersect(const NDIndex< Dim > &) const
Definition: NDIndexInlines.h:31
BareField::rightGuard
unsigned rightGuard(unsigned d) const
Definition: BareField.h:149
BareField::begin_if
iterator_if begin_if()
Definition: BareField.h:100
BareField::iterator_if
ac_id_larray::iterator iterator_if
Definition: BareField.h:92
BareField::end_if
iterator_if end_if()
Definition: BareField.h:101
BareField::getGuardCellSizes
const GuardCellSizes< Dim > & getGuardCellSizes() const
Definition: BareField.h:147
BareField::fillGuardCells
virtual void fillGuardCells(bool reallyFill=true) const
Definition: BareField.hpp:297
BareField::begin
iterator begin() const
Definition: BareField.h:234
BareField::leftGuard
unsigned leftGuard(unsigned d) const
Definition: BareField.h:148
LField
Definition: LField.h:58
LField::getOwned
const NDIndex< Dim > & getOwned() const
Definition: LField.h:99
LField::getAllocated
const NDIndex< Dim > & getAllocated() const
Definition: LField.h:98
LField::begin
const iterator & begin() const
Definition: LField.h:110
UniformCartesian::getVertexBelow
NDIndex< Dim > getVertexBelow(const Vektor< MFLOAT, Dim > &) const
Definition: UniformCartesian.hpp:1006
UniformCartesian::getSurfaceNormalFields
void getSurfaceNormalFields(Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > **) const
Definition: UniformCartesian.hpp:1176
UniformCartesian::getDeltaCellField
Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Vert > & getDeltaCellField(Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Vert > &) const
Definition: UniformCartesian.hpp:1145
UniformCartesian::set_origin
void set_origin(const Vektor< MFLOAT, Dim > &o)
Definition: UniformCartesian.hpp:412
UniformCartesian::getSurfaceNormals
Vektor< MFLOAT, Dim > * getSurfaceNormals(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:1158
UniformCartesian::getCellVolumeField
Field< MFLOAT, Dim, UniformCartesian< Dim, MFLOAT >, Cell > & getCellVolumeField(Field< MFLOAT, Dim, UniformCartesian< Dim, MFLOAT >, Cell > &) const
Definition: UniformCartesian.hpp:928
UniformCartesian::getVertexPosition
Vektor< MFLOAT, Dim > getVertexPosition(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:1025
UniformCartesian::getDeltaVertexField
Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > & getDeltaVertexField(Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > &) const
Definition: UniformCartesian.hpp:1118
UniformCartesian::get_meshSpacing
MFLOAT get_meshSpacing(unsigned d) const
Definition: UniformCartesian.hpp:461
UniformCartesian::getSurfaceNormalField
Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > & getSurfaceNormalField(Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > &, unsigned) const
Definition: UniformCartesian.hpp:1233
UniformCartesian::get_volume
MFLOAT get_volume() const
Definition: UniformCartesian.hpp:471
UniformCartesian::getDeltaVertex
Vektor< MFLOAT, Dim > getDeltaVertex(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:1105
UniformCartesian::getCellVolume
MFLOAT getCellVolume(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:915
UniformCartesian::getCellRangeVolume
MFLOAT getCellRangeVolume(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:962
UniformCartesian::getCellPositionField
Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > & getCellPositionField(Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Cell > &) const
Definition: UniformCartesian.hpp:1083
UniformCartesian::getVertexPositionField
Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Vert > & getVertexPositionField(Field< Vektor< MFLOAT, Dim >, Dim, UniformCartesian< Dim, MFLOAT >, Vert > &) const
Definition: UniformCartesian.hpp:1043
UniformCartesian::set_meshSpacing
void set_meshSpacing(MFLOAT *const del)
Definition: UniformCartesian.hpp:428
UniformCartesian::getNearestVertex
NDIndex< Dim > getNearestVertex(const Vektor< MFLOAT, Dim > &) const
Definition: UniformCartesian.hpp:987
UniformCartesian::getVertRangeVolume
MFLOAT getVertRangeVolume(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:938
UniformCartesian::getCellPosition
Vektor< MFLOAT, Dim > getCellPosition(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:1065
UniformCartesian::print
void print(std::ostream &)
Definition: UniformCartesian.hpp:883
UniformCartesian::get_origin
Vektor< MFLOAT, Dim > get_origin() const
Definition: UniformCartesian.hpp:420
UniformCartesian::getSurfaceNormal
Vektor< MFLOAT, Dim > getSurfaceNormal(const NDIndex< Dim > &, unsigned) const
Definition: UniformCartesian.hpp:1202
UniformCartesian::initialize
void initialize(const NDIndex< Dim > &ndi)
Definition: UniformCartesian.hpp:221
UniformCartesian::getDeltaCell
Vektor< MFLOAT, Dim > getDeltaCell(const NDIndex< Dim > &) const
Definition: UniformCartesian.hpp:1132
Index
Definition: Index.h:237
Index::stride
int stride() const
Definition: IndexInlines.h:121
Index::length
unsigned int length() const
Definition: IndexInlines.h:131
Index::first
int first() const
Definition: IndexInlines.h:116
Cell
Definition: Centering.h:32
Vert
Definition: Centering.h:43
Edge
Definition: Centering.h:50
OpUMeshExtrapolate::OpUMeshExtrapolate
OpUMeshExtrapolate(T &o, T &s)
Definition: UniformCartesian.hpp:485
Inform
Definition: Inform.h:42
Generated on Wed Aug 25 2021 16:40:59 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.1