00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "Utility/PAssert.h"
00031 #include "Utility/IpplInfo.h"
00032 #include "Field/BareField.h"
00033 #include "Field/BrickExpression.h"
00034 #include "Field/LField.h"
00035 #include "Field/Field.h"
00036 #include "Field/Assign.h"
00037 #include "Field/AssignDefs.h"
00038
00039
00040
00041
00042 template <unsigned Dim, class MFLOAT>
00043 void
00044 Cartesian<Dim,MFLOAT>::
00045 setup()
00046 {
00047 hasSpacingFields = false;
00048 }
00049
00050
00051
00052
00053 template <unsigned Dim, class MFLOAT>
00054 Cartesian<Dim,MFLOAT>::
00055 Cartesian(const NDIndex<Dim>& ndi)
00056 {
00057 int d,i;
00058 for (d=0; d<Dim; d++)
00059 gridSizes[d] = ndi[d].length();
00060 setup();
00061 for (d=0; d<Dim; d++) {
00062 MeshBC[2*d] = Reflective;
00063 MeshBC[2*d+1] = Reflective;
00064 origin(d) = ndi[d].first();
00065
00066 for (i=0; i < gridSizes[d]-1; i++) {
00067 (meshSpacing[d])[i] = ndi[d].stride();
00068 (meshPosition[d])[i] = MFLOAT(i);
00069 }
00070 (meshPosition[d])[gridSizes[d]-1] = MFLOAT(gridSizes[d]-1);
00071 }
00072 set_Dvc();
00073 }
00074
00075 template <unsigned Dim, class MFLOAT>
00076 Cartesian<Dim,MFLOAT>::
00077 Cartesian(const NDIndex<Dim>& ndi, MFLOAT** const delX)
00078 {
00079 unsigned int d;
00080 for (d=0; d<Dim; d++)
00081 gridSizes[d] = ndi[d].length();
00082 setup();
00083 for (d=0; d<Dim; d++) {
00084 MeshBC[2*d] = Reflective;
00085 MeshBC[2*d+1] = Reflective;
00086 origin(d) = ndi[d].first();
00087 }
00088 set_meshSpacing(delX);
00089 set_Dvc();
00090 }
00091
00092 template <unsigned Dim, class MFLOAT>
00093 Cartesian<Dim,MFLOAT>::
00094 Cartesian(const NDIndex<Dim>& ndi, MFLOAT** const delX,
00095 const Vektor<MFLOAT,Dim>& orig)
00096 {
00097 int d;
00098 for (d=0; d<Dim; d++)
00099 gridSizes[d] = ndi[d].length();
00100 setup();
00101 for (d=0; d<Dim; d++) {
00102 MeshBC[2*d] = Reflective;
00103 MeshBC[2*d+1] = Reflective;
00104 }
00105 set_origin(orig);
00106 set_meshSpacing(delX);
00107 set_Dvc();
00108 }
00109
00110 template <unsigned Dim, class MFLOAT>
00111 Cartesian<Dim,MFLOAT>::
00112 Cartesian(const NDIndex<Dim>& ndi, MFLOAT** const delX,
00113 const Vektor<MFLOAT,Dim>& orig, MeshBC_E* const mbc)
00114 {
00115 int d;
00116 for (d=0; d<Dim; d++)
00117 gridSizes[d] = ndi[d].length();
00118 setup();
00119 set_origin(orig);
00120 set_MeshBC(mbc);
00121 set_meshSpacing(delX);
00122 set_Dvc();
00123 }
00124
00125
00126
00127
00128
00129 template <unsigned Dim, class MFLOAT>
00130 Cartesian<Dim,MFLOAT>::
00131 Cartesian(const Index& I)
00132 {
00133 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00134 gridSizes[0] = I.length();
00135 setup();
00136 origin(0) = I.first();
00137 int i;
00138
00139 for (i=0; i < gridSizes[0]-1; i++) {
00140 (meshSpacing[0])[i] = I.stride();
00141 (meshPosition[0])[i] = MFLOAT(i);
00142 }
00143 (meshPosition[0])[gridSizes[0]-1] = MFLOAT(gridSizes[0]-1);
00144 for (int d=0; d<Dim; d++) {
00145 MeshBC[2*d] = Reflective;
00146 MeshBC[2*d+1] = Reflective;
00147 }
00148 set_Dvc();
00149 }
00150
00151 template <unsigned Dim, class MFLOAT>
00152 Cartesian<Dim,MFLOAT>::
00153 Cartesian(const Index& I, MFLOAT** const delX)
00154 {
00155 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00156 gridSizes[0] = I.length();
00157 setup();
00158 origin(0) = I.first();
00159 for (int d=0; d<Dim; d++) {
00160 MeshBC[2*d] = Reflective;
00161 MeshBC[2*d+1] = Reflective;
00162 }
00163 set_meshSpacing(delX);
00164 set_Dvc();
00165 }
00166
00167 template <unsigned Dim, class MFLOAT>
00168 Cartesian<Dim,MFLOAT>::
00169 Cartesian(const Index& I, MFLOAT** const delX,
00170 const Vektor<MFLOAT,Dim>& orig)
00171 {
00172 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00173 setup();
00174 gridSizes[0] = I.length();
00175 for (int d=0; d<Dim; d++) {
00176 MeshBC[2*d] = Reflective;
00177 MeshBC[2*d+1] = Reflective;
00178 }
00179 set_origin(orig);
00180 set_meshSpacing(delX);
00181 set_Dvc();
00182 }
00183
00184 template <unsigned Dim, class MFLOAT>
00185 Cartesian<Dim,MFLOAT>::
00186 Cartesian(const Index& I, MFLOAT** const delX,
00187 const Vektor<MFLOAT,Dim>& orig, MeshBC_E* const mbc)
00188 {
00189 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00190 setup();
00191 gridSizes[0] = I.length();
00192 set_origin(orig);
00193 set_MeshBC(mbc);
00194 set_meshSpacing(delX);
00195 set_Dvc();
00196 }
00197
00198
00199 template <unsigned Dim, class MFLOAT>
00200 Cartesian<Dim,MFLOAT>::
00201 Cartesian(const Index& I, const Index& J)
00202 {
00203 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00204 gridSizes[0] = I.length();
00205 gridSizes[1] = J.length();
00206 setup();
00207 origin(0) = I.first();
00208 origin(1) = J.first();
00209 int i;
00210
00211 for (i=0; i < gridSizes[0]-1; i++) {
00212 (meshSpacing[0])[i] = I.stride();
00213 (meshPosition[0])[i] = MFLOAT(i);
00214 }
00215 (meshPosition[0])[gridSizes[0]-1] = MFLOAT(gridSizes[0]-1);
00216 for (i=0; i < gridSizes[1]-1; i++) {
00217 (meshSpacing[1])[i] = J.stride();
00218 (meshPosition[1])[i] = MFLOAT(i);
00219 }
00220 (meshPosition[1])[gridSizes[1]-1] = MFLOAT(gridSizes[1]-1);
00221 for (int d=0; d<Dim; d++) {
00222 MeshBC[2*d] = Reflective;
00223 MeshBC[2*d+1] = Reflective;
00224 }
00225 set_Dvc();
00226 }
00227
00228 template <unsigned Dim, class MFLOAT>
00229 Cartesian<Dim,MFLOAT>::
00230 Cartesian(const Index& I, const Index& J, MFLOAT** const delX)
00231 {
00232 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00233 gridSizes[0] = I.length();
00234 gridSizes[1] = J.length();
00235 setup();
00236 origin(0) = I.first();
00237 origin(1) = J.first();
00238 for (int d=0; d<Dim; d++) {
00239 MeshBC[2*d] = Reflective;
00240 MeshBC[2*d+1] = Reflective;
00241 }
00242 set_meshSpacing(delX);
00243 set_Dvc();
00244 }
00245
00246 template <unsigned Dim, class MFLOAT>
00247 Cartesian<Dim,MFLOAT>::
00248 Cartesian(const Index& I, const Index& J, MFLOAT** const delX,
00249 const Vektor<MFLOAT,Dim>& orig)
00250 {
00251 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00252 gridSizes[0] = I.length();
00253 gridSizes[1] = J.length();
00254 setup();
00255 for (int d=0; d<Dim; d++) {
00256 MeshBC[2*d] = Reflective;
00257 MeshBC[2*d+1] = Reflective;
00258 }
00259 set_origin(orig);
00260 set_meshSpacing(delX);
00261 set_Dvc();
00262 }
00263
00264 template <unsigned Dim, class MFLOAT>
00265 Cartesian<Dim,MFLOAT>::
00266 Cartesian(const Index& I, const Index& J, MFLOAT** const delX,
00267 const Vektor<MFLOAT,Dim>& orig, MeshBC_E* const mbc)
00268 {
00269 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00270 gridSizes[0] = I.length();
00271 gridSizes[1] = J.length();
00272 setup();
00273 set_origin(orig);
00274 set_MeshBC(mbc);
00275 set_meshSpacing(delX);
00276 set_Dvc();
00277 }
00278
00279
00280 template <unsigned Dim, class MFLOAT>
00281 Cartesian<Dim,MFLOAT>::
00282 Cartesian(const Index& I, const Index& J, const Index& K)
00283 {
00284 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00285 gridSizes[0] = I.length();
00286 gridSizes[1] = J.length();
00287 gridSizes[2] = K.length();
00288
00289 setup();
00290 origin(0) = I.first();
00291 origin(1) = J.first();
00292 origin(2) = K.first();
00293 int i;
00294
00295 for (i=0; i < gridSizes[0]-1; i++) {
00296 (meshSpacing[0])[i] = I.stride();
00297 (meshPosition[0])[i] = MFLOAT(i);
00298 }
00299 (meshPosition[0])[gridSizes[0]-1] = MFLOAT(gridSizes[0]-1);
00300 for (i=0; i < gridSizes[1]-1; i++) {
00301 (meshSpacing[1])[i] = J.stride();
00302 (meshPosition[1])[i] = MFLOAT(i);
00303 }
00304 (meshPosition[1])[gridSizes[1]-1] = MFLOAT(gridSizes[1]-1);
00305 for (i=0; i < gridSizes[2]-1; i++) {
00306 (meshSpacing[2])[i] = K.stride();
00307 (meshPosition[2])[i] = MFLOAT(i);
00308 }
00309 (meshPosition[2])[gridSizes[2]-1] = MFLOAT(gridSizes[2]-1);
00310
00311 for (int d=0; d<Dim; d++) {
00312 MeshBC[2*d] = Reflective;
00313 MeshBC[2*d+1] = Reflective;
00314 }
00315 set_Dvc();
00316 }
00317
00318 template <unsigned Dim, class MFLOAT>
00319 Cartesian<Dim,MFLOAT>::
00320 Cartesian(const Index& I, const Index& J, const Index& K,
00321 MFLOAT** const delX)
00322 {
00323 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00324 gridSizes[0] = I.length();
00325 gridSizes[1] = J.length();
00326 gridSizes[2] = K.length();
00327 setup();
00328 origin(0) = I.first();
00329 origin(1) = J.first();
00330 origin(2) = K.first();
00331 for (int d=0; d<Dim; d++) {
00332 MeshBC[2*d] = Reflective;
00333 MeshBC[2*d+1] = Reflective;
00334 }
00335 set_meshSpacing(delX);
00336 set_Dvc();
00337 }
00338
00339 template <unsigned Dim, class MFLOAT>
00340 Cartesian<Dim,MFLOAT>::
00341 Cartesian(const Index& I, const Index& J, const Index& K,
00342 MFLOAT** const delX, const Vektor<MFLOAT,Dim>& orig)
00343 {
00344 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00345 gridSizes[0] = I.length();
00346 gridSizes[1] = J.length();
00347 gridSizes[2] = K.length();
00348 setup();
00349 for (int d=0; d<Dim; d++) {
00350 MeshBC[2*d] = Reflective;
00351 MeshBC[2*d+1] = Reflective;
00352 }
00353 set_origin(orig);
00354 set_meshSpacing(delX);
00355 set_Dvc();
00356 }
00357
00358 template <unsigned Dim, class MFLOAT>
00359 Cartesian<Dim,MFLOAT>::
00360 Cartesian(const Index& I, const Index& J, const Index& K,
00361 MFLOAT** const delX, const Vektor<MFLOAT,Dim>& orig,
00362 MeshBC_E* const mbc)
00363 {
00364 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00365 gridSizes[0] = I.length();
00366 gridSizes[1] = J.length();
00367 gridSizes[2] = K.length();
00368 setup();
00369 set_origin(orig);
00370 set_MeshBC(mbc);
00371 set_meshSpacing(delX);
00372 set_Dvc();
00373 }
00374
00375
00376
00377
00378 template <unsigned Dim, class MFLOAT>
00379 void
00380 Cartesian<Dim,MFLOAT>::
00381 initialize(const NDIndex<Dim>& ndi)
00382 {
00383 int d,i;
00384 for (d=0; d<Dim; d++)
00385 gridSizes[d] = ndi[d].length();
00386 setup();
00387 for (d=0; d<Dim; d++) {
00388 MeshBC[2*d] = Reflective;
00389 MeshBC[2*d+1] = Reflective;
00390 origin(d) = ndi[d].first();
00391
00392 for (i=0; i < gridSizes[d]-1; i++) {
00393 (meshSpacing[d])[i] = ndi[d].stride();
00394 (meshPosition[d])[i] = MFLOAT(i);
00395 }
00396 (meshPosition[d])[gridSizes[d]-1] = MFLOAT(gridSizes[d]-1);
00397 }
00398 set_Dvc();
00399 }
00400
00401 template <unsigned Dim, class MFLOAT>
00402 void
00403 Cartesian<Dim,MFLOAT>::
00404 initialize(const NDIndex<Dim>& ndi, MFLOAT** const delX)
00405 {
00406 unsigned int d;
00407 for (d=0; d<Dim; d++)
00408 gridSizes[d] = ndi[d].length();
00409 setup();
00410 for (d=0; d<Dim; d++) {
00411 MeshBC[2*d] = Reflective;
00412 MeshBC[2*d+1] = Reflective;
00413 origin(d) = ndi[d].first();
00414 }
00415 set_meshSpacing(delX);
00416 set_Dvc();
00417 }
00418
00419 template <unsigned Dim, class MFLOAT>
00420 void
00421 Cartesian<Dim,MFLOAT>::
00422 initialize(const NDIndex<Dim>& ndi, MFLOAT** const delX,
00423 const Vektor<MFLOAT,Dim>& orig)
00424 {
00425 int d;
00426 for (d=0; d<Dim; d++)
00427 gridSizes[d] = ndi[d].length();
00428 setup();
00429 for (d=0; d<Dim; d++) {
00430 MeshBC[2*d] = Reflective;
00431 MeshBC[2*d+1] = Reflective;
00432 }
00433 set_origin(orig);
00434 set_meshSpacing(delX);
00435 set_Dvc();
00436 }
00437
00438 template <unsigned Dim, class MFLOAT>
00439 void
00440 Cartesian<Dim,MFLOAT>::
00441 initialize(const NDIndex<Dim>& ndi, MFLOAT** const delX,
00442 const Vektor<MFLOAT,Dim>& orig, MeshBC_E* const mbc)
00443 {
00444 int d;
00445 for (d=0; d<Dim; d++)
00446 gridSizes[d] = ndi[d].length();
00447 setup();
00448 set_origin(orig);
00449 set_MeshBC(mbc);
00450 set_meshSpacing(delX);
00451 set_Dvc();
00452 }
00453
00454
00455
00456
00457
00458 template <unsigned Dim, class MFLOAT>
00459 void
00460 Cartesian<Dim,MFLOAT>::
00461 initialize(const Index& I)
00462 {
00463 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00464 gridSizes[0] = I.length();
00465 setup();
00466 origin(0) = I.first();
00467 int i;
00468
00469 for (i=0; i < gridSizes[0]-1; i++) {
00470 (meshSpacing[0])[i] = I.stride();
00471 (meshPosition[0])[i] = MFLOAT(i);
00472 }
00473 (meshPosition[0])[gridSizes[0]-1] = MFLOAT(gridSizes[0]-1);
00474 for (int d=0; d<Dim; d++) {
00475 MeshBC[2*d] = Reflective;
00476 MeshBC[2*d+1] = Reflective;
00477 }
00478 set_Dvc();
00479 }
00480
00481 template <unsigned Dim, class MFLOAT>
00482 void
00483 Cartesian<Dim,MFLOAT>::
00484 initialize(const Index& I, MFLOAT** const delX)
00485 {
00486 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00487 gridSizes[0] = I.length();
00488 setup();
00489 origin(0) = I.first();
00490 for (int d=0; d<Dim; d++) {
00491 MeshBC[2*d] = Reflective;
00492 MeshBC[2*d+1] = Reflective;
00493 }
00494 set_meshSpacing(delX);
00495 set_Dvc();
00496 }
00497
00498 template <unsigned Dim, class MFLOAT>
00499 void
00500 Cartesian<Dim,MFLOAT>::
00501 initialize(const Index& I, MFLOAT** const delX,
00502 const Vektor<MFLOAT,Dim>& orig)
00503 {
00504 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00505 setup();
00506 gridSizes[0] = I.length();
00507 for (int d=0; d<Dim; d++) {
00508 MeshBC[2*d] = Reflective;
00509 MeshBC[2*d+1] = Reflective;
00510 }
00511 set_origin(orig);
00512 set_meshSpacing(delX);
00513 set_Dvc();
00514 }
00515
00516 template <unsigned Dim, class MFLOAT>
00517 void
00518 Cartesian<Dim,MFLOAT>::
00519 initialize(const Index& I, MFLOAT** const delX,
00520 const Vektor<MFLOAT,Dim>& orig, MeshBC_E* const mbc)
00521 {
00522 PInsist(Dim==1,"Number of Index arguments does not match mesh dimension!!");
00523 setup();
00524 gridSizes[0] = I.length();
00525 set_origin(orig);
00526 set_MeshBC(mbc);
00527 set_meshSpacing(delX);
00528 set_Dvc();
00529 }
00530
00531
00532 template <unsigned Dim, class MFLOAT>
00533 void
00534 Cartesian<Dim,MFLOAT>::
00535 initialize(const Index& I, const Index& J)
00536 {
00537 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00538 gridSizes[0] = I.length();
00539 gridSizes[1] = J.length();
00540 setup();
00541 origin(0) = I.first();
00542 origin(1) = J.first();
00543 int i;
00544
00545 for (i=0; i < gridSizes[0]-1; i++) {
00546 (meshSpacing[0])[i] = I.stride();
00547 (meshPosition[0])[i] = MFLOAT(i);
00548 }
00549 (meshPosition[0])[gridSizes[0]-1] = MFLOAT(gridSizes[0]-1);
00550 for (i=0; i < gridSizes[1]-1; i++) {
00551 (meshSpacing[1])[i] = J.stride();
00552 (meshPosition[1])[i] = MFLOAT(i);
00553 }
00554 (meshPosition[1])[gridSizes[1]-1] = MFLOAT(gridSizes[1]-1);
00555 for (int d=0; d<Dim; d++) {
00556 MeshBC[2*d] = Reflective;
00557 MeshBC[2*d+1] = Reflective;
00558 }
00559 set_Dvc();
00560 }
00561
00562 template <unsigned Dim, class MFLOAT>
00563 void
00564 Cartesian<Dim,MFLOAT>::
00565 initialize(const Index& I, const Index& J, MFLOAT** const delX)
00566 {
00567 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00568 gridSizes[0] = I.length();
00569 gridSizes[1] = J.length();
00570 setup();
00571 origin(0) = I.first();
00572 origin(1) = J.first();
00573 for (int d=0; d<Dim; d++) {
00574 MeshBC[2*d] = Reflective;
00575 MeshBC[2*d+1] = Reflective;
00576 }
00577 set_meshSpacing(delX);
00578 set_Dvc();
00579 }
00580
00581 template <unsigned Dim, class MFLOAT>
00582 void
00583 Cartesian<Dim,MFLOAT>::
00584 initialize(const Index& I, const Index& J, MFLOAT** const delX,
00585 const Vektor<MFLOAT,Dim>& orig)
00586 {
00587 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00588 gridSizes[0] = I.length();
00589 gridSizes[1] = J.length();
00590 setup();
00591 for (int d=0; d<Dim; d++) {
00592 MeshBC[2*d] = Reflective;
00593 MeshBC[2*d+1] = Reflective;
00594 }
00595 set_origin(orig);
00596 set_meshSpacing(delX);
00597 set_Dvc();
00598 }
00599
00600 template <unsigned Dim, class MFLOAT>
00601 void
00602 Cartesian<Dim,MFLOAT>::
00603 initialize(const Index& I, const Index& J, MFLOAT** const delX,
00604 const Vektor<MFLOAT,Dim>& orig, MeshBC_E* const mbc)
00605 {
00606 PInsist(Dim==2,"Number of Index arguments does not match mesh dimension!!");
00607 gridSizes[0] = I.length();
00608 gridSizes[1] = J.length();
00609 setup();
00610 set_origin(orig);
00611 set_MeshBC(mbc);
00612 set_meshSpacing(delX);
00613 set_Dvc();
00614 }
00615
00616
00617 template <unsigned Dim, class MFLOAT>
00618 void
00619 Cartesian<Dim,MFLOAT>::
00620 initialize(const Index& I, const Index& J, const Index& K)
00621 {
00622 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00623 gridSizes[0] = I.length();
00624 gridSizes[1] = J.length();
00625 gridSizes[2] = K.length();
00626
00627 setup();
00628 origin(0) = I.first();
00629 origin(1) = J.first();
00630 origin(2) = K.first();
00631 int i;
00632
00633 for (i=0; i < gridSizes[0]-1; i++) {
00634 (meshSpacing[0])[i] = I.stride();
00635 (meshPosition[0])[i] = MFLOAT(i);
00636 }
00637 (meshPosition[0])[gridSizes[0]-1] = MFLOAT(gridSizes[0]-1);
00638 for (i=0; i < gridSizes[1]-1; i++) {
00639 (meshSpacing[1])[i] = J.stride();
00640 (meshPosition[1])[i] = MFLOAT(i);
00641 }
00642 (meshPosition[1])[gridSizes[1]-1] = MFLOAT(gridSizes[1]-1);
00643 for (i=0; i < gridSizes[2]-1; i++) {
00644 (meshSpacing[2])[i] = K.stride();
00645 (meshPosition[2])[i] = MFLOAT(i);
00646 }
00647 (meshPosition[2])[gridSizes[2]-1] = MFLOAT(gridSizes[2]-1);
00648
00649 for (int d=0; d<Dim; d++) {
00650 MeshBC[2*d] = Reflective;
00651 MeshBC[2*d+1] = Reflective;
00652 }
00653 set_Dvc();
00654 }
00655
00656 template <unsigned Dim, class MFLOAT>
00657 void
00658 Cartesian<Dim,MFLOAT>::
00659 initialize(const Index& I, const Index& J, const Index& K,
00660 MFLOAT** const delX)
00661 {
00662 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00663 gridSizes[0] = I.length();
00664 gridSizes[1] = J.length();
00665 gridSizes[2] = K.length();
00666 setup();
00667 origin(0) = I.first();
00668 origin(1) = J.first();
00669 origin(2) = K.first();
00670 for (int d=0; d<Dim; d++) {
00671 MeshBC[2*d] = Reflective;
00672 MeshBC[2*d+1] = Reflective;
00673 }
00674 set_meshSpacing(delX);
00675 set_Dvc();
00676 }
00677
00678 template <unsigned Dim, class MFLOAT>
00679 void
00680 Cartesian<Dim,MFLOAT>::
00681 initialize(const Index& I, const Index& J, const Index& K,
00682 MFLOAT** const delX, const Vektor<MFLOAT,Dim>& orig)
00683 {
00684 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00685 gridSizes[0] = I.length();
00686 gridSizes[1] = J.length();
00687 gridSizes[2] = K.length();
00688 setup();
00689 for (int d=0; d<Dim; d++) {
00690 MeshBC[2*d] = Reflective;
00691 MeshBC[2*d+1] = Reflective;
00692 }
00693 set_origin(orig);
00694 set_meshSpacing(delX);
00695 set_Dvc();
00696 }
00697
00698 template <unsigned Dim, class MFLOAT>
00699 void
00700 Cartesian<Dim,MFLOAT>::
00701 initialize(const Index& I, const Index& J, const Index& K,
00702 MFLOAT** const delX, const Vektor<MFLOAT,Dim>& orig,
00703 MeshBC_E* const mbc)
00704 {
00705 PInsist(Dim==3,"Number of Index arguments does not match mesh dimension!!");
00706 gridSizes[0] = I.length();
00707 gridSizes[1] = J.length();
00708 gridSizes[2] = K.length();
00709 setup();
00710 set_origin(orig);
00711 set_MeshBC(mbc);
00712 set_meshSpacing(delX);
00713 set_Dvc();
00714 }
00715
00716
00717
00718
00719
00720 template<unsigned Dim, class MFLOAT>
00721 void Cartesian<Dim,MFLOAT>::
00722 set_origin(const Vektor<MFLOAT,Dim>& o)
00723 {
00724 origin = o;
00725 for (unsigned d=0; d<Dim; ++d) {
00726 (meshPosition[d])[0] = o(d);
00727 for (unsigned vert=1; vert<gridSizes[d]; ++vert) {
00728 (meshPosition[d])[vert] = (meshPosition[d])[vert-1] +
00729 (meshSpacing[d])[vert-1];
00730 }
00731 }
00732
00733 for (unsigned face=0; face < 2*Dim; ++face) updateMeshSpacingGuards(face);
00734 this->notifyOfChange();
00735 }
00736
00737 template<unsigned Dim, class MFLOAT>
00738 Vektor<MFLOAT,Dim> Cartesian<Dim,MFLOAT>::
00739 get_origin() const
00740 {
00741 return origin;
00742 }
00743
00744
00745 template<unsigned Dim, class MFLOAT>
00746 void Cartesian<Dim,MFLOAT>::
00747 set_meshSpacing(MFLOAT** const del)
00748 {
00749 unsigned d, cell, face;
00750
00751 for (d=0;d<Dim;++d) {
00752 (meshPosition[d])[0] = origin(d);
00753 for (cell=0; cell < gridSizes[d]-1; cell++) {
00754 (meshSpacing[d])[cell] = del[d][cell];
00755 (meshPosition[d])[cell+1] = (meshPosition[d])[cell] + del[d][cell];
00756 }
00757 }
00758
00759 for (face=0; face < 2*Dim; ++face) updateMeshSpacingGuards(face);
00760
00761 if (hasSpacingFields) storeSpacingFields();
00762 this->notifyOfChange();
00763 }
00764
00765
00766 template<unsigned Dim, class MFLOAT>
00767 void Cartesian<Dim,MFLOAT>::
00768 set_Dvc()
00769 {
00770 unsigned d;
00771 MFLOAT coef = 1.0;
00772 for (d=1;d<Dim;++d) coef *= 0.5;
00773
00774 for (d=0;d<Dim;++d) {
00775 MFLOAT dvc = coef;
00776 for (unsigned b=0; b<(1<<Dim); ++b) {
00777 int s = ( b&(1<<d) ) ? 1 : -1;
00778 Dvc[b](d) = s*dvc;
00779 }
00780 }
00781 }
00782
00783
00784 template<unsigned Dim, class MFLOAT>
00785 void Cartesian<Dim,MFLOAT>::
00786 get_meshSpacing(unsigned d, MFLOAT* spacings) const
00787 {
00788 PAssert(d<Dim);
00789 for (int cell=0; cell < gridSizes[d]-1; cell++)
00790 spacings[cell] = (*(meshSpacing[d].find(cell))).second;
00791 return;
00792 }
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00805
00806
00807
00808
00809
00810 template<class T>
00811 struct OpMeshPeriodic
00812 {
00813 #ifdef IPPL_PURIFY
00814 OpMeshPeriodic() {}
00815 OpMeshPeriodic(const OpMeshPeriodic<T> &) {}
00816 OpMeshPeriodic<T>& operator=(const OpMeshPeriodic<T> &) { return *this; }
00817 #endif
00818 };
00819 template<class T>
00820 inline void PETE_apply(OpMeshPeriodic<T> e, T& a, T b) { a = b; }
00821
00822
00823 template<class T>
00824 struct OpMeshExtrapolate
00825 {
00826 OpMeshExtrapolate(T& o, T& s) : Offset(o), Slope(s) {}
00827 T Offset, Slope;
00828 };
00829
00830
00831 template<class T>
00832 void PETE_apply(OpMeshExtrapolate<T> e, T& a, T b)
00833 {
00834 a = b*e.Slope+e.Offset;
00835 }
00836
00838
00839
00840
00841
00842 template<unsigned Dim, class MFLOAT>
00843 void Cartesian<Dim,MFLOAT>::
00844 storeSpacingFields()
00845 {
00846
00847 e_dim_tag et[Dim];
00848 for (int d=0; d<Dim; d++) et[d] = PARALLEL;
00849 storeSpacingFields(et, -1);
00850 }
00851
00852 template<unsigned Dim, class MFLOAT>
00853 void Cartesian<Dim,MFLOAT>::
00854 storeSpacingFields(e_dim_tag p1, int vnodes)
00855 {
00856 e_dim_tag et[1];
00857 et[0] = p1;
00858 storeSpacingFields(et, vnodes);
00859 }
00860
00861 template<unsigned Dim, class MFLOAT>
00862 void Cartesian<Dim,MFLOAT>::
00863 storeSpacingFields(e_dim_tag p1, e_dim_tag p2, int vnodes)
00864 {
00865 e_dim_tag et[2];
00866 et[0] = p1;
00867 et[1] = p2;
00868 storeSpacingFields(et, vnodes);
00869 }
00870
00871 template<unsigned Dim, class MFLOAT>
00872 void Cartesian<Dim,MFLOAT>::
00873 storeSpacingFields(e_dim_tag p1, e_dim_tag p2, e_dim_tag p3, int vnodes)
00874 {
00875 e_dim_tag et[3];
00876 et[0] = p1;
00877 et[1] = p2;
00878 et[2] = p3;
00879 storeSpacingFields(et, vnodes);
00880 }
00881
00882
00883
00884 template<unsigned Dim, class MFLOAT>
00885 void Cartesian<Dim,MFLOAT>::
00886 storeSpacingFields(e_dim_tag* et, int vnodes)
00887 {
00888 int d;
00889 int currentLocation[Dim];
00890 NDIndex<Dim> cells, verts;
00891 for (d=0; d<Dim; d++) {
00892 cells[d] = Index(gridSizes[d]-1);
00893 verts[d] = Index(gridSizes[d]);
00894 }
00895 if (!hasSpacingFields) {
00896
00897 FlCell = new FieldLayout<Dim>(cells, et, vnodes);
00898
00899 VertSpacings =
00900 new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlCell,GuardCellSizes<Dim>(1));
00901 FlVert = new FieldLayout<Dim>(verts, et, vnodes);
00902
00903 CellSpacings =
00904 new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlVert,GuardCellSizes<Dim>(1));
00905 }
00906
00907 BareField<Vektor<MFLOAT,Dim>,Dim>& vertSpacings = *VertSpacings;
00908 Vektor<MFLOAT,Dim> vertexSpacing;
00909 vertSpacings.Uncompress();
00910 typename BareField<Vektor<MFLOAT,Dim>,Dim>::iterator cfi,
00911 cfi_end = vertSpacings.end();
00912 for (cfi = vertSpacings.begin(); cfi != cfi_end; ++cfi) {
00913 cfi.GetCurrentLocation(currentLocation);
00914 for (d=0; d<Dim; d++)
00915 vertexSpacing(d) = (*(meshSpacing[d].find(currentLocation[d]))).second;
00916 *cfi = vertexSpacing;
00917 }
00918
00919 BareField<Vektor<MFLOAT,Dim>,Dim>& cellSpacings = *CellSpacings;
00920 Vektor<MFLOAT,Dim> cellSpacing;
00921 cellSpacings.Uncompress();
00922 typename BareField<Vektor<MFLOAT,Dim>,Dim>::iterator vfi,
00923 vfi_end = cellSpacings.end();
00924 for (vfi = cellSpacings.begin(); vfi != vfi_end; ++vfi) {
00925 vfi.GetCurrentLocation(currentLocation);
00926 for (d=0; d<Dim; d++)
00927 cellSpacing(d) = 0.5 * ((meshSpacing[d])[currentLocation[d]] +
00928 (meshSpacing[d])[currentLocation[d]-1]);
00929 *vfi = cellSpacing;
00930 }
00931
00932
00933
00934
00935
00936 cellSpacings.fillGuardCells();
00937 vertSpacings.fillGuardCells();
00938
00939
00940
00941
00942 Vektor<MFLOAT,Dim> v0,v1; v0 = 0.0; v1 = 1.0;
00943 int face;
00944 typedef Vektor<MFLOAT,Dim> T;
00945 typename BareField<T,Dim>::iterator_if cfill_i;
00946 typename BareField<T,Dim>::iterator_if vfill_i;
00947 int coffset, voffset;
00948 MeshBC_E bct;
00949
00950 for (face=0; face < 2*Dim; face++) {
00951
00952 NDIndex<Dim> cSlab = AddGuardCells(verts,cellSpacings.getGuardCellSizes());
00953 NDIndex<Dim> vSlab = AddGuardCells(cells,vertSpacings.getGuardCellSizes());
00954
00955 d = face/2;
00956
00957
00958
00959
00960 if ( face & 1 ) {
00961 cSlab[d] = Index(verts[d].max() + 1,
00962 verts[d].max() + cellSpacings.rightGuard(d));
00963 vSlab[d] = Index(cells[d].max() + 1,
00964 cells[d].max() + vertSpacings.rightGuard(d));
00965 } else {
00966 cSlab[d] = Index(verts[d].min() - cellSpacings.leftGuard(d),
00967 verts[d].min() - 1);
00968 vSlab[d] = Index(cells[d].min() - vertSpacings.leftGuard(d),
00969 cells[d].min() - 1);
00970 }
00971
00972 switch (MeshBC[face]) {
00973 case Periodic:
00974 bct = Periodic;
00975 if ( face & 1 ) {
00976 coffset = -verts[d].length();
00977 voffset = -cells[d].length();
00978 } else {
00979 coffset = verts[d].length();
00980 voffset = cells[d].length();
00981 }
00982 break;
00983 case Reflective:
00984 bct = Reflective;
00985 if ( face & 1 ) {
00986 coffset = 2*verts[d].max() + 1;
00987 voffset = 2*cells[d].max() + 1 - 1;
00988 } else {
00989 coffset = 2*verts[d].min() - 1;
00990 voffset = 2*cells[d].min() - 1 + 1;
00991 }
00992 break;
00993 case NoBC:
00994 bct = NoBC;
00995 if ( face & 1 ) {
00996 coffset = 2*verts[d].max() + 1;
00997 voffset = 2*cells[d].max() + 1 - 1;
00998 } else {
00999 coffset = 2*verts[d].min() - 1;
01000 voffset = 2*cells[d].min() - 1 + 1;
01001 }
01002 break;
01003 default:
01004 ERRORMSG("Cartesian::storeSpacingFields(): unknown MeshBC type" << endl);
01005 break;
01006 }
01007
01008
01009
01010 for (cfill_i=cellSpacings.begin_if();
01011 cfill_i!=cellSpacings.end_if(); ++cfill_i)
01012 {
01013
01014
01015 LField<T,Dim> &fill = *(*cfill_i).second;
01016
01017 const NDIndex<Dim> &fill_alloc = fill.getAllocated();
01018
01019
01020
01021
01022 if ( cSlab.touches( fill_alloc ) )
01023 {
01024
01025 NDIndex<Dim> dest = cSlab.intersect( fill_alloc );
01026
01027
01028
01029
01030
01031
01032 NDIndex<Dim> src = dest;
01033
01034
01035 src[d] = coffset - src[d];
01036
01037
01038
01039 typename BareField<T,Dim>::iterator_if from_i;
01040 for (from_i=cellSpacings.begin_if();
01041 from_i!=cellSpacings.end_if(); ++from_i)
01042 {
01043
01044 LField<T,Dim> &from = *(*from_i).second;
01045 const NDIndex<Dim> &from_owned = from.getOwned();
01046 const NDIndex<Dim> &from_alloc = from.getAllocated();
01047
01048 if ( src.touches( from_owned ) )
01049 {
01050 NDIndex<Dim> from_it = src.intersect( from_alloc );
01051 NDIndex<Dim> cfill_it = dest.plugBase( from_it );
01052
01053 typedef typename LField<T,Dim>::iterator LFI;
01054 LFI lhs = fill.begin(cfill_it);
01055 LFI rhs = from.begin(from_it);
01056
01057 if (bct == Periodic) {
01058 BrickExpression<Dim,LFI,LFI,OpMeshPeriodic<T> >
01059 (lhs,rhs,OpMeshPeriodic<T>()).apply();
01060 } else {
01061 if (bct == Reflective) {
01062 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01063 (lhs,rhs,OpMeshExtrapolate<T>(v0,v1)).apply();
01064 } else {
01065 if (bct == NoBC) {
01066 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01067 (lhs,rhs,OpMeshExtrapolate<T>(v0,v0)).apply();
01068 }
01069 }
01070 }
01071 }
01072 }
01073 }
01074 }
01075
01076 for (vfill_i=vertSpacings.begin_if();
01077 vfill_i!=vertSpacings.end_if(); ++vfill_i)
01078 {
01079
01080
01081 LField<T,Dim> &fill = *(*vfill_i).second;
01082
01083 const NDIndex<Dim> &fill_alloc = fill.getAllocated();
01084
01085
01086
01087
01088 if ( vSlab.touches( fill_alloc ) )
01089 {
01090
01091 NDIndex<Dim> dest = vSlab.intersect( fill_alloc );
01092
01093
01094
01095
01096
01097
01098 NDIndex<Dim> src = dest;
01099
01100
01101 src[d] = voffset - src[d];
01102
01103
01104
01105 typename BareField<T,Dim>::iterator_if from_i;
01106 for (from_i=vertSpacings.begin_if();
01107 from_i!=vertSpacings.end_if(); ++from_i)
01108 {
01109
01110 LField<T,Dim> &from = *(*from_i).second;
01111 const NDIndex<Dim> &from_owned = from.getOwned();
01112 const NDIndex<Dim> &from_alloc = from.getAllocated();
01113
01114 if ( src.touches( from_owned ) )
01115 {
01116 NDIndex<Dim> from_it = src.intersect( from_alloc );
01117 NDIndex<Dim> vfill_it = dest.plugBase( from_it );
01118
01119 typedef typename LField<T,Dim>::iterator LFI;
01120 LFI lhs = fill.begin(vfill_it);
01121 LFI rhs = from.begin(from_it);
01122
01123 if (bct == Periodic) {
01124 BrickExpression<Dim,LFI,LFI,OpMeshPeriodic<T> >
01125 (lhs,rhs,OpMeshPeriodic<T>()).apply();
01126 } else {
01127 if (bct == Reflective) {
01128 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01129 (lhs,rhs,OpMeshExtrapolate<T>(v0,v1)).apply();
01130 } else {
01131 if (bct == NoBC) {
01132 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01133 (lhs,rhs,OpMeshExtrapolate<T>(v0,v0)).apply();
01134 }
01135 }
01136 }
01137 }
01138 }
01139 }
01140 }
01141
01142 }
01143
01144 hasSpacingFields = true;
01145 }
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156 template<unsigned Dim, class MFLOAT>
01157 void Cartesian<Dim,MFLOAT>::
01158 storeSpacingFields(e_dim_tag p1,
01159 unsigned vnodes1,
01160 bool recurse,
01161 int vnodes) {
01162 e_dim_tag et[1];
01163 et[0] = p1;
01164 unsigned vnodesPerDirection[Dim];
01165 vnodesPerDirection[0] = vnodes1;
01166 storeSpacingFields(et, vnodesPerDirection, recurse, vnodes);
01167 }
01168
01169 template<unsigned Dim, class MFLOAT>
01170 void Cartesian<Dim,MFLOAT>::
01171 storeSpacingFields(e_dim_tag p1, e_dim_tag p2,
01172 unsigned vnodes1, unsigned vnodes2,
01173 bool recurse,int vnodes) {
01174 e_dim_tag et[2];
01175 et[0] = p1;
01176 et[1] = p2;
01177 unsigned vnodesPerDirection[Dim];
01178 vnodesPerDirection[0] = vnodes1;
01179 vnodesPerDirection[1] = vnodes2;
01180 storeSpacingFields(et, vnodesPerDirection, recurse, vnodes);
01181 }
01182
01183 template<unsigned Dim, class MFLOAT>
01184 void Cartesian<Dim,MFLOAT>::
01185 storeSpacingFields(e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
01186 unsigned vnodes1, unsigned vnodes2, unsigned vnodes3,
01187 bool recurse, int vnodes) {
01188 e_dim_tag et[3];
01189 et[0] = p1;
01190 et[1] = p2;
01191 et[2] = p3;
01192 unsigned vnodesPerDirection[Dim];
01193 vnodesPerDirection[0] = vnodes1;
01194 vnodesPerDirection[1] = vnodes2;
01195 vnodesPerDirection[2] = vnodes3;
01196 storeSpacingFields(et, vnodesPerDirection, recurse, vnodes);
01197 }
01198
01199
01200
01201
01202
01203 template<unsigned Dim, class MFLOAT>
01204 void Cartesian<Dim,MFLOAT>::
01205 storeSpacingFields(e_dim_tag *p,
01206 unsigned* vnodesPerDirection,
01207 bool recurse, int vnodes) {
01208 int d;
01209 int currentLocation[Dim];
01210 NDIndex<Dim> cells, verts;
01211 for (d=0; d<Dim; d++) {
01212 cells[d] = Index(gridSizes[d]-1);
01213 verts[d] = Index(gridSizes[d]);
01214 }
01215 if (!hasSpacingFields) {
01216
01217 FlCell =
01218 new FieldLayout<Dim>(cells, p, vnodesPerDirection, recurse, vnodes);
01219
01220 VertSpacings =
01221 new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlCell,GuardCellSizes<Dim>(1));
01222 FlVert =
01223 new FieldLayout<Dim>(verts, p, vnodesPerDirection, recurse, vnodes);
01224
01225 CellSpacings =
01226 new BareField<Vektor<MFLOAT,Dim>,Dim>(*FlVert,GuardCellSizes<Dim>(1));
01227 }
01228
01229 BareField<Vektor<MFLOAT,Dim>,Dim>& vertSpacings = *VertSpacings;
01230 Vektor<MFLOAT,Dim> vertexSpacing;
01231 vertSpacings.Uncompress();
01232 typename BareField<Vektor<MFLOAT,Dim>,Dim>::iterator cfi,
01233 cfi_end = vertSpacings.end();
01234 for (cfi = vertSpacings.begin(); cfi != cfi_end; ++cfi) {
01235 cfi.GetCurrentLocation(currentLocation);
01236 for (d=0; d<Dim; d++)
01237 vertexSpacing(d) = (*(meshSpacing[d].find(currentLocation[d]))).second;
01238 *cfi = vertexSpacing;
01239 }
01240
01241 BareField<Vektor<MFLOAT,Dim>,Dim>& cellSpacings = *CellSpacings;
01242 Vektor<MFLOAT,Dim> cellSpacing;
01243 cellSpacings.Uncompress();
01244 typename BareField<Vektor<MFLOAT,Dim>,Dim>::iterator vfi,
01245 vfi_end = cellSpacings.end();
01246 for (vfi = cellSpacings.begin(); vfi != vfi_end; ++vfi) {
01247 vfi.GetCurrentLocation(currentLocation);
01248 for (d=0; d<Dim; d++)
01249 cellSpacing(d) = 0.5 * ((meshSpacing[d])[currentLocation[d]] +
01250 (meshSpacing[d])[currentLocation[d]-1]);
01251 *vfi = cellSpacing;
01252 }
01253
01254
01255
01256
01257
01258 cellSpacings.fillGuardCells();
01259 vertSpacings.fillGuardCells();
01260
01261
01262
01263
01264 Vektor<MFLOAT,Dim> v0,v1; v0 = 0.0; v1 = 1.0;
01265 int face;
01266 typedef Vektor<MFLOAT,Dim> T;
01267 typename BareField<T,Dim>::iterator_if cfill_i;
01268 typename BareField<T,Dim>::iterator_if vfill_i;
01269 int coffset, voffset;
01270 MeshBC_E bct;
01271
01272 for (face=0; face < 2*Dim; face++) {
01273
01274 NDIndex<Dim> cSlab = AddGuardCells(verts,cellSpacings.getGuardCellSizes());
01275 NDIndex<Dim> vSlab = AddGuardCells(cells,vertSpacings.getGuardCellSizes());
01276
01277 d = face/2;
01278
01279
01280
01281
01282 if ( face & 1 ) {
01283 cSlab[d] = Index(verts[d].max() + 1,
01284 verts[d].max() + cellSpacings.rightGuard(d));
01285 vSlab[d] = Index(cells[d].max() + 1,
01286 cells[d].max() + vertSpacings.rightGuard(d));
01287 } else {
01288 cSlab[d] = Index(verts[d].min() - cellSpacings.leftGuard(d),
01289 verts[d].min() - 1);
01290 vSlab[d] = Index(cells[d].min() - vertSpacings.leftGuard(d),
01291 cells[d].min() - 1);
01292 }
01293
01294 switch (MeshBC[face]) {
01295 case Periodic:
01296 bct = Periodic;
01297 if ( face & 1 ) {
01298 coffset = -verts[d].length();
01299 voffset = -cells[d].length();
01300 } else {
01301 coffset = verts[d].length();
01302 voffset = cells[d].length();
01303 }
01304 break;
01305 case Reflective:
01306 bct = Reflective;
01307 if ( face & 1 ) {
01308 coffset = 2*verts[d].max() + 1;
01309 voffset = 2*cells[d].max() + 1 - 1;
01310 } else {
01311 coffset = 2*verts[d].min() - 1;
01312 voffset = 2*cells[d].min() - 1 + 1;
01313 }
01314 break;
01315 case NoBC:
01316 bct = NoBC;
01317 if ( face & 1 ) {
01318 coffset = 2*verts[d].max() + 1;
01319 voffset = 2*cells[d].max() + 1 - 1;
01320 } else {
01321 coffset = 2*verts[d].min() - 1;
01322 voffset = 2*cells[d].min() - 1 + 1;
01323 }
01324 break;
01325 default:
01326 ERRORMSG("Cartesian::storeSpacingFields(): unknown MeshBC type" << endl);
01327 break;
01328 }
01329
01330
01331
01332 for (cfill_i=cellSpacings.begin_if();
01333 cfill_i!=cellSpacings.end_if(); ++cfill_i)
01334 {
01335
01336
01337 LField<T,Dim> &fill = *(*cfill_i).second;
01338
01339 const NDIndex<Dim> &fill_alloc = fill.getAllocated();
01340
01341
01342
01343
01344 if ( cSlab.touches( fill_alloc ) )
01345 {
01346
01347 NDIndex<Dim> dest = cSlab.intersect( fill_alloc );
01348
01349
01350
01351
01352
01353
01354 NDIndex<Dim> src = dest;
01355
01356
01357 src[d] = coffset - src[d];
01358
01359
01360
01361 typename BareField<T,Dim>::iterator_if from_i;
01362 for (from_i=cellSpacings.begin_if();
01363 from_i!=cellSpacings.end_if(); ++from_i)
01364 {
01365
01366 LField<T,Dim> &from = *(*from_i).second;
01367 const NDIndex<Dim> &from_owned = from.getOwned();
01368 const NDIndex<Dim> &from_alloc = from.getAllocated();
01369
01370 if ( src.touches( from_owned ) )
01371 {
01372 NDIndex<Dim> from_it = src.intersect( from_alloc );
01373 NDIndex<Dim> cfill_it = dest.plugBase( from_it );
01374
01375 typedef typename LField<T,Dim>::iterator LFI;
01376 LFI lhs = fill.begin(cfill_it);
01377 LFI rhs = from.begin(from_it);
01378
01379 if (bct == Periodic) {
01380 BrickExpression<Dim,LFI,LFI,OpMeshPeriodic<T> >
01381 (lhs,rhs,OpMeshPeriodic<T>()).apply();
01382 } else {
01383 if (bct == Reflective) {
01384 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01385 (lhs,rhs,OpMeshExtrapolate<T>(v0,v1)).apply();
01386 } else {
01387 if (bct == NoBC) {
01388 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01389 (lhs,rhs,OpMeshExtrapolate<T>(v0,v0)).apply();
01390 }
01391 }
01392 }
01393 }
01394 }
01395 }
01396 }
01397
01398 for (vfill_i=vertSpacings.begin_if();
01399 vfill_i!=vertSpacings.end_if(); ++vfill_i)
01400 {
01401
01402
01403 LField<T,Dim> &fill = *(*vfill_i).second;
01404
01405 const NDIndex<Dim> &fill_alloc = fill.getAllocated();
01406
01407
01408
01409
01410 if ( vSlab.touches( fill_alloc ) )
01411 {
01412
01413 NDIndex<Dim> dest = vSlab.intersect( fill_alloc );
01414
01415
01416
01417
01418
01419
01420 NDIndex<Dim> src = dest;
01421
01422
01423 src[d] = voffset - src[d];
01424
01425
01426
01427 typename BareField<T,Dim>::iterator_if from_i;
01428 for (from_i=vertSpacings.begin_if();
01429 from_i!=vertSpacings.end_if(); ++from_i)
01430 {
01431
01432 LField<T,Dim> &from = *(*from_i).second;
01433 const NDIndex<Dim> &from_owned = from.getOwned();
01434 const NDIndex<Dim> &from_alloc = from.getAllocated();
01435
01436 if ( src.touches( from_owned ) )
01437 {
01438 NDIndex<Dim> from_it = src.intersect( from_alloc );
01439 NDIndex<Dim> vfill_it = dest.plugBase( from_it );
01440
01441 typedef typename LField<T,Dim>::iterator LFI;
01442 LFI lhs = fill.begin(vfill_it);
01443 LFI rhs = from.begin(from_it);
01444
01445 if (bct == Periodic) {
01446 BrickExpression<Dim,LFI,LFI,OpMeshPeriodic<T> >
01447 (lhs,rhs,OpMeshPeriodic<T>()).apply();
01448 } else {
01449 if (bct == Reflective) {
01450 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01451 (lhs,rhs,OpMeshExtrapolate<T>(v0,v1)).apply();
01452 } else {
01453 if (bct == NoBC) {
01454 BrickExpression<Dim,LFI,LFI,OpMeshExtrapolate<T> >
01455 (lhs,rhs,OpMeshExtrapolate<T>(v0,v0)).apply();
01456 }
01457 }
01458 }
01459 }
01460 }
01461 }
01462 }
01463
01464 }
01465
01466 hasSpacingFields = true;
01467 }
01468
01469
01470
01471
01472
01473
01474 template< unsigned Dim, class MFLOAT >
01475 void
01476 Cartesian<Dim,MFLOAT>::
01477 print(ostream& out)
01478 {
01479 unsigned int d;
01480 out << "======Cartesian<" << Dim << ",MFLOAT>==begin======" << endl;
01481 for (d=0; d < Dim; d++)
01482 out << "gridSizes[" << d << "] = " << gridSizes[d] << endl;
01483 out << "origin = " << origin << endl;
01484 for (d=0; d < Dim; d++) {
01485 out << "--------meshSpacing[" << d << "]---------" << endl;
01486 typename map<int,MFLOAT>::iterator mi;
01487 for (mi=meshSpacing[d].begin(); mi != meshSpacing[d].end(); ++mi) {
01488 out << "meshSpacing[" << d << "][" << (*mi).first << "] = "
01489 << (*mi).second << endl;
01490 }
01491 }
01492 for (unsigned b=0; b < (1<<Dim); b++)
01493 out << "Dvc[" << b << "] = " << Dvc[b] << endl;
01494 for (d=0; d < Dim; d++)
01495 out << "MeshBC[" << 2*d << "] = " << Mesh<Dim>::MeshBC_E_Names[MeshBC[2*d]]
01496 << " ; MeshBC[" << 2*d+1 << "] = " << Mesh<Dim>::MeshBC_E_Names[MeshBC[2*d+1]]
01497 << endl;
01498 out << "======Cartesian<" << Dim << ",MFLOAT>==end========" << endl;
01499 }
01500
01501
01502
01503
01504
01505
01506 template <unsigned Dim, class MFLOAT>
01507 MFLOAT
01508 Cartesian<Dim,MFLOAT>::
01509 getCellVolume(const NDIndex<Dim>& ndi) const
01510 {
01511 MFLOAT volume = 1.0;
01512 int d;
01513 for (d=0; d<Dim; d++)
01514 if (ndi[d].length() != 1) {
01515 ERRORMSG("Cartesian::getCellVolume() error: arg is not a NDIndex"
01516 << "specifying a single element" << endl);
01517 }
01518 else {
01519 volume *= (*(meshSpacing[d].find(ndi[d].first()))).second;
01520 }
01521 return volume;
01522 }
01523
01524 template <unsigned Dim, class MFLOAT>
01525 Field<MFLOAT,Dim,Cartesian<Dim,MFLOAT>,Cell>&
01526 Cartesian<Dim,MFLOAT>::
01527 getCellVolumeField(Field<MFLOAT,Dim,Cartesian<Dim,MFLOAT>,Cell>& volumes) const
01528 {
01529
01530
01531
01532 volumes = 1.0;
01533 int d;
01534 int currentLocation[Dim];
01535 volumes.Uncompress();
01536
01537 typename Field<MFLOAT,Dim,Cartesian<Dim,MFLOAT>,Cell>::iterator fi,
01538 fi_end=volumes.end();
01539 for (fi = volumes.begin(); fi != fi_end; ++fi) {
01540 fi.GetCurrentLocation(currentLocation);
01541 for (d=0; d<Dim; d++)
01542 *fi *= (*(meshSpacing[d].find(currentLocation[d]))).second;
01543 }
01544 return volumes;
01545 }
01546
01547 template <unsigned Dim, class MFLOAT>
01548 MFLOAT
01549 Cartesian<Dim,MFLOAT>::
01550 getVertRangeVolume(const NDIndex<Dim>& ndi) const
01551 {
01552
01553 Vektor<MFLOAT,Dim> v0, v1;
01554 int d, i0, i1;
01555 for (d=0; d<Dim; d++) {
01556 i0 = ndi[d].first();
01557 if ( (i0 < -(int(gridSizes[d])-1)/2) ||
01558 (i0 > 3*(int(gridSizes[d])-1)/2) )
01559 ERRORMSG("Cartesian::getVertRangeVolume() error: " << ndi
01560 << " is an NDIndex ranging outside the mesh and guard layers;"
01561 << " not allowed." << endl);
01562 v0(d) = (*(meshPosition[d].find(i0))).second;
01563 i1 = ndi[d].last();
01564 if ( (i1 < -(int(gridSizes[d])-1)/2) ||
01565 (i1 > 3*(int(gridSizes[d])-1)/2) )
01566 ERRORMSG("Cartesian::getVertRangeVolume() error: " << ndi
01567 << " is an NDIndex ranging outside the mesh and guard layers;"
01568 << " not allowed." << endl);
01569 v1(d) = (*(meshPosition[d].find(i1))).second;
01570 }
01571
01572 MFLOAT volume = 1.0;
01573 for (d=0; d<Dim; d++) volume *= abs(v1(d) - v0(d));
01574 return volume;
01575 }
01576
01577 template <unsigned Dim, class MFLOAT>
01578 MFLOAT
01579 Cartesian<Dim,MFLOAT>::
01580 getCellRangeVolume(const NDIndex<Dim>& ndi) const
01581 {
01582
01583 Vektor<MFLOAT,Dim> v0, v1;
01584 int d, i0, i1;
01585 for (d=0; d<Dim; d++) {
01586 i0 = ndi[d].first();
01587 if ( (i0 < -(int(gridSizes[d])-1)/2) ||
01588 (i0 > 3*(int(gridSizes[d])-1)/2) )
01589 ERRORMSG("Cartesian::getCellRangeVolume() error: " << ndi
01590 << " is an NDIndex ranging outside the mesh and guard layers;"
01591 << " not allowed." << endl);
01592 v0(d) = (*(meshPosition[d].find(i0))).second;
01593 i1 = ndi[d].last()+1;
01594 if ( (i1 < -(int(gridSizes[d])-1)/2) ||
01595 (i1 > 3*(int(gridSizes[d])-1)/2) )
01596 ERRORMSG("Cartesian::getCellRangeVolume() error: " << ndi
01597 << " is an NDIndex ranging outside the mesh and guard layers;"
01598 << " not allowed." << endl);
01599 v1(d) = (*(meshPosition[d].find(i1))).second;
01600 }
01601
01602 MFLOAT volume = 1.0;
01603 for (d=0; d<Dim; d++) volume *= abs(v1(d) - v0(d));
01604 return volume;
01605 }
01606
01607
01608 template <unsigned Dim, class MFLOAT>
01609 NDIndex<Dim>
01610 Cartesian<Dim,MFLOAT>::
01611 getNearestVertex(const Vektor<MFLOAT,Dim>& x) const
01612 {
01613 int d;
01614 Vektor<MFLOAT,Dim> boxMin, boxMax;
01615 for (d=0; d<Dim; d++) {
01616 int gs = (int(gridSizes[d])-1)/2;
01617 boxMin(d) = (*(meshPosition[d].find(-gs))).second;
01618 boxMax(d) = (*(meshPosition[d].find(3*gs))).second;
01619 }
01620 for (d=0; d<Dim; d++)
01621 if ( (x(d) < boxMin(d)) || (x(d) > boxMax(d)) )
01622 ERRORMSG("Cartesian::getNearestVertex() - input point is outside"
01623 << " mesh boundary and guard layers; not allowed." << endl);
01624
01625
01626
01627 MFLOAT xVertexBelow, xVertexAbove, xVertex;
01628 int vertBelow, vertAbove, vertNearest[Dim];
01629 for (d=0; d<Dim; d++) {
01630 vertBelow = -(int(gridSizes[d])-1)/2;
01631 vertAbove = 3*(int(gridSizes[d])-1)/2;
01632 xVertexBelow = (*(meshPosition[d].find(vertBelow))).second;
01633 xVertexAbove = (*(meshPosition[d].find(vertAbove))).second;
01634
01635 if (x(d) < xVertexBelow) {
01636 vertNearest[d] = vertBelow;
01637 continue;
01638 }
01639 if (x(d) > xVertexAbove) {
01640 vertNearest[d] = vertAbove;
01641 continue;
01642 }
01643 while (vertAbove > vertBelow+1) {
01644 vertNearest[d] = (vertAbove+vertBelow)/2;
01645 xVertex = (*(meshPosition[d].find(vertNearest[d]))).second;
01646 if (x(d) > xVertex) {
01647 vertBelow = vertNearest[d];
01648 xVertexBelow = xVertex;
01649 }
01650 else if (x(d) < xVertex) {
01651 vertAbove = vertNearest[d];
01652 xVertexAbove = xVertex;
01653 }
01654 else {
01655 vertAbove = vertBelow;
01656 }
01657 }
01658 if (vertAbove != vertBelow) {
01659 if ((x(d)-xVertexBelow)<(xVertexAbove-x(d))) {
01660 vertNearest[d] = vertBelow;
01661 }
01662 else {
01663 vertNearest[d] = vertAbove;
01664 }
01665 }
01666 }
01667
01668
01669 NDIndex<Dim> ndi;
01670 for (d=0; d<Dim; d++) ndi[d] = Index(vertNearest[d],vertNearest[d],1);
01671
01672 return ndi;
01673 }
01674
01675 template <unsigned Dim, class MFLOAT>
01676 NDIndex<Dim>
01677 Cartesian<Dim,MFLOAT>::
01678 getVertexBelow(const Vektor<MFLOAT,Dim>& x) const
01679 {
01680 int d;
01681 Vektor<MFLOAT,Dim> boxMin, boxMax;
01682 for (d=0; d<Dim; d++) {
01683 int gs = (int(gridSizes[d]) - 1)/2;
01684 boxMin(d) = (*(meshPosition[d].find(-gs))).second;
01685 boxMax(d) = (*(meshPosition[d].find(3*gs))).second;
01686 }
01687 for (d=0; d<Dim; d++)
01688 if ( (x(d) < boxMin(d)) || (x(d) > boxMax(d)) )
01689 ERRORMSG("Cartesian::getVertexBelow() - input point is outside"
01690 << " mesh boundary and guard layers; not allowed." << endl);
01691
01692
01693 MFLOAT xVertexBelow, xVertexAbove, xVertex;
01694 int vertBelow, vertAbove, vertNearest[Dim];
01695 for (d=0; d<Dim; d++) {
01696 vertBelow = -(int(gridSizes[d])-1)/2;
01697 vertAbove = 3*(int(gridSizes[d])-1)/2;
01698 xVertexBelow = (*(meshPosition[d].find(vertBelow))).second;
01699 xVertexAbove = (*(meshPosition[d].find(vertAbove))).second;
01700
01701 if (x(d) < xVertexBelow) {
01702 vertNearest[d] = vertBelow;
01703 continue;
01704 }
01705 if (x(d) > xVertexAbove) {
01706 vertNearest[d] = vertAbove;
01707 continue;
01708 }
01709 while (vertAbove > vertBelow+1) {
01710 vertNearest[d] = (vertAbove+vertBelow)/2;
01711 xVertex = (*(meshPosition[d].find(vertNearest[d]))).second;
01712 if (x(d) > xVertex) {
01713 vertBelow = vertNearest[d];
01714 xVertexBelow = xVertex;
01715 }
01716 else if (x(d) < xVertex) {
01717 vertAbove = vertNearest[d];
01718 xVertexAbove = xVertex;
01719 }
01720 else {
01721 vertAbove = vertBelow;
01722 }
01723 }
01724 if (vertAbove != vertBelow) {
01725 vertNearest[d] = vertBelow;
01726 }
01727 }
01728
01729
01730 NDIndex<Dim> ndi;
01731 for (d=0; d<Dim; d++) ndi[d] = Index(vertNearest[d],vertNearest[d],1);
01732
01733 return ndi;
01734 }
01735
01736 template <unsigned Dim, class MFLOAT>
01737 Vektor<MFLOAT,Dim>
01738 Cartesian<Dim,MFLOAT>::
01739 getVertexPosition(const NDIndex<Dim>& ndi) const
01740 {
01741 int d, i;
01742 Vektor<MFLOAT,Dim> vertexPosition;
01743 for (d=0; d<Dim; d++) {
01744 if (ndi[d].length() != 1)
01745 ERRORMSG("Cartesian::getVertexPosition() error: " << ndi
01746 << " is not an NDIndex specifying a single element" << endl);
01747 i = ndi[d].first();
01748 if ( (i < -(int(gridSizes[d])-1)/2) ||
01749 (i > 3*(int(gridSizes[d])-1)/2) )
01750 ERRORMSG("Cartesian::getVertexPosition() error: " << ndi
01751 << " is an NDIndex outside the mesh and guard layers;"
01752 << " not allowed." << endl);
01753 vertexPosition(d) = (*(meshPosition[d].find(i))).second;
01754 }
01755 return vertexPosition;
01756 }
01757
01758 template <unsigned Dim, class MFLOAT>
01759 Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Vert>&
01760 Cartesian<Dim,MFLOAT>::
01761 getVertexPositionField(Field<Vektor<MFLOAT,Dim>,Dim,
01762 Cartesian<Dim,MFLOAT>,Vert>& vertexPositions) const
01763 {
01764 int d;
01765 int currentLocation[Dim];
01766 Vektor<MFLOAT,Dim> vertexPosition;
01767 vertexPositions.Uncompress();
01768 typename Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Vert>::iterator fi,
01769 fi_end = vertexPositions.end();
01770 for (fi = vertexPositions.begin(); fi != fi_end; ++fi) {
01771
01772 fi.GetCurrentLocation(currentLocation);
01773 for (d=0; d<Dim; d++) {
01774 vertexPosition(d) = (*(meshPosition[d].find(currentLocation[d]))).second;
01775 }
01776 *fi = vertexPosition;
01777 }
01778 return vertexPositions;
01779 }
01780
01781
01782 template <unsigned Dim, class MFLOAT>
01783 Vektor<MFLOAT,Dim>
01784 Cartesian<Dim,MFLOAT>::
01785 getCellPosition(const NDIndex<Dim>& ndi) const
01786 {
01787 int d, i;
01788 Vektor<MFLOAT,Dim> cellPosition;
01789 for (d=0; d<Dim; d++) {
01790 if (ndi[d].length() != 1)
01791 ERRORMSG("Cartesian::getCellPosition() error: " << ndi
01792 << " is not an NDIndex specifying a single element" << endl);
01793 i = ndi[d].first();
01794 if ( (i < -(int(gridSizes[d])-1)/2) ||
01795 (i >= 3*(int(gridSizes[d])-1)/2) )
01796 ERRORMSG("Cartesian::getCellPosition() error: " << ndi
01797 << " is an NDIndex outside the mesh and guard layers;"
01798 << " not allowed." << endl);
01799 cellPosition(d) = 0.5 * ( (*(meshPosition[d].find(i))).second +
01800 (*(meshPosition[d].find(i+1))).second );
01801 }
01802 return cellPosition;
01803 }
01804
01805 template <unsigned Dim, class MFLOAT>
01806 Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Cell>&
01807 Cartesian<Dim,MFLOAT>::
01808 getCellPositionField(Field<Vektor<MFLOAT,Dim>,Dim,
01809 Cartesian<Dim,MFLOAT>,Cell>& cellPositions) const
01810 {
01811 int d;
01812 int currentLocation[Dim];
01813 Vektor<MFLOAT,Dim> cellPosition;
01814 cellPositions.Uncompress();
01815 typename Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Cell>::iterator fi,
01816 fi_end = cellPositions.end();
01817 for (fi = cellPositions.begin(); fi != fi_end; ++fi) {
01818
01819 fi.GetCurrentLocation(currentLocation);
01820 for (d=0; d<Dim; d++) {
01821 cellPosition(d) =
01822 0.5 * ( (*(meshPosition[d].find(currentLocation[d]))).second +
01823 (*(meshPosition[d].find(currentLocation[d]+1))).second );
01824 }
01825 *fi = cellPosition;
01826 }
01827 return cellPositions;
01828 }
01829
01830
01831 template <unsigned Dim, class MFLOAT>
01832 Vektor<MFLOAT,Dim>
01833 Cartesian<Dim,MFLOAT>::
01834 getDeltaVertex(const NDIndex<Dim>& ndi) const
01835 {
01836
01837 Vektor<MFLOAT,Dim> vertexVertexSpacing(0);
01838
01839 for (int d=0; d<Dim; d++) {
01840
01841 int a = ndi[d].first();
01842 int b = ndi[d].last();
01843 if (b < a) {
01844 int tmpa = a; a = b; b = tmpa;
01845 }
01846
01847
01848 if (a < -((int(gridSizes[d])-1)/2) || b >= 3*(int(gridSizes[d])-1)/2) {
01849 ERRORMSG("Cartesian::getDeltaVertex() error: " << ndi
01850 << " is an NDIndex ranging outside"
01851 << " the mesh and guard layers region; not allowed."
01852 << endl);
01853 }
01854
01855
01856
01857 while (a <= b)
01858 vertexVertexSpacing[d] += (*(meshSpacing[d].find(a++))).second;
01859 }
01860
01861 return vertexVertexSpacing;
01862 }
01863
01864
01865 template <unsigned Dim, class MFLOAT>
01866 Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Cell>&
01867 Cartesian<Dim,MFLOAT>::
01868 getDeltaVertexField(Field<Vektor<MFLOAT,Dim>,Dim,
01869 Cartesian<Dim,MFLOAT>,Cell>& vertexSpacings) const
01870 {
01871 int currentLocation[Dim];
01872 Vektor<MFLOAT,Dim> vertexVertexSpacing;
01873 vertexSpacings.Uncompress();
01874 typename Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Cell>::iterator fi,
01875 fi_end = vertexSpacings.end();
01876 for (fi = vertexSpacings.begin(); fi != fi_end; ++fi) {
01877 fi.GetCurrentLocation(currentLocation);
01878 for (unsigned int d=0; d<Dim; d++)
01879 vertexVertexSpacing[d]=(*(meshSpacing[d].find(currentLocation[d]))).second;
01880 *fi = vertexVertexSpacing;
01881 }
01882 return vertexSpacings;
01883 }
01884
01885
01886 template <unsigned Dim, class MFLOAT>
01887 Vektor<MFLOAT,Dim>
01888 Cartesian<Dim,MFLOAT>::
01889 getDeltaCell(const NDIndex<Dim>& ndi) const
01890 {
01891
01892 Vektor<MFLOAT,Dim> cellCellSpacing(0);
01893
01894 for (int d=0; d<Dim; d++) {
01895
01896 int a = ndi[d].first();
01897 int b = ndi[d].last();
01898 if (b < a) {
01899 int tmpa = a; a = b; b = tmpa;
01900 }
01901
01902
01903 if (a <= -(int(gridSizes[d])-1)/2 || b >= 3*(int(gridSizes[d])-1)/2) {
01904 ERRORMSG("Cartesian::getDeltaCell() error: " << ndi
01905 << " is an NDIndex ranging outside"
01906 << " the mesh and guard layers region; not allowed."
01907 << endl);
01908 }
01909
01910
01911 while (a <= b) {
01912 cellCellSpacing[d] += ((*(meshSpacing[d].find(a))).second +
01913 (*(meshSpacing[d].find(a-1))).second) * 0.5;
01914 a++;
01915 }
01916 }
01917 return cellCellSpacing;
01918 }
01919
01920
01921 template <unsigned Dim, class MFLOAT>
01922 Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Vert>&
01923 Cartesian<Dim,MFLOAT>::
01924 getDeltaCellField(Field<Vektor<MFLOAT,Dim>,Dim,
01925 Cartesian<Dim,MFLOAT>,Vert>& cellSpacings) const
01926 {
01927 int currentLocation[Dim];
01928 Vektor<MFLOAT,Dim> cellCellSpacing;
01929 cellSpacings.Uncompress();
01930 typename Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Vert>::iterator fi,
01931 fi_end = cellSpacings.end();
01932 for (fi = cellSpacings.begin(); fi != fi_end; ++fi) {
01933 fi.GetCurrentLocation(currentLocation);
01934 for (unsigned int d=0; d<Dim; d++)
01935 cellCellSpacing[d]+=((*(meshSpacing[d].find(currentLocation[d]))).second +
01936 (*(meshSpacing[d].find(currentLocation[d]-1))).second) * 0.5;
01937 *fi = cellCellSpacing;
01938 }
01939 return cellSpacings;
01940 }
01941
01942 template <unsigned Dim, class MFLOAT>
01943 Vektor<MFLOAT,Dim>*
01944 Cartesian<Dim,MFLOAT>::
01945 getSurfaceNormals(const NDIndex<Dim>& ndi) const
01946 {
01947 Vektor<MFLOAT,Dim>* surfaceNormals = new Vektor<MFLOAT,Dim>[2*Dim];
01948 int d, i;
01949 for (d=0; d<Dim; d++) {
01950 for (i=0; i<Dim; i++) {
01951 surfaceNormals[2*d](i) = 0.0;
01952 surfaceNormals[2*d+1](i) = 0.0;
01953 }
01954 surfaceNormals[2*d](d) = -1.0;
01955 surfaceNormals[2*d+1](d) = 1.0;
01956 }
01957 return surfaceNormals;
01958 }
01959
01960 template <unsigned Dim, class MFLOAT>
01961 void
01962 Cartesian<Dim,MFLOAT>::
01963 getSurfaceNormalFields(Field<Vektor<MFLOAT,Dim>, Dim,
01964 Cartesian<Dim,MFLOAT>,Cell>**
01965 surfaceNormalsFields ) const
01966 {
01967 Vektor<MFLOAT,Dim>* surfaceNormals = new Vektor<MFLOAT,Dim>[2*Dim];
01968 int d, i;
01969 for (d=0; d<Dim; d++) {
01970 for (i=0; i<Dim; i++) {
01971 surfaceNormals[2*d](i) = 0.0;
01972 surfaceNormals[2*d+1](i) = 0.0;
01973 }
01974 surfaceNormals[2*d](d) = -1.0;
01975 surfaceNormals[2*d+1](d) = 1.0;
01976 }
01977 for (d=0; d<2*Dim; d++) assign((*(surfaceNormalsFields[d])),
01978 surfaceNormals[d]);
01979
01980 }
01981
01982
01983
01984
01985
01986 template <unsigned Dim, class MFLOAT>
01987 Vektor<MFLOAT,Dim>
01988 Cartesian<Dim,MFLOAT>::
01989 getSurfaceNormal(const NDIndex<Dim>& ndi, unsigned face) const
01990 {
01991 Vektor<MFLOAT,Dim> surfaceNormal;
01992 unsigned int d;
01993
01994
01995
01996
01997 if ( face & 1 ) {
01998 for (d=0; d<Dim; d++) {
01999 if ((face/2) == d) {
02000 surfaceNormal(face) = -1.0;
02001 } else {
02002 surfaceNormal(face) = 0.0;
02003 }
02004 }
02005 } else {
02006 for (d=0; d<Dim; d++) {
02007 if ((face/2) == d) {
02008 surfaceNormal(face) = 1.0;
02009 } else {
02010 surfaceNormal(face) = 0.0;
02011 }
02012 }
02013 }
02014 return surfaceNormal;
02015 }
02016
02017 template <unsigned Dim, class MFLOAT>
02018 Field<Vektor<MFLOAT,Dim>,Dim,Cartesian<Dim,MFLOAT>,Cell>&
02019 Cartesian<Dim,MFLOAT>::
02020 getSurfaceNormalField(Field<Vektor<MFLOAT,Dim>, Dim,
02021 Cartesian<Dim,MFLOAT>,Cell>& surfaceNormalField,
02022 unsigned face) const
02023 {
02024 Vektor<MFLOAT,Dim> surfaceNormal;
02025 unsigned int d;
02026
02027
02028
02029
02030 if ( face & 1 ) {
02031 for (d=0; d<Dim; d++) {
02032 if ((face/2) == d) {
02033 surfaceNormal(face) = -1.0;
02034 } else {
02035 surfaceNormal(face) = 0.0;
02036 }
02037 }
02038 } else {
02039 for (d=0; d<Dim; d++) {
02040 if ((face/2) == d) {
02041 surfaceNormal(face) = 1.0;
02042 } else {
02043 surfaceNormal(face) = 0.0;
02044 }
02045 }
02046 }
02047 surfaceNormalField = surfaceNormal;
02048 return surfaceNormalField;
02049 }
02050
02051
02052
02053
02054
02055 template <unsigned Dim, class MFLOAT>
02056 void
02057 Cartesian<Dim,MFLOAT>::
02058 set_MeshBC(unsigned face, MeshBC_E meshBCType)
02059 {
02060 MeshBC[face] = meshBCType;
02061 updateMeshSpacingGuards(face);
02062
02063 if (hasSpacingFields) storeSpacingFields();
02064 }
02065
02066 template <unsigned Dim, class MFLOAT>
02067 void
02068 Cartesian<Dim,MFLOAT>::
02069 set_MeshBC(MeshBC_E* meshBCTypes)
02070 {
02071 for (int face=0; face < 2*Dim; face++) {
02072 MeshBC[face] = meshBCTypes[face];
02073 updateMeshSpacingGuards(face);
02074 }
02075
02076 if (hasSpacingFields) storeSpacingFields();
02077 }
02078
02079 template <unsigned Dim, class MFLOAT>
02080 void
02081 Cartesian<Dim,MFLOAT>::
02082 updateMeshSpacingGuards(int face)
02083 {
02084
02085
02086
02087 int d = face/2;
02088 int cell, guardLayer;
02089
02090
02091
02092
02093 if ( face & 1 ) {
02094
02095 switch (MeshBC[d*2]) {
02096 case Periodic:
02097 for (guardLayer = 0; guardLayer <= (gridSizes[d]-1)/2; guardLayer++) {
02098 cell = gridSizes[d] - 1 + guardLayer;
02099 (meshSpacing[d])[cell] = (meshSpacing[d])[guardLayer];
02100 (meshPosition[d])[cell+1] = (meshPosition[d])[cell] +
02101 (meshSpacing[d])[cell];
02102 }
02103 break;
02104 case Reflective:
02105 for (guardLayer = 0; guardLayer <= (gridSizes[d]-1)/2; guardLayer++) {
02106 cell = gridSizes[d] - 1 + guardLayer;
02107 (meshSpacing[d])[cell] = (meshSpacing[d])[cell - guardLayer - 1];
02108 (meshPosition[d])[cell+1] = (meshPosition[d])[cell] +
02109 (meshSpacing[d])[cell];
02110 }
02111 break;
02112 case NoBC:
02113 for (guardLayer = 0; guardLayer <= (gridSizes[d]-1)/2; guardLayer++) {
02114 cell = gridSizes[d] - 1 + guardLayer;
02115 (meshSpacing[d])[cell] = 0;
02116 (meshPosition[d])[cell+1] = (meshPosition[d])[cell] +
02117 (meshSpacing[d])[cell];
02118 }
02119 break;
02120 default:
02121 ERRORMSG("Cartesian::updateMeshSpacingGuards(): unknown MeshBC type"
02122 << endl);
02123 break;
02124 }
02125 }
02126 else {
02127
02128 switch (MeshBC[d]) {
02129 case Periodic:
02130 for (guardLayer = 0; guardLayer <= (gridSizes[d]-1)/2; guardLayer++) {
02131 cell = -1 - guardLayer;
02132 (meshSpacing[d])[cell] = (meshSpacing[d])[gridSizes[d] + cell];
02133 (meshPosition[d])[cell] = (meshPosition[d])[cell+1] -
02134 (meshSpacing[d])[cell];
02135 }
02136 break;
02137 case Reflective:
02138 for (guardLayer = 0; guardLayer <= (gridSizes[d]-1)/2; guardLayer++) {
02139 cell = -1 - guardLayer;
02140 (meshSpacing[d])[cell] = (meshSpacing[d])[-cell - 1];
02141 (meshPosition[d])[cell] = (meshPosition[d])[cell+1] -
02142 (meshSpacing[d])[cell];
02143 }
02144 break;
02145 case NoBC:
02146 for (guardLayer = 0; guardLayer <= (gridSizes[d]-1)/2; guardLayer++) {
02147 cell = -1 - guardLayer;
02148 (meshSpacing[d])[cell] = 0;
02149 (meshPosition[d])[cell] = (meshPosition[d])[cell+1] -
02150 (meshSpacing[d])[cell];
02151 }
02152 break;
02153 default:
02154 ERRORMSG("Cartesian::updateMeshSpacingGuards(): unknown MeshBC type"
02155 << endl);
02156 break;
02157 }
02158 }
02159 }
02160
02161
02162
02163 template <unsigned Dim, class MFLOAT>
02164 MeshBC_E
02165 Cartesian<Dim,MFLOAT>::
02166 get_MeshBC(unsigned face) const
02167 {
02168 MeshBC_E mb;
02169 mb = MeshBC[face];
02170 return mb;
02171 }
02172
02173 template <unsigned Dim, class MFLOAT>
02174 MeshBC_E*
02175 Cartesian<Dim,MFLOAT>::
02176 get_MeshBC() const
02177 {
02178 MeshBC_E* mb = new MeshBC_E[2*Dim];
02179 for (int b=0; b < 2*Dim; b++) mb[b] = MeshBC[b];
02180 return mb;
02181 }
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197 template < class T, class MFLOAT >
02198 Field<T,1U,Cartesian<1U,MFLOAT>,Cell>&
02199 Div(Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02200 Field<T,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02201 {
02202 PAssert(x.get_mesh().hasSpacingFields);
02203 BareField<Vektor<MFLOAT,1U>,1U>& vertSpacings = *(x.get_mesh().VertSpacings);
02204 const NDIndex<1U>& domain = r.getDomain();
02205 Index I = domain[0];
02206 r[I] =
02207 dot(x[I ], x.get_mesh().Dvc[0]/vertSpacings[I]) +
02208 dot(x[I+1], x.get_mesh().Dvc[1]/vertSpacings[I]);
02209 return r;
02210 }
02211
02212 template < class T, class MFLOAT >
02213 Field<T,2U,Cartesian<2U,MFLOAT>,Cell>&
02214 Div(Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& x,
02215 Field<T,2U,Cartesian<2U,MFLOAT>,Cell>& r)
02216 {
02217 PAssert(x.get_mesh().hasSpacingFields);
02218 BareField<Vektor<MFLOAT,2U>,2U>& vertSpacings = *(x.get_mesh().VertSpacings);
02219 const NDIndex<2U>& domain = r.getDomain();
02220 Index I = domain[0];
02221 Index J = domain[1];
02222 r[I][J] =
02223 dot(x[I ][J ], x.get_mesh().Dvc[0]/vertSpacings[I][J]) +
02224 dot(x[I+1][J ], x.get_mesh().Dvc[1]/vertSpacings[I][J]) +
02225 dot(x[I ][J+1], x.get_mesh().Dvc[2]/vertSpacings[I][J]) +
02226 dot(x[I+1][J+1], x.get_mesh().Dvc[3]/vertSpacings[I][J]);
02227 return r;
02228 }
02229
02230 template < class T, class MFLOAT >
02231 Field<T,3U,Cartesian<3U,MFLOAT>,Cell>&
02232 Div(Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& x,
02233 Field<T,3U,Cartesian<3U,MFLOAT>,Cell>& r)
02234 {
02235 PAssert(x.get_mesh().hasSpacingFields);
02236 BareField<Vektor<MFLOAT,3U>,3U>& vertSpacings = *(x.get_mesh().VertSpacings);
02237 const NDIndex<3U>& domain = r.getDomain();
02238 Index I = domain[0];
02239 Index J = domain[1];
02240 Index K = domain[2];
02241 r[I][J][K] =
02242 dot(x[I ][J ][K ], x.get_mesh().Dvc[0]/vertSpacings[I][J][K]) +
02243 dot(x[I+1][J ][K ], x.get_mesh().Dvc[1]/vertSpacings[I][J][K]) +
02244 dot(x[I ][J+1][K ], x.get_mesh().Dvc[2]/vertSpacings[I][J][K]) +
02245 dot(x[I+1][J+1][K ], x.get_mesh().Dvc[3]/vertSpacings[I][J][K]) +
02246 dot(x[I ][J ][K+1], x.get_mesh().Dvc[4]/vertSpacings[I][J][K]) +
02247 dot(x[I+1][J ][K+1], x.get_mesh().Dvc[5]/vertSpacings[I][J][K]) +
02248 dot(x[I ][J+1][K+1], x.get_mesh().Dvc[6]/vertSpacings[I][J][K]) +
02249 dot(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]/vertSpacings[I][J][K]);
02250 return r;
02251 }
02252
02253
02254
02255 template < class T, class MFLOAT >
02256 Field<T,1U,Cartesian<1U,MFLOAT>,Vert>&
02257 Div(Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& x,
02258 Field<T,1U,Cartesian<1U,MFLOAT>,Vert>& r)
02259 {
02260 PAssert(x.get_mesh().hasSpacingFields);
02261 BareField<Vektor<MFLOAT,1U>,1U>& cellSpacings = *(x.get_mesh().CellSpacings);
02262 const NDIndex<1U>& domain = r.getDomain();
02263 Index I = domain[0];
02264 r[I] =
02265 dot(x[I-1], x.get_mesh().Dvc[0]/cellSpacings[I]) +
02266 dot(x[I ], x.get_mesh().Dvc[1]/cellSpacings[I]);
02267 return r;
02268 }
02269
02270 template < class T, class MFLOAT >
02271 Field<T,2U,Cartesian<2U,MFLOAT>,Vert>&
02272 Div(Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& x,
02273 Field<T,2U,Cartesian<2U,MFLOAT>,Vert>& r)
02274 {
02275 PAssert(x.get_mesh().hasSpacingFields);
02276 BareField<Vektor<MFLOAT,2U>,2U>& cellSpacings = *(x.get_mesh().CellSpacings);
02277 const NDIndex<2U>& domain = r.getDomain();
02278 Index I = domain[0];
02279 Index J = domain[1];
02280 r[I][J] =
02281 dot(x[I-1][J-1], x.get_mesh().Dvc[0]/cellSpacings[I][J]) +
02282 dot(x[I ][J-1], x.get_mesh().Dvc[1]/cellSpacings[I][J]) +
02283 dot(x[I-1][J ], x.get_mesh().Dvc[2]/cellSpacings[I][J]) +
02284 dot(x[I ][J ], x.get_mesh().Dvc[3]/cellSpacings[I][J]);
02285 return r;
02286 }
02287
02288 template < class T, class MFLOAT >
02289 Field<T,3U,Cartesian<3U,MFLOAT>,Vert>&
02290 Div(Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& x,
02291 Field<T,3U,Cartesian<3U,MFLOAT>,Vert>& r)
02292 {
02293 PAssert(x.get_mesh().hasSpacingFields);
02294 BareField<Vektor<MFLOAT,3U>,3U>& cellSpacings = *(x.get_mesh().CellSpacings);
02295 const NDIndex<3U>& domain = r.getDomain();
02296 Index I = domain[0];
02297 Index J = domain[1];
02298 Index K = domain[2];
02299 r[I][J][K] =
02300 dot(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]/cellSpacings[I][J][K]) +
02301 dot(x[I ][J-1][K-1], x.get_mesh().Dvc[1]/cellSpacings[I][J][K]) +
02302 dot(x[I-1][J ][K-1], x.get_mesh().Dvc[2]/cellSpacings[I][J][K]) +
02303 dot(x[I ][J ][K-1], x.get_mesh().Dvc[3]/cellSpacings[I][J][K]) +
02304 dot(x[I-1][J-1][K ], x.get_mesh().Dvc[4]/cellSpacings[I][J][K]) +
02305 dot(x[I ][J-1][K ], x.get_mesh().Dvc[5]/cellSpacings[I][J][K]) +
02306 dot(x[I-1][J ][K ], x.get_mesh().Dvc[6]/cellSpacings[I][J][K]) +
02307 dot(x[I ][J ][K ], x.get_mesh().Dvc[7]/cellSpacings[I][J][K]);
02308 return r;
02309 }
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328 template < class T, class MFLOAT >
02329 Field<T,1U,Cartesian<1U,MFLOAT>,Vert>&
02330 Div(Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02331 Field<T,1U,Cartesian<1U,MFLOAT>,Vert>& r)
02332 {
02333 PAssert(x.get_mesh().hasSpacingFields);
02334 BareField<Vektor<MFLOAT,1U>,1U>& vS = *(x.get_mesh().VertSpacings);
02335 const NDIndex<1U>& domain = r.getDomain();
02336 Index I = domain[0];
02337 Vektor<MFLOAT,1U> idx;
02338 idx[0] = 1.0;
02339 r[I] = dot(idx, (x[I+1] - x[I-1])/(vS[I ] + vS[I-1]));
02340 return r;
02341 }
02342
02343 template < class T, class MFLOAT >
02344 Field<T,2U,Cartesian<2U,MFLOAT>,Vert>&
02345 Div(Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& x,
02346 Field<T,2U,Cartesian<2U,MFLOAT>,Vert>& r)
02347 {
02348 PAssert(x.get_mesh().hasSpacingFields);
02349 BareField<Vektor<MFLOAT,2U>,2U>& vS = *(x.get_mesh().VertSpacings);
02350 const NDIndex<2U>& domain = r.getDomain();
02351 Index I = domain[0];
02352 Index J = domain[1];
02353 Vektor<MFLOAT,2U> idx,idy;
02354 idx[0] = 1.0;
02355 idx[1] = 0.0;
02356 idy[0] = 0.0;
02357 idy[1] = 1.0;
02358 r[I][J] =
02359 dot(idx, (x[I+1][J ] - x[I-1][J ])/(vS[I ][J ] + vS[I-1][J ])) +
02360 dot(idy, (x[I ][J+1] - x[I ][J-1])/(vS[I ][J ] + vS[I ][J-1]));
02361 return r;
02362 }
02363
02364 template < class T, class MFLOAT >
02365 Field<T,3U,Cartesian<3U,MFLOAT>,Vert>&
02366 Div(Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& x,
02367 Field<T,3U,Cartesian<3U,MFLOAT>,Vert>& r)
02368 {
02369 PAssert(x.get_mesh().hasSpacingFields);
02370 BareField<Vektor<MFLOAT,3U>,3U>& vS = *(x.get_mesh().VertSpacings);
02371 const NDIndex<3U>& domain = r.getDomain();
02372 Index I = domain[0];
02373 Index J = domain[1];
02374 Index K = domain[2];
02375 Vektor<MFLOAT,3U> idx,idy,idz;
02376 idx[0] = 1.0;
02377 idx[1] = 0.0;
02378 idx[2] = 0.0;
02379 idy[0] = 0.0;
02380 idy[1] = 1.0;
02381 idy[2] = 0.0;
02382 idz[0] = 0.0;
02383 idz[1] = 0.0;
02384 idz[2] = 1.0;
02385 r[I][J][K] =
02386 dot(idx, ((x[I+1][J ][K ] - x[I-1][J ][K ])/
02387 (vS[I ][J ][K ] + vS[I-1][J ][K ]))) +
02388 dot(idy, ((x[I ][J+1][K ] - x[I ][J-1][K ])/
02389 (vS[I ][J ][K ] + vS[I ][J-1][K ]))) +
02390 dot(idz, ((x[I ][J ][K+1] - x[I ][J ][K-1])/
02391 (vS[I ][J ][K ] + vS[I ][J ][K-1])));
02392 return r;
02393 }
02394
02395
02396
02397
02398
02399
02400
02401 template < class T, class MFLOAT >
02402 Field<T,1U,Cartesian<1U,MFLOAT>,Cell>&
02403 Div(Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& x,
02404 Field<T,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02405 {
02406 PAssert(x.get_mesh().hasSpacingFields);
02407 BareField<Vektor<MFLOAT,1U>,1U>& cS = *(x.get_mesh().CellSpacings);
02408 const NDIndex<1U>& domain = r.getDomain();
02409 Index I = domain[0];
02410 Vektor<MFLOAT,1U> idx;
02411 idx[0] = 1.0;
02412 r[I] = dot(idx, (x[I+1] - x[I-1])/(cS[I ] + cS[I-1]));
02413 return r;
02414 }
02415
02416 template < class T, class MFLOAT >
02417 Field<T,2U,Cartesian<2U,MFLOAT>,Cell>&
02418 Div(Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& x,
02419 Field<T,2U,Cartesian<2U,MFLOAT>,Cell>& r)
02420 {
02421 PAssert(x.get_mesh().hasSpacingFields);
02422 BareField<Vektor<MFLOAT,2U>,2U>& cS = *(x.get_mesh().CellSpacings);
02423 const NDIndex<2U>& domain = r.getDomain();
02424 Index I = domain[0];
02425 Index J = domain[1];
02426 Vektor<MFLOAT,2U> idx,idy;
02427 idx[0] = 1.0;
02428 idx[1] = 0.0;
02429 idy[0] = 0.0;
02430 idy[1] = 1.0;
02431 r[I][J] =
02432 dot(idx, (x[I+1][J ] - x[I-1][J ])/(cS[I ][J ] + cS[I-1][J ])) +
02433 dot(idy, (x[I ][J+1] - x[I ][J-1])/(cS[I ][J ] + cS[I ][J-1]));
02434 return r;
02435 }
02436
02437 template < class T, class MFLOAT >
02438 Field<T,3U,Cartesian<3U,MFLOAT>,Cell>&
02439 Div(Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& x,
02440 Field<T,3U,Cartesian<3U,MFLOAT>,Cell>& r)
02441 {
02442 PAssert(x.get_mesh().hasSpacingFields);
02443 BareField<Vektor<MFLOAT,3U>,3U>& cS = *(x.get_mesh().CellSpacings);
02444 const NDIndex<3U>& domain = r.getDomain();
02445 Index I = domain[0];
02446 Index J = domain[1];
02447 Index K = domain[2];
02448 Vektor<MFLOAT,3U> idx,idy,idz;
02449 idx[0] = 1.0;
02450 idx[1] = 0.0;
02451 idx[2] = 0.0;
02452 idy[0] = 0.0;
02453 idy[1] = 1.0;
02454 idy[2] = 0.0;
02455 idz[0] = 0.0;
02456 idz[1] = 0.0;
02457 idz[2] = 1.0;
02458 r[I][J][K] =
02459 dot(idx, ((x[I+1][J ][K ] - x[I-1][J ][K ])/
02460 (cS[I ][J ][K ] + cS[I-1][J ][K ]))) +
02461 dot(idy, ((x[I ][J+1][K ] - x[I ][J-1][K ])/
02462 (cS[I ][J ][K ] + cS[I ][J-1][K ]))) +
02463 dot(idz, ((x[I ][J ][K+1] - x[I ][J ][K-1])/
02464 (cS[I ][J ][K ] + cS[I ][J ][K-1])));
02465 return r;
02466 }
02467
02468
02469
02470 template < class T, class MFLOAT >
02471 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>&
02472 Div(Field<Tenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02473 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02474 {
02475 PAssert(x.get_mesh().hasSpacingFields);
02476 BareField<Vektor<MFLOAT,1U>,1U>& vertSpacings = *(x.get_mesh().VertSpacings);
02477 const NDIndex<1U>& domain = r.getDomain();
02478 Index I = domain[0];
02479 r[I] =
02480 dot(x[I ], x.get_mesh().Dvc[0]/vertSpacings[I]) +
02481 dot(x[I+1], x.get_mesh().Dvc[1]/vertSpacings[I]);
02482 return r;
02483 }
02484
02485 template < class T, class MFLOAT >
02486 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>&
02487 Div(Field<Tenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& x,
02488 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& r)
02489 {
02490 PAssert(x.get_mesh().hasSpacingFields);
02491 BareField<Vektor<MFLOAT,2U>,2U>& vertSpacings = *(x.get_mesh().VertSpacings);
02492 const NDIndex<2U>& domain = r.getDomain();
02493 Index I = domain[0];
02494 Index J = domain[1];
02495 r[I][J] =
02496 dot(x[I ][J ], x.get_mesh().Dvc[0]/vertSpacings[I][J]) +
02497 dot(x[I+1][J ], x.get_mesh().Dvc[1]/vertSpacings[I][J]) +
02498 dot(x[I ][J+1], x.get_mesh().Dvc[2]/vertSpacings[I][J]) +
02499 dot(x[I+1][J+1], x.get_mesh().Dvc[3]/vertSpacings[I][J]);
02500 return r;
02501 }
02502
02503 template < class T, class MFLOAT >
02504 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>&
02505 Div(Field<Tenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& x,
02506 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& r)
02507 {
02508 PAssert(x.get_mesh().hasSpacingFields);
02509 BareField<Vektor<MFLOAT,3U>,3U>& vertSpacings = *(x.get_mesh().VertSpacings);
02510 const NDIndex<3U>& domain = r.getDomain();
02511 Index I = domain[0];
02512 Index J = domain[1];
02513 Index K = domain[2];
02514 r[I][J][K] =
02515 dot(x[I ][J ][K ], x.get_mesh().Dvc[0]/vertSpacings[I][J][K]) +
02516 dot(x[I+1][J ][K ], x.get_mesh().Dvc[1]/vertSpacings[I][J][K]) +
02517 dot(x[I ][J+1][K ], x.get_mesh().Dvc[2]/vertSpacings[I][J][K]) +
02518 dot(x[I+1][J+1][K ], x.get_mesh().Dvc[3]/vertSpacings[I][J][K]) +
02519 dot(x[I ][J ][K+1], x.get_mesh().Dvc[4]/vertSpacings[I][J][K]) +
02520 dot(x[I+1][J ][K+1], x.get_mesh().Dvc[5]/vertSpacings[I][J][K]) +
02521 dot(x[I ][J+1][K+1], x.get_mesh().Dvc[6]/vertSpacings[I][J][K]) +
02522 dot(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]/vertSpacings[I][J][K]);
02523 return r;
02524 }
02525
02526
02527
02528 template < class T, class MFLOAT >
02529 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>&
02530 Div(Field<SymTenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02531 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02532 {
02533 PAssert(x.get_mesh().hasSpacingFields);
02534 BareField<Vektor<MFLOAT,1U>,1U>& vertSpacings = *(x.get_mesh().VertSpacings);
02535 const NDIndex<1U>& domain = r.getDomain();
02536 Index I = domain[0];
02537 r[I] =
02538 dot(x[I ], x.get_mesh().Dvc[0]/vertSpacings[I]) +
02539 dot(x[I+1], x.get_mesh().Dvc[1]/vertSpacings[I]);
02540 return r;
02541 }
02542
02543 template < class T, class MFLOAT >
02544 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>&
02545 Div(Field<SymTenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& x,
02546 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& r)
02547 {
02548 PAssert(x.get_mesh().hasSpacingFields);
02549 BareField<Vektor<MFLOAT,2U>,2U>& vertSpacings = *(x.get_mesh().VertSpacings);
02550 const NDIndex<2U>& domain = r.getDomain();
02551 Index I = domain[0];
02552 Index J = domain[1];
02553 r[I][J] =
02554 dot(x[I ][J ], x.get_mesh().Dvc[0]/vertSpacings[I][J]) +
02555 dot(x[I+1][J ], x.get_mesh().Dvc[1]/vertSpacings[I][J]) +
02556 dot(x[I ][J+1], x.get_mesh().Dvc[2]/vertSpacings[I][J]) +
02557 dot(x[I+1][J+1], x.get_mesh().Dvc[3]/vertSpacings[I][J]);
02558 return r;
02559 }
02560
02561 template < class T, class MFLOAT >
02562 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>&
02563 Div(Field<SymTenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& x,
02564 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& r)
02565 {
02566 PAssert(x.get_mesh().hasSpacingFields);
02567 BareField<Vektor<MFLOAT,3U>,3U>& vertSpacings = *(x.get_mesh().VertSpacings);
02568 const NDIndex<3U>& domain = r.getDomain();
02569 Index I = domain[0];
02570 Index J = domain[1];
02571 Index K = domain[2];
02572 r[I][J][K] =
02573 dot(x[I ][J ][K ], x.get_mesh().Dvc[0]/vertSpacings[I][J][K]) +
02574 dot(x[I+1][J ][K ], x.get_mesh().Dvc[1]/vertSpacings[I][J][K]) +
02575 dot(x[I ][J+1][K ], x.get_mesh().Dvc[2]/vertSpacings[I][J][K]) +
02576 dot(x[I+1][J+1][K ], x.get_mesh().Dvc[3]/vertSpacings[I][J][K]) +
02577 dot(x[I ][J ][K+1], x.get_mesh().Dvc[4]/vertSpacings[I][J][K]) +
02578 dot(x[I+1][J ][K+1], x.get_mesh().Dvc[5]/vertSpacings[I][J][K]) +
02579 dot(x[I ][J+1][K+1], x.get_mesh().Dvc[6]/vertSpacings[I][J][K]) +
02580 dot(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]/vertSpacings[I][J][K]);
02581 return r;
02582 }
02583
02584
02585
02586
02587 template < class T, class MFLOAT >
02588 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>&
02589 Div(Field<Tenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& x,
02590 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& r)
02591 {
02592 PAssert(x.get_mesh().hasSpacingFields);
02593 BareField<Vektor<MFLOAT,1U>,1U>& cellSpacings = *(x.get_mesh().CellSpacings);
02594 const NDIndex<1U>& domain = r.getDomain();
02595 Index I = domain[0];
02596 r[I] =
02597 dot(x[I-1], x.get_mesh().Dvc[0]/cellSpacings[I]) +
02598 dot(x[I ], x.get_mesh().Dvc[1]/cellSpacings[I]);
02599 return r;
02600 }
02601
02602 template < class T, class MFLOAT >
02603 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>&
02604 Div(Field<Tenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& x,
02605 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& r)
02606 {
02607 PAssert(x.get_mesh().hasSpacingFields);
02608 BareField<Vektor<MFLOAT,2U>,2U>& cellSpacings = *(x.get_mesh().CellSpacings);
02609 const NDIndex<2U>& domain = r.getDomain();
02610 Index I = domain[0];
02611 Index J = domain[1];
02612 r[I][J] =
02613 dot(x[I-1][J-1], x.get_mesh().Dvc[0]/cellSpacings[I][J]) +
02614 dot(x[I ][J-1], x.get_mesh().Dvc[1]/cellSpacings[I][J]) +
02615 dot(x[I-1][J ], x.get_mesh().Dvc[2]/cellSpacings[I][J]) +
02616 dot(x[I ][J ], x.get_mesh().Dvc[3]/cellSpacings[I][J]);
02617 return r;
02618 }
02619
02620 template < class T, class MFLOAT >
02621 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>&
02622 Div(Field<Tenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& x,
02623 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& r)
02624 {
02625 PAssert(x.get_mesh().hasSpacingFields);
02626 BareField<Vektor<MFLOAT,3U>,3U>& cellSpacings = *(x.get_mesh().CellSpacings);
02627 const NDIndex<3U>& domain = r.getDomain();
02628 Index I = domain[0];
02629 Index J = domain[1];
02630 Index K = domain[2];
02631 r[I][J][K] =
02632 dot(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]/cellSpacings[I][J][K]) +
02633 dot(x[I ][J-1][K-1], x.get_mesh().Dvc[1]/cellSpacings[I][J][K]) +
02634 dot(x[I-1][J ][K-1], x.get_mesh().Dvc[2]/cellSpacings[I][J][K]) +
02635 dot(x[I ][J ][K-1], x.get_mesh().Dvc[3]/cellSpacings[I][J][K]) +
02636 dot(x[I-1][J-1][K ], x.get_mesh().Dvc[4]/cellSpacings[I][J][K]) +
02637 dot(x[I ][J-1][K ], x.get_mesh().Dvc[5]/cellSpacings[I][J][K]) +
02638 dot(x[I-1][J ][K ], x.get_mesh().Dvc[6]/cellSpacings[I][J][K]) +
02639 dot(x[I ][J ][K ], x.get_mesh().Dvc[7]/cellSpacings[I][J][K]);
02640 return r;
02641 }
02642
02643
02644
02645
02646 template < class T, class MFLOAT >
02647 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>&
02648 Div(Field<SymTenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& x,
02649 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& r)
02650 {
02651 PAssert(x.get_mesh().hasSpacingFields);
02652 BareField<Vektor<MFLOAT,1U>,1U>& cellSpacings = *(x.get_mesh().CellSpacings);
02653 const NDIndex<1U>& domain = r.getDomain();
02654 Index I = domain[0];
02655 r[I] =
02656 dot(x[I-1], x.get_mesh().Dvc[0]/cellSpacings[I]) +
02657 dot(x[I ], x.get_mesh().Dvc[1]/cellSpacings[I]);
02658 return r;
02659 }
02660
02661 template < class T, class MFLOAT >
02662 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>&
02663 Div(Field<SymTenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& x,
02664 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& r)
02665 {
02666 PAssert(x.get_mesh().hasSpacingFields);
02667 BareField<Vektor<MFLOAT,2U>,2U>& cellSpacings = *(x.get_mesh().CellSpacings);
02668 const NDIndex<2U>& domain = r.getDomain();
02669 Index I = domain[0];
02670 Index J = domain[1];
02671 r[I][J] =
02672 dot(x[I-1][J-1], x.get_mesh().Dvc[0]/cellSpacings[I][J]) +
02673 dot(x[I ][J-1], x.get_mesh().Dvc[1]/cellSpacings[I][J]) +
02674 dot(x[I-1][J ], x.get_mesh().Dvc[2]/cellSpacings[I][J]) +
02675 dot(x[I ][J ], x.get_mesh().Dvc[3]/cellSpacings[I][J]);
02676 return r;
02677 }
02678
02679 template < class T, class MFLOAT >
02680 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>&
02681 Div(Field<SymTenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& x,
02682 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& r)
02683 {
02684 PAssert(x.get_mesh().hasSpacingFields);
02685 BareField<Vektor<MFLOAT,3U>,3U>& cellSpacings = *(x.get_mesh().CellSpacings);
02686 const NDIndex<3U>& domain = r.getDomain();
02687 Index I = domain[0];
02688 Index J = domain[1];
02689 Index K = domain[2];
02690 r[I][J][K] =
02691 dot(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]/cellSpacings[I][J][K]) +
02692 dot(x[I ][J-1][K-1], x.get_mesh().Dvc[1]/cellSpacings[I][J][K]) +
02693 dot(x[I-1][J ][K-1], x.get_mesh().Dvc[2]/cellSpacings[I][J][K]) +
02694 dot(x[I ][J ][K-1], x.get_mesh().Dvc[3]/cellSpacings[I][J][K]) +
02695 dot(x[I-1][J-1][K ], x.get_mesh().Dvc[4]/cellSpacings[I][J][K]) +
02696 dot(x[I ][J-1][K ], x.get_mesh().Dvc[5]/cellSpacings[I][J][K]) +
02697 dot(x[I-1][J ][K ], x.get_mesh().Dvc[6]/cellSpacings[I][J][K]) +
02698 dot(x[I ][J ][K ], x.get_mesh().Dvc[7]/cellSpacings[I][J][K]);
02699 return r;
02700 }
02701
02702
02703
02704
02705
02706
02707 template < class T, class MFLOAT >
02708 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>&
02709 Grad(Field<T,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02710 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02711 {
02712 PAssert(x.get_mesh().hasSpacingFields);
02713 BareField<Vektor<MFLOAT,1U>,1U>& vertSpacings = *(x.get_mesh().VertSpacings);
02714 const NDIndex<1U>& domain = r.getDomain();
02715 Index I = domain[0];
02716 r[I] = (x[I ]*x.get_mesh().Dvc[0] +
02717 x[I+1]*x.get_mesh().Dvc[1])/vertSpacings[I];
02718 return r;
02719 }
02720
02721 template < class T, class MFLOAT >
02722 Field<Vektor<T,2u>,2U,Cartesian<2U,MFLOAT>,Cell>&
02723 Grad(Field<T,2U,Cartesian<2U,MFLOAT>,Vert>& x,
02724 Field<Vektor<T,2u>,2U,Cartesian<2U,MFLOAT>,Cell>& r)
02725 {
02726 PAssert(x.get_mesh().hasSpacingFields);
02727 BareField<Vektor<MFLOAT,2U>,2U>& vertSpacings = *(x.get_mesh().VertSpacings);
02728 const NDIndex<2U>& domain = r.getDomain();
02729 Index I = domain[0];
02730 Index J = domain[1];
02731 r[I][J] = (x[I ][J ]*x.get_mesh().Dvc[0] +
02732 x[I+1][J ]*x.get_mesh().Dvc[1] +
02733 x[I ][J+1]*x.get_mesh().Dvc[2] +
02734 x[I+1][J+1]*x.get_mesh().Dvc[3])/vertSpacings[I][J];
02735 return r;
02736 }
02737
02738 template < class T, class MFLOAT >
02739 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>&
02740 Grad(Field<T,3U,Cartesian<3U,MFLOAT>,Vert>& x,
02741 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& r)
02742 {
02743 PAssert(x.get_mesh().hasSpacingFields);
02744 BareField<Vektor<MFLOAT,3U>,3U>& vertSpacings = *(x.get_mesh().VertSpacings);
02745 const NDIndex<3U>& domain = r.getDomain();
02746 Index I = domain[0];
02747 Index J = domain[1];
02748 Index K = domain[2];
02749 r[I][J][K] = (x[I ][J ][K ]*x.get_mesh().Dvc[0] +
02750 x[I+1][J ][K ]*x.get_mesh().Dvc[1] +
02751 x[I ][J+1][K ]*x.get_mesh().Dvc[2] +
02752 x[I+1][J+1][K ]*x.get_mesh().Dvc[3] +
02753 x[I ][J ][K+1]*x.get_mesh().Dvc[4] +
02754 x[I+1][J ][K+1]*x.get_mesh().Dvc[5] +
02755 x[I ][J+1][K+1]*x.get_mesh().Dvc[6] +
02756 x[I+1][J+1][K+1]*x.get_mesh().Dvc[7])/vertSpacings[I][J][K];
02757 return r;
02758 }
02759
02760
02761
02762
02763
02764 template < class T, class MFLOAT >
02765 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>&
02766 Grad(Field<T,1U,Cartesian<1U,MFLOAT>,Cell>& x,
02767 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& r)
02768 {
02769 PAssert(x.get_mesh().hasSpacingFields);
02770 BareField<Vektor<MFLOAT,1U>,1U>& cellSpacings = *(x.get_mesh().CellSpacings);
02771 const NDIndex<1U>& domain = r.getDomain();
02772 Index I = domain[0];
02773 r[I] = (x[I-1]*x.get_mesh().Dvc[0] +
02774 x[I ]*x.get_mesh().Dvc[1])/cellSpacings[I];
02775 return r;
02776 }
02777
02778 template < class T, class MFLOAT >
02779 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>&
02780 Grad(Field<T,2U,Cartesian<2U,MFLOAT>,Cell>& x,
02781 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& r)
02782 {
02783 PAssert(x.get_mesh().hasSpacingFields);
02784 BareField<Vektor<MFLOAT,2U>,2U>& cellSpacings = *(x.get_mesh().CellSpacings);
02785 const NDIndex<2U>& domain = r.getDomain();
02786 Index I = domain[0];
02787 Index J = domain[1];
02788 r[I][J] = (x[I-1][J-1]*x.get_mesh().Dvc[0] +
02789 x[I ][J-1]*x.get_mesh().Dvc[1] +
02790 x[I-1][J ]*x.get_mesh().Dvc[2] +
02791 x[I ][J ]*x.get_mesh().Dvc[3])/cellSpacings[I][J];
02792 return r;
02793 }
02794
02795 template < class T, class MFLOAT >
02796 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>&
02797 Grad(Field<T,3U,Cartesian<3U,MFLOAT>,Cell>& x,
02798 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& r)
02799 {
02800 PAssert(x.get_mesh().hasSpacingFields);
02801 BareField<Vektor<MFLOAT,3U>,3U>& cellSpacings = *(x.get_mesh().CellSpacings);
02802 const NDIndex<3U>& domain = r.getDomain();
02803 Index I = domain[0];
02804 Index J = domain[1];
02805 Index K = domain[2];
02806 Vektor<MFLOAT,3U> dvc[1<<3U];
02807 for (int d=0; d < 1<<3U; d++) dvc[d] = x.get_mesh().Dvc[d];
02808 r[I][J][K] = (x[I-1][J-1][K-1]*dvc[0] +
02809 x[I ][J-1][K-1]*dvc[1] +
02810 x[I-1][J ][K-1]*dvc[2] +
02811 x[I ][J ][K-1]*dvc[3] +
02812 x[I-1][J-1][K ]*dvc[4] +
02813 x[I ][J-1][K ]*dvc[5] +
02814 x[I-1][J ][K ]*dvc[6] +
02815 x[I ][J ][K ]*dvc[7])/
02816 cellSpacings[I][J][K];
02817 return r;
02818 }
02819
02820
02821
02822
02823
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835 template < class T, class MFLOAT >
02836 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>&
02837 Grad(Field<T,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02838 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& r)
02839 {
02840 PAssert(x.get_mesh().hasSpacingFields);
02841 BareField<Vektor<MFLOAT,1U>,1U>& vertSpacings =
02842 *(x.get_mesh().VertSpacings);
02843 const NDIndex<1U>& domain = r.getDomain();
02844 Index I = domain[0];
02845
02846 Vektor<MFLOAT,1U> idx;
02847 idx[0] = 1.0;
02848
02849 r[I] = idx*((x[I+1] - x[I-1])/(vertSpacings[I] + vertSpacings[I-1]));
02850 return r;
02851 }
02852
02853 template < class T, class MFLOAT >
02854 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>&
02855 Grad(Field<T,2U,Cartesian<2U,MFLOAT>,Vert>& x,
02856 Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& r)
02857 {
02858 PAssert(x.get_mesh().hasSpacingFields);
02859 BareField<Vektor<MFLOAT,2U>,2U>& vertSpacings =
02860 *(x.get_mesh().VertSpacings);
02861 const NDIndex<2U>& domain = r.getDomain();
02862 Index I = domain[0];
02863 Index J = domain[1];
02864
02865 Vektor<MFLOAT,2U> idx,idy;
02866 idx[0] = 1.0;
02867 idx[1] = 0.0;
02868 idy[0] = 0.0;
02869 idy[1] = 1.0;
02870
02871 r[I][J] =
02872 idx*((x[I+1][J ] - x[I-1][J ])/
02873 (vertSpacings[I][J] + vertSpacings[I-1][J])) +
02874 idy*((x[I ][J+1] - x[I ][J-1])/
02875 (vertSpacings[I][J] + vertSpacings[I][J-1]));
02876 return r;
02877 }
02878
02879 template < class T, class MFLOAT >
02880 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>&
02881 Grad(Field<T,3U,Cartesian<3U,MFLOAT>,Vert>& x,
02882 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& r)
02883 {
02884 PAssert(x.get_mesh().hasSpacingFields);
02885 BareField<Vektor<MFLOAT,3U>,3U>& vertSpacings =
02886 *(x.get_mesh().VertSpacings);
02887 const NDIndex<3U>& domain = r.getDomain();
02888 Index I = domain[0];
02889 Index J = domain[1];
02890 Index K = domain[2];
02891
02892 Vektor<MFLOAT,3U> idx,idy,idz;
02893 idx[0] = 1.0;
02894 idx[1] = 0.0;
02895 idx[2] = 0.0;
02896 idy[0] = 0.0;
02897 idy[1] = 1.0;
02898 idy[2] = 0.0;
02899 idz[0] = 0.0;
02900 idz[1] = 0.0;
02901 idz[2] = 1.0;
02902
02903 r[I][J][K] =
02904 idx*((x[I+1][J ][K ] - x[I-1][J ][K ])/
02905 (vertSpacings[I][J][K] + vertSpacings[I-1][J][K])) +
02906 idy*((x[I ][J+1][K ] - x[I ][J-1][K ])/
02907 (vertSpacings[I][J][K] + vertSpacings[I][J-1][K])) +
02908 idz*((x[I ][J ][K+1] - x[I ][J ][K-1])/
02909 (vertSpacings[I][J][K] + vertSpacings[I][J][K-1]));
02910 return r;
02911 }
02912
02913
02914
02915 template < class T, class MFLOAT >
02916 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>&
02917 Grad(Field<T,1U,Cartesian<1U,MFLOAT>,Cell>& x,
02918 Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02919 {
02920 PAssert(x.get_mesh().hasSpacingFields);
02921 BareField<Vektor<MFLOAT,1U>,1U>& cellSpacings =
02922 *(x.get_mesh().CellSpacings);
02923 const NDIndex<1U>& domain = r.getDomain();
02924 Index I = domain[0];
02925
02926 Vektor<MFLOAT,1U> idx;
02927 idx[0] = 1.0;
02928
02929 r[I] = idx*((x[I+1] - x[I-1])/(cellSpacings[I] + cellSpacings[I+1]));
02930 return r;
02931 }
02932
02933 template < class T, class MFLOAT >
02934 Field<Vektor<T,2u>,2U,Cartesian<2U,MFLOAT>,Cell>&
02935 Grad(Field<T,2U,Cartesian<2U,MFLOAT>,Cell>& x,
02936 Field<Vektor<T,2u>,2U,Cartesian<2U,MFLOAT>,Cell>& r)
02937 {
02938 PAssert(x.get_mesh().hasSpacingFields);
02939 BareField<Vektor<MFLOAT,2U>,2U>& cellSpacings =
02940 *(x.get_mesh().CellSpacings);
02941 const NDIndex<2U>& domain = r.getDomain();
02942 Index I = domain[0];
02943 Index J = domain[1];
02944
02945 Vektor<MFLOAT,2U> idx,idy;
02946 idx[0] = 1.0;
02947 idx[1] = 0.0;
02948 idy[0] = 0.0;
02949 idy[1] = 1.0;
02950
02951 r[I][J] =
02952 idx*((x[I+1][J ] - x[I-1][J ])/
02953 (cellSpacings[I][J] + cellSpacings[I+1][J])) +
02954 idy*((x[I ][J+1] - x[I ][J-1])/
02955 (cellSpacings[I][J] + cellSpacings[I][J+1]));
02956 return r;
02957 }
02958
02959 template < class T, class MFLOAT >
02960 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>&
02961 Grad(Field<T,3U,Cartesian<3U,MFLOAT>,Cell>& x,
02962 Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& r)
02963 {
02964 PAssert(x.get_mesh().hasSpacingFields);
02965 BareField<Vektor<MFLOAT,3U>,3U>& cellSpacings =
02966 *(x.get_mesh().CellSpacings);
02967 const NDIndex<3U>& domain = r.getDomain();
02968 Index I = domain[0];
02969 Index J = domain[1];
02970 Index K = domain[2];
02971
02972 Vektor<MFLOAT,3U> idx,idy,idz;
02973 idx[0] = 1.0;
02974 idx[1] = 0.0;
02975 idx[2] = 0.0;
02976 idy[0] = 0.0;
02977 idy[1] = 1.0;
02978 idy[2] = 0.0;
02979 idz[0] = 0.0;
02980 idz[1] = 0.0;
02981 idz[2] = 1.0;
02982
02983 r[I][J][K] =
02984 idx*((x[I+1][J ][K ] - x[I-1][J ][K ])/
02985 (cellSpacings[I][J][K] + cellSpacings[I+1][J][K])) +
02986 idy*((x[I ][J+1][K ] - x[I ][J-1][K ])/
02987 (cellSpacings[I][J][K] + cellSpacings[I][J+1][K])) +
02988 idz*((x[I ][J ][K+1] - x[I ][J ][K-1])/
02989 (cellSpacings[I][J][K] + cellSpacings[I][J][K+1]));
02990 return r;
02991 }
02992
02993
02994
02995 template < class T, class MFLOAT >
02996 Field<Tenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>&
02997 Grad(Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& x,
02998 Field<Tenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& r)
02999 {
03000 PAssert(x.get_mesh().hasSpacingFields);
03001 BareField<Vektor<MFLOAT,1U>,1U>& vS = *(x.get_mesh().VertSpacings);
03002 const NDIndex<1U>& domain = r.getDomain();
03003 Index I = domain[0];
03004 r[I] =
03005 outerProduct(x[I ], x.get_mesh().Dvc[0]/vS[I]) +
03006 outerProduct(x[I+1], x.get_mesh().Dvc[1]/vS[I]);
03007 return r;
03008 }
03009
03010 template < class T, class MFLOAT >
03011 Field<Tenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>&
03012 Grad(Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& x,
03013 Field<Tenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& r)
03014 {
03015 PAssert(x.get_mesh().hasSpacingFields);
03016 BareField<Vektor<MFLOAT,2U>,2U>& vS = *(x.get_mesh().VertSpacings);
03017 const NDIndex<2U>& domain = r.getDomain();
03018 Index I = domain[0];
03019 Index J = domain[1];
03020 r[I][J] =
03021 outerProduct(x[I ][J ], x.get_mesh().Dvc[0]/vS[I][J]) +
03022 outerProduct(x[I+1][J ], x.get_mesh().Dvc[1]/vS[I][J]) +
03023 outerProduct(x[I ][J+1], x.get_mesh().Dvc[2]/vS[I][J]) +
03024 outerProduct(x[I+1][J+1], x.get_mesh().Dvc[3]/vS[I][J]);
03025 return r;
03026 }
03027
03028 template < class T, class MFLOAT >
03029 Field<Tenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>&
03030 Grad(Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& x,
03031 Field<Tenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& r)
03032 {
03033 PAssert(x.get_mesh().hasSpacingFields);
03034 BareField<Vektor<MFLOAT,3U>,3U>& vS = *(x.get_mesh().VertSpacings);
03035 const NDIndex<3U>& domain = r.getDomain();
03036 Index I = domain[0];
03037 Index J = domain[1];
03038 Index K = domain[2];
03039 r[I][J][K] =
03040 outerProduct(x[I ][J ][K ], x.get_mesh().Dvc[0]/vS[I][J][K]) +
03041 outerProduct(x[I+1][J ][K ], x.get_mesh().Dvc[1]/vS[I][J][K]) +
03042 outerProduct(x[I ][J+1][K ], x.get_mesh().Dvc[2]/vS[I][J][K]) +
03043 outerProduct(x[I+1][J+1][K ], x.get_mesh().Dvc[3]/vS[I][J][K]) +
03044 outerProduct(x[I ][J ][K+1], x.get_mesh().Dvc[4]/vS[I][J][K]) +
03045 outerProduct(x[I+1][J ][K+1], x.get_mesh().Dvc[5]/vS[I][J][K]) +
03046 outerProduct(x[I ][J+1][K+1], x.get_mesh().Dvc[6]/vS[I][J][K]) +
03047 outerProduct(x[I+1][J+1][K+1], x.get_mesh().Dvc[7]/vS[I][J][K]);
03048
03049 return r;
03050 }
03051
03052
03053
03054
03055 template < class T, class MFLOAT >
03056 Field<Tenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>&
03057 Grad(Field<Vektor<T,1U>,1U,Cartesian<1U,MFLOAT>,Cell>& x,
03058 Field<Tenzor<T,1U>,1U,Cartesian<1U,MFLOAT>,Vert>& r)
03059 {
03060 PAssert(x.get_mesh().hasSpacingFields);
03061 BareField<Vektor<MFLOAT,1U>,1U>& cellSpacings = *(x.get_mesh().CellSpacings);
03062 const NDIndex<1U>& domain = r.getDomain();
03063 Index I = domain[0];
03064 r[I] =
03065 outerProduct(x[I-1], x.get_mesh().Dvc[0]/cellSpacings[I-1]) +
03066 outerProduct(x[I ], x.get_mesh().Dvc[1]/cellSpacings[I-1]);
03067 return r;
03068 }
03069
03070 template < class T, class MFLOAT >
03071 Field<Tenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>&
03072 Grad(Field<Vektor<T,2U>,2U,Cartesian<2U,MFLOAT>,Cell>& x,
03073 Field<Tenzor<T,2U>,2U,Cartesian<2U,MFLOAT>,Vert>& r)
03074 {
03075 PAssert(x.get_mesh().hasSpacingFields);
03076 BareField<Vektor<MFLOAT,2U>,2U>& cS = *(x.get_mesh().CellSpacings);
03077 const NDIndex<2U>& domain = r.getDomain();
03078 Index I = domain[0];
03079 Index J = domain[1];
03080 r[I][J] =
03081 outerProduct(x[I-1][J-1], x.get_mesh().Dvc[0]/cS[I-1][J-1]) +
03082 outerProduct(x[I ][J-1], x.get_mesh().Dvc[1]/cS[I-1][J-1]) +
03083 outerProduct(x[I-1][J ], x.get_mesh().Dvc[2]/cS[I-1][J-1]) +
03084 outerProduct(x[I ][J ], x.get_mesh().Dvc[3]/cS[I-1][J-1]);
03085 return r;
03086 }
03087
03088 template < class T, class MFLOAT >
03089 Field<Tenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>&
03090 Grad(Field<Vektor<T,3U>,3U,Cartesian<3U,MFLOAT>,Cell>& x,
03091 Field<Tenzor<T,3U>,3U,Cartesian<3U,MFLOAT>,Vert>& r)
03092 {
03093 PAssert(x.get_mesh().hasSpacingFields);
03094 BareField<Vektor<MFLOAT,3U>,3U>& cS = *(x.get_mesh().CellSpacings);
03095 const NDIndex<3U>& domain = r.getDomain();
03096 Index I = domain[0];
03097 Index J = domain[1];
03098 Index K = domain[2];
03099 r[I][J][K] =
03100 outerProduct(x[I-1][J-1][K-1], x.get_mesh().Dvc[0]/cS[I-1][J-1][K-1]) +
03101 outerProduct(x[I ][J-1][K-1], x.get_mesh().Dvc[1]/cS[I-1][J-1][K-1]) +
03102 outerProduct(x[I-1][J ][K-1], x.get_mesh().Dvc[2]/cS[I-1][J-1][K-1]) +
03103 outerProduct(x[I ][J ][K-1], x.get_mesh().Dvc[3]/cS[I-1][J-1][K-1]) +
03104 outerProduct(x[I-1][J-1][K ], x.get_mesh().Dvc[4]/cS[I-1][J-1][K-1]) +
03105 outerProduct(x[I ][J-1][K ], x.get_mesh().Dvc[5]/cS[I-1][J-1][K-1]) +
03106 outerProduct(x[I-1][J ][K ], x.get_mesh().Dvc[6]/cS[I-1][J-1][K-1]) +
03107 outerProduct(x[I ][J ][K ], x.get_mesh().Dvc[7]/cS[I-1][J-1][K-1]);
03108 return r;
03109 }
03110
03111
03112
03113
03114
03115
03116 template < class T1, class T2, class MFLOAT >
03117 Field<T1,1U,Cartesian<1U,MFLOAT>,Vert>&
03118 Average(Field<T1,1U,Cartesian<1U,MFLOAT>,Cell>& x,
03119 Field<T2,1U,Cartesian<1U,MFLOAT>,Cell>& w,
03120 Field<T1,1U,Cartesian<1U,MFLOAT>,Vert>& r)
03121 {
03122 const NDIndex<1U>& domain = r.getDomain();
03123 Index I = domain[0];
03124 r[I] = (x[I-1]*w[I-1] + x[I ]*w[I ])/(w[I-1] + w[I ]);
03125 return r;
03126 }
03127
03128 template < class T1, class T2, class MFLOAT >
03129 Field<T1,2U,Cartesian<2U,MFLOAT>,Vert>&
03130 Average(Field<T1,2U,Cartesian<2U,MFLOAT>,Cell>& x,
03131 Field<T2,2U,Cartesian<2U,MFLOAT>,Cell>& w,
03132 Field<T1,2U,Cartesian<2U,MFLOAT>,Vert>& r)
03133 {
03134 const NDIndex<2U>& domain = r.getDomain();
03135 Index I = domain[0];
03136 Index J = domain[1];
03137
03138 r[I][J] = (x[I-1][J-1] * w[I-1][J-1] +
03139 x[I ][J-1] * w[I ][J-1] +
03140 x[I-1][J ] * w[I-1][J ] +
03141 x[I ][J ] * w[I ][J ])/
03142 (w[I-1][J-1] +
03143 w[I ][J-1] +
03144 w[I-1][J ] +
03145 w[I ][J ]);
03146 return r;
03147 }
03148
03149 template < class T1, class T2, class MFLOAT >
03150 Field<T1,3U,Cartesian<3U,MFLOAT>,Vert>&
03151 Average(Field<T1,3U,Cartesian<3U,MFLOAT>,Cell>& x,
03152 Field<T2,3U,Cartesian<3U,MFLOAT>,Cell>& w,
03153 Field<T1,3U,Cartesian<3U,MFLOAT>,Vert>& r)
03154 {
03155 const NDIndex<3U>& domain = r.getDomain();
03156 Index I = domain[0];
03157 Index J = domain[1];
03158 Index K = domain[2];
03159
03160 r[I][J][K] = (x[I-1][J-1][K-1] * w[I-1][J-1][K-1] +
03161 x[I ][J-1][K-1] * w[I ][J-1][K-1] +
03162 x[I-1][J ][K-1] * w[I-1][J ][K-1] +
03163 x[I ][J ][K-1] * w[I ][J ][K-1] +
03164 x[I-1][J-1][K ] * w[I-1][J-1][K ] +
03165 x[I ][J-1][K ] * w[I ][J-1][K ] +
03166 x[I-1][J ][K ] * w[I-1][J ][K ] +
03167 x[I ][J ][K ] * w[I ][J ][K ] )/
03168 (w[I-1][J-1][K-1] +
03169 w[I ][J-1][K-1] +
03170 w[I-1][J ][K-1] +
03171 w[I ][J ][K-1] +
03172 w[I-1][J-1][K ] +
03173 w[I ][J-1][K ] +
03174 w[I-1][J ][K ] +
03175 w[I ][J ][K ]);
03176 return r;
03177 }
03178
03179
03180
03181
03182 template < class T1, class T2, class MFLOAT >
03183 Field<T1,1U,Cartesian<1U,MFLOAT>,Cell>&
03184 Average(Field<T1,1U,Cartesian<1U,MFLOAT>,Vert>& x,
03185 Field<T2,1U,Cartesian<1U,MFLOAT>,Vert>& w,
03186 Field<T1,1U,Cartesian<1U,MFLOAT>,Cell>& r)
03187 {
03188 const NDIndex<1U>& domain = r.getDomain();
03189 Index I = domain[0];
03190 r[I] = (x[I+1]*w[I+1] + x[I ]*w[I ])/(w[I+1] + w[I ]);
03191 return r;
03192 }
03193
03194 template < class T1, class T2, class MFLOAT >
03195 Field<T1,2U,Cartesian<2U,MFLOAT>,Cell>&
03196 Average(Field<T1,2U,Cartesian<2U,MFLOAT>,Vert>& x,
03197 Field<T2,2U,Cartesian<2U,MFLOAT>,Vert>& w,
03198 Field<T1,2U,Cartesian<2U,MFLOAT>,Cell>& r)
03199 {
03200 const NDIndex<2U>& domain = r.getDomain();
03201 Index I = domain[0];
03202 Index J = domain[1];
03203
03204 r[I][J] = (x[I+1][J+1] * w[I+1][J+1] +
03205 x[I ][J+1] * w[I ][J+1] +
03206 x[I+1][J ] * w[I+1][J ] +
03207 x[I ][J ] * w[I ][J ])/
03208 (w[I+1][J+1] +
03209 w[I ][J+1] +
03210 w[I+1][J ] +
03211 w[I ][J ]);
03212 return r;
03213 }
03214
03215 template < class T1, class T2, class MFLOAT >
03216 Field<T1,3U,Cartesian<3U,MFLOAT>,Cell>&
03217 Average(Field<T1,3U,Cartesian<3U,MFLOAT>,Vert>& x,
03218 Field<T2,3U,Cartesian<3U,MFLOAT>,Vert>& w,
03219 Field<T1,3U,Cartesian<3U,MFLOAT>,Cell>& r)
03220 {
03221 const NDIndex<3U>& domain = r.getDomain();
03222 Index I = domain[0];
03223 Index J = domain[1];
03224 Index K = domain[2];
03225
03226 r[I][J][K] = (x[I+1][J+1][K+1] * w[I+1][J+1][K+1] +
03227 x[I ][J+1][K+1] * w[I ][J+1][K+1] +
03228 x[I+1][J ][K+1] * w[I+1][J ][K+1] +
03229 x[I ][J ][K+1] * w[I ][J ][K+1] +
03230 x[I+1][J+1][K ] * w[I+1][J+1][K ] +
03231 x[I ][J+1][K ] * w[I ][J+1][K ] +
03232 x[I+1][J ][K ] * w[I+1][J ][K ] +
03233 x[I ][J ][K ] * w[I ][J ][K ])/
03234 (w[I+1][J+1][K+1] +
03235 w[I ][J+1][K+1] +
03236 w[I+1][J ][K+1] +
03237 w[I ][J ][K+1] +
03238 w[I+1][J+1][K ] +
03239 w[I ][J+1][K ] +
03240 w[I+1][J ][K ] +
03241 w[I ][J ][K ]);
03242 return r;
03243 }
03244
03245
03246
03247
03248 template < class T1, class MFLOAT >
03249 Field<T1,1U,Cartesian<1U,MFLOAT>,Vert>&
03250 Average(Field<T1,1U,Cartesian<1U,MFLOAT>,Cell>& x,
03251 Field<T1,1U,Cartesian<1U,MFLOAT>,Vert>& r)
03252 {
03253 const NDIndex<1U>& domain = r.getDomain();
03254 Index I = domain[0];
03255 r[I] = 0.5*(x[I-1] + x[I ]);
03256 return r;
03257 }
03258
03259 template < class T1, class MFLOAT >
03260 Field<T1,2U,Cartesian<2U,MFLOAT>,Vert>&
03261 Average(Field<T1,2U,Cartesian<2U,MFLOAT>,Cell>& x,
03262 Field<T1,2U,Cartesian<2U,MFLOAT>,Vert>& r)
03263 {
03264 const NDIndex<2U>& domain = r.getDomain();
03265 Index I = domain[0];
03266 Index J = domain[1];
03267 r[I][J] = 0.25*(x[I-1][J-1] + x[I ][J-1] + x[I-1][J ] + x[I ][J ]);
03268 return r;
03269 }
03270
03271 template < class T1, class MFLOAT >
03272 Field<T1,3U,Cartesian<3U,MFLOAT>,Vert>&
03273 Average(Field<T1,3U,Cartesian<3U,MFLOAT>,Cell>& x,
03274 Field<T1,3U,Cartesian<3U,MFLOAT>,Vert>& r)
03275 {
03276 const NDIndex<3U>& domain = r.getDomain();
03277 Index I = domain[0];
03278 Index J = domain[1];
03279 Index K = domain[2];
03280 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] +
03281 x[I ][J ][K-1] + x[I-1][J-1][K ] + x[I ][J-1][K ] +
03282 x[I-1][J ][K ] + x[I ][J ][K ]);
03283 return r;
03284 }
03285
03286
03287
03288
03289 template < class T1, class MFLOAT >
03290 Field<T1,1U,Cartesian<1U,MFLOAT>,Cell>&
03291 Average(Field<T1,1U,Cartesian<1U,MFLOAT>,Vert>& x,
03292 Field<T1,1U,Cartesian<1U,MFLOAT>,Cell>& r)
03293 {
03294 const NDIndex<1U>& domain = r.getDomain();
03295 Index I = domain[0];
03296 r[I] = 0.5*(x[I+1] + x[I ]);
03297 return r;
03298 }
03299
03300 template < class T1, class MFLOAT >
03301 Field<T1,2U,Cartesian<2U,MFLOAT>,Cell>&
03302 Average(Field<T1,2U,Cartesian<2U,MFLOAT>,Vert>& x,
03303 Field<T1,2U,Cartesian<2U,MFLOAT>,Cell>& r)
03304 {
03305 const NDIndex<2U>& domain = r.getDomain();
03306 Index I = domain[0];
03307 Index J = domain[1];
03308 r[I][J] = 0.25*(x[I+1][J+1] + x[I ][J+1] + x[I+1][J ] + x[I ][J ]);
03309 return r;
03310 }
03311
03312 template < class T1, class MFLOAT >
03313 Field<T1,3U,Cartesian<3U,MFLOAT>,Cell>&
03314 Average(Field<T1,3U,Cartesian<3U,MFLOAT>,Vert>& x,
03315 Field<T1,3U,Cartesian<3U,MFLOAT>,Cell>& r)
03316 {
03317 const NDIndex<3U>& domain = r.getDomain();
03318 Index I = domain[0];
03319 Index J = domain[1];
03320 Index K = domain[2];
03321 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] +
03322 x[I ][J ][K+1] + x[I+1][J+1][K ] + x[I ][J+1][K ] +
03323 x[I+1][J ][K ] + x[I ][J ][K ]);
03324 return r;
03325 }
03326
03327
03328
03329
03330