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