00001 #ifdef CH_LANG_CC
00002
00003
00004
00005
00006
00007
00008
00009 #endif
00010
00011 #ifndef _PROBLEMDOMAIN_H_
00012 #define _PROBLEMDOMAIN_H_
00013
00014 #ifndef WRAPPER
00015 #include <iostream>
00016
00017 #include "Vector.H"
00018 #include "IntVect.H"
00019 #include "Box.H"
00020 #include "Misc.H"
00021 #endif
00022
00023 #include "SPACE.H"
00024 #include "NamespaceHeader.H"
00025
00026
00027
00028
00029
00030
00031
00032
00033 class ShiftIterator
00034 {
00035 public:
00036
00037 ShiftIterator();
00038
00039
00040
00041
00042
00043
00044 ShiftIterator(const bool* a_isPeriodic);
00045
00046
00047 ShiftIterator(const ShiftIterator& a_shiftIt);
00048
00049
00050 ~ShiftIterator();
00051
00052
00053 ShiftIterator& operator=(const ShiftIterator& a_src);
00054
00055
00056 void computeShifts(const bool* a_isPeriodic);
00057
00058
00059 inline IntVect operator()() const;
00060
00061
00062 IntVect i() const
00063 {
00064 return this->operator()();
00065 }
00066
00067 inline void operator++();
00068
00069
00070 void incr()
00071 {
00072 ++(*this);
00073 }
00074
00075
00076 inline bool ok() const;
00077
00078
00079 int index() const
00080 {
00081 return m_index;
00082 }
00083
00084 const IntVect& operator[](int index) const
00085 {
00086 return m_shift_vectors[index];
00087 }
00088
00089 inline void reset();
00090
00091
00092 inline void begin();
00093
00094
00095
00096
00097
00098 void end();
00099
00100 private:
00101
00102 int m_index;
00103
00104 Vector<IntVect> m_shift_vectors;
00105
00106 };
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 class ProblemDomain
00131 {
00132 public:
00133
00134
00135
00136
00137
00138
00139 ProblemDomain ();
00140
00141
00142
00143
00144
00145 ProblemDomain(const Box& a_domBox);
00146
00147
00148
00149
00150
00151
00152
00153 ProblemDomain(const Box& a_domBox, const bool* a_isPeriodic);
00154
00155
00156
00157
00158
00159
00160 ProblemDomain (const IntVect& small,
00161 const IntVect& big);
00162
00163
00164
00165
00166
00167
00168
00169
00170 ProblemDomain (const IntVect& small,
00171 const IntVect& big,
00172 const bool* a_isPeriodic);
00173
00174
00175
00176
00177
00178
00179 ProblemDomain (const IntVect& small,
00180 const int* vec_len);
00181
00182
00183
00184
00185
00186
00187
00188
00189 ProblemDomain (const IntVect& small,
00190 const int* vec_len,
00191 const bool* a_isPeriodic);
00192
00193
00194
00195
00196 ProblemDomain (const ProblemDomain& a_src);
00197
00198
00199
00200
00201
00202 void define(const Box& a_domBox);
00203
00204
00205
00206
00207
00208
00209
00210 void define(const Box& a_domBox, const bool* a_isPeriodic);
00211
00212
00213
00214
00215
00216
00217 void define (const IntVect& small,
00218 const IntVect& big);
00219
00220
00221
00222
00223
00224
00225
00226
00227 void define (const IntVect& small,
00228 const IntVect& big,
00229 const bool* a_isPeriodic);
00230
00231
00232
00233
00234
00235
00236 void define (const IntVect& small,
00237 const int* vec_len);
00238
00239
00240
00241
00242
00243
00244
00245
00246 void define (const IntVect& small,
00247 const int* vec_len,
00248 const bool* a_isPeriodic);
00249
00250
00251
00252
00253 void define (const ProblemDomain& a_src);
00254
00255
00256
00257
00258
00259
00260 const Box& domainBox() const;
00261
00262
00263
00264
00265 bool isPeriodic(int a_dir) const;
00266
00267
00268
00269
00270 bool isPeriodic() const;
00271
00272
00273
00274
00275
00276
00277 ShiftIterator shiftIterator() const;
00278
00279
00280
00281
00282 bool isEmpty () const;
00283
00284
00285
00286
00287 int size (const int& a_idir) const;
00288
00289
00290
00291
00292 IntVect size () const;
00293
00294
00295
00296
00297
00298
00299
00300
00301 bool contains (const IntVect& p) const;
00302
00303
00304
00305
00306
00307
00308 bool image (IntVect& p) const;
00309
00310
00311
00312
00313
00314
00315 bool contains (const Box& b) const;
00316
00317
00318 bool contains_box(const Box& b) const
00319 {
00320 return contains(b);
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 bool intersects (const Box& a_box) const;
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342 bool intersectsNotEmpty (const Box& a_box) const;
00343
00344
00345
00346
00347 bool periodicAdjacent(const Box& a_box) const;
00348
00349
00350
00351
00352 void insertImages(std::list<Box>& a_list, const Box& a_box) const;
00353
00354
00355
00356
00357
00358
00359 bool intersects(const Box& box1, const Box& box2) const;
00360
00361
00362
00363
00364 bool operator== (const ProblemDomain& a_otherDomain) const;
00365
00366
00367
00368
00369 bool operator!= (const ProblemDomain& a_otherDomain) const;
00370
00371
00372
00373
00374 friend void operator &=(Box& a_box, const ProblemDomain& a_probdomain);
00375
00376
00377
00378
00379 friend Box operator & (const Box& a_box, const ProblemDomain& a_probdomain);
00380
00381
00382
00383
00384
00385
00386 ProblemDomain& operator= (const ProblemDomain& b);
00387
00388
00389
00390
00391 void setPeriodic(int a_dir, bool a_isPeriodic);
00392
00393
00394
00395
00396 inline ProblemDomain& grow (int i);
00397
00398
00399
00400
00401 friend inline ProblemDomain grow(const ProblemDomain& pd,
00402 int i);
00403
00404
00405
00406
00407
00408
00409 inline ProblemDomain& grow(const IntVect& v);
00410
00411
00412
00413
00414
00415
00416 friend inline ProblemDomain grow (const ProblemDomain& pd,
00417 const IntVect& v);
00418
00419
00420
00421
00422
00423
00424 inline ProblemDomain& grow(int idir, int n_cell);
00425
00426
00427
00428
00429
00430
00431 inline ProblemDomain& growLo(int idir, int n_cell=1);
00432
00433
00434
00435
00436
00437 inline ProblemDomain& growHi(int idir, int n_cell=1);
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448 friend Box bdryLo (const ProblemDomain& a_pd,
00449 int a_dir,
00450 int a_len=1);
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461 friend Box bdryHi (const ProblemDomain& a_pd,
00462 int a_dir,
00463 int a_len=1);
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491 friend Box adjCellLo (const ProblemDomain& a_pd,
00492 int a_dir,
00493 int a_len=1);
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522 friend Box adjCellHi (const ProblemDomain& a_pd,
00523 int a_dir,
00524 int a_len=1);
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 Box operator& (const Box& a_b) const;
00536
00537
00538
00539
00540
00541
00542
00543
00544 ProblemDomain& refine (int a_refinement_ratio);
00545
00546
00547
00548
00549
00550
00551
00552 friend ProblemDomain refine (const ProblemDomain& a_probdomain,
00553 int a_refinement_ratio);
00554
00555
00556
00557
00558
00559
00560 ProblemDomain& refine (const IntVect& a_refinement_ratio);
00561
00562
00563
00564
00565
00566
00567
00568
00569 friend ProblemDomain refine (const ProblemDomain& a_probdomain,
00570 const IntVect& a_refinement_ratio);
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 ProblemDomain& coarsen (int a_refinement_ratio);
00582
00583
00584
00585
00586
00587
00588
00589 friend ProblemDomain coarsen (const ProblemDomain& a_probdomain,
00590 int a_refinement_ratio);
00591
00592
00593
00594
00595
00596
00597 ProblemDomain& coarsen (const IntVect& refinement_ratio);
00598
00599
00600
00601
00602
00603
00604
00605
00606 friend ProblemDomain coarsen (const ProblemDomain& a_probdomain,
00607 const IntVect& a_refinement_ratio);
00608
00609
00610
00611 void shift(const IntVect& a_shift)
00612 {
00613 m_domainBox.shift(a_shift);
00614 }
00615
00616 bool operator<(const ProblemDomain& rhs) const
00617 {
00618 return m_domainBox < rhs.m_domainBox;
00619 }
00620
00621
00622
00623
00624
00625 friend std::ostream& operator<< (std::ostream& os,
00626 const ProblemDomain& bx);
00627
00628
00629
00630
00631 friend std::istream& operator>> (std::istream& is,
00632 ProblemDomain& bx);
00633
00634 void shiftIt(Box& a_box, int shiftIndex) const ;
00635 void unshiftIt(Box& a_box, int shiftIndex) const ;
00636
00637
00638
00639
00640 void dumpOn (std::ostream& strm) const;
00641
00642
00643 protected:
00644 friend class HDF5Handle;
00645
00646
00647
00648
00649 bool m_isPeriodic[SpaceDim];
00650
00651
00652
00653
00654 Box m_domainBox;
00655
00656
00657
00658
00659 ShiftIterator m_shiftIt;
00660 };
00661
00662 class ImageIterator
00663 {
00664 public:
00665 ImageIterator(const ProblemDomain& a_domain)
00666 {
00667 define(a_domain);
00668 }
00669
00670 void define(const ProblemDomain& a_domain);
00671
00672 void begin(const Box& a_box)
00673 {
00674 m_counter=-1;
00675 m_box = a_box;
00676 this->operator++();
00677 }
00678
00679 void operator++();
00680
00681 bool ok()
00682 {
00683 return m_shifter[m_counter] != IntVect::Zero;
00684 }
00685
00686 const Box& box() const
00687 {
00688 return m_current;
00689 }
00690
00691 const ProblemDomain& domain() const
00692 {
00693 return m_domain;
00694 }
00695
00696 void checkDefine(const ProblemDomain& a_domain)
00697 {
00698 if (!(m_domain == a_domain)) define(a_domain);
00699 }
00700
00701 protected:
00702 ProblemDomain m_domain;
00703 Box m_quadrant[D_TERM6(3,*3,*3,*3,*3,*3)];
00704 IntVect m_shifter[D_TERM6(3,*3,*3,*3,*3,*3)];
00705 Box m_box;
00706 Box m_current;
00707 int m_counter;
00708 };
00709
00710
00711
00712
00713 #ifndef WRAPPER
00714
00715 inline
00716 ShiftIterator::ShiftIterator()
00717 : m_index(100), m_shift_vectors()
00718 {
00719 }
00720
00721 inline
00722 ShiftIterator::ShiftIterator(const ShiftIterator& a_src)
00723 {
00724 m_index = a_src.m_index;
00725 m_shift_vectors = a_src.m_shift_vectors;
00726 }
00727
00728 inline
00729 ShiftIterator&
00730 ShiftIterator::operator=(const ShiftIterator& a_src)
00731 {
00732 m_index = a_src.m_index;
00733 m_shift_vectors = a_src.m_shift_vectors;
00734 return *this;
00735 }
00736
00737 inline
00738 IntVect
00739 ShiftIterator::operator()() const
00740 {
00741 CH_assert(ok());
00742 return m_shift_vectors[m_index];
00743 }
00744
00745 inline
00746 void
00747 ShiftIterator::operator++()
00748 {
00749 m_index++;
00750 }
00751
00752 inline
00753 bool
00754 ShiftIterator::ok() const
00755 {
00756 return (m_index < m_shift_vectors.size());
00757 }
00758
00759 inline
00760 void
00761 ShiftIterator::reset()
00762 {
00763 m_index = 0;
00764 }
00765
00766 inline
00767 void
00768 ShiftIterator::begin()
00769 {
00770 m_index = 0;
00771 }
00772
00773 inline
00774 void
00775 ShiftIterator::end()
00776 {
00777 m_index = m_shift_vectors.size();
00778 }
00779
00780 inline
00781 ProblemDomain::ProblemDomain ()
00782 {
00783
00784 for (int dir=0; dir<SpaceDim; dir++)
00785 {
00786 m_isPeriodic[dir] = false;
00787 }
00788
00789 }
00790
00791 inline
00792 ProblemDomain::ProblemDomain (const ProblemDomain& b)
00793 : m_domainBox(b.m_domainBox), m_shiftIt(b.m_shiftIt)
00794 {
00795 for (int dir=0; dir<SpaceDim; dir++)
00796 {
00797 m_isPeriodic[dir] = b.m_isPeriodic[dir];
00798 }
00799 }
00800
00801 inline void
00802 ProblemDomain::define(const ProblemDomain& b)
00803
00804 {
00805 m_domainBox=(b.m_domainBox);
00806 m_shiftIt = (b.m_shiftIt);
00807 for (int dir=0; dir<SpaceDim; dir++)
00808 {
00809 m_isPeriodic[dir] = b.m_isPeriodic[dir];
00810 }
00811 }
00812
00813 inline
00814 bool
00815 ProblemDomain::operator== (const ProblemDomain& a_otherDomain) const
00816 {
00817 bool result = true;
00818
00819 if (m_domainBox != a_otherDomain.m_domainBox)
00820 {
00821 result = false;
00822 }
00823 else
00824 {
00825 for (int dir=0; dir<SpaceDim; dir++)
00826 {
00827 if (m_isPeriodic[dir] != a_otherDomain.m_isPeriodic[dir])
00828 {
00829 result = false;
00830 break;
00831 }
00832 }
00833 }
00834
00835 return result;
00836 }
00837
00838 inline
00839 bool
00840 ProblemDomain::operator!= (const ProblemDomain& a_otherDomain) const
00841 {
00842 return !(*this == a_otherDomain);
00843 }
00844
00845 inline
00846 ProblemDomain&
00847 ProblemDomain::operator= (const ProblemDomain& b)
00848 {
00849 m_domainBox = b.m_domainBox;
00850 for (int dir=0; dir<SpaceDim; dir++)
00851 {
00852 m_isPeriodic[dir] = b.m_isPeriodic[dir];
00853 }
00854 m_shiftIt = b.m_shiftIt;
00855 return *this;
00856 }
00857
00858 inline void
00859 ProblemDomain::shiftIt(Box& a_box, int a_shiftIndex) const
00860 {
00861 a_box.shift(m_shiftIt[a_shiftIndex]*m_domainBox.size());
00862 }
00863
00864 inline void
00865 ProblemDomain::unshiftIt(Box& a_box, int a_shiftIndex) const
00866 {
00867 a_box.shift(- m_shiftIt[a_shiftIndex]*m_domainBox.size());
00868 }
00869
00870
00871 inline
00872 const Box&
00873 ProblemDomain::domainBox() const
00874 {
00875 return m_domainBox;
00876 }
00877
00878 inline
00879 bool
00880 ProblemDomain::isPeriodic(int a_dir) const
00881 {
00882 return m_isPeriodic[a_dir];
00883 }
00884
00885 inline
00886 bool
00887 ProblemDomain::isPeriodic() const
00888 {
00889 return D_TERM6(m_isPeriodic[0], ||
00890 m_isPeriodic[1], ||
00891 m_isPeriodic[2], ||
00892 m_isPeriodic[3], ||
00893 m_isPeriodic[4], ||
00894 m_isPeriodic[5]);
00895 }
00896
00897 inline
00898 ProblemDomain&
00899 ProblemDomain::grow(int i)
00900 {
00901 m_domainBox.grow(i);
00902 return *this;
00903 }
00904
00905 inline
00906 ProblemDomain
00907 grow(const ProblemDomain& pd, int i)
00908 {
00909 ProblemDomain newPd(pd);
00910 newPd.grow(i);
00911 return newPd;
00912 }
00913
00914 inline
00915 ProblemDomain&
00916 ProblemDomain::grow(const IntVect& v)
00917 {
00918 m_domainBox.grow(v);
00919 return *this;
00920 }
00921
00922 inline
00923 ProblemDomain
00924 grow(const ProblemDomain& pd, const IntVect& v)
00925 {
00926 ProblemDomain newPd(pd);
00927 newPd.grow(v);
00928 return newPd;
00929 }
00930
00931 inline
00932 ProblemDomain&
00933 ProblemDomain::grow(int idir, int n_cell)
00934 {
00935 m_domainBox.grow(idir, n_cell);
00936 return *this;
00937 }
00938
00939 inline
00940 ProblemDomain&
00941 ProblemDomain::growLo(int idir, int n_cell)
00942 {
00943 m_domainBox.growLo(idir, n_cell);
00944 return *this;
00945 }
00946
00947 inline
00948 ProblemDomain&
00949 ProblemDomain::growHi(int idir, int n_cell)
00950 {
00951 m_domainBox.growHi(idir, n_cell);
00952 return *this;
00953 }
00954
00955 inline
00956 ShiftIterator
00957 ProblemDomain::shiftIterator() const
00958 {
00959 return m_shiftIt;
00960 }
00961
00962 inline
00963 bool
00964 ProblemDomain::isEmpty () const
00965 {
00966 return (m_domainBox.isEmpty());
00967 }
00968
00969 inline
00970 int
00971 ProblemDomain::size (const int& a_idir) const
00972 {
00973 return (m_domainBox.size(a_idir));
00974 }
00975
00976 inline
00977 IntVect
00978 ProblemDomain::size () const
00979 {
00980 return (m_domainBox.size());
00981 }
00982
00983 inline
00984 bool
00985 ProblemDomain::contains (const IntVect& p) const
00986 {
00987
00988
00989 return ( !isEmpty()
00990 && (D_TERM6((m_isPeriodic[0]
00991 || ( p[0] >= m_domainBox.smallEnd(0)
00992 && p[0] <= m_domainBox.bigEnd (0))),
00993 && (m_isPeriodic[1]
00994 || ( p[1] >= m_domainBox.smallEnd(1)
00995 && p[1] <= m_domainBox.bigEnd (1))),
00996 && (m_isPeriodic[2]
00997 || ( p[2] >= m_domainBox.smallEnd(2)
00998 && p[2] <= m_domainBox.bigEnd (2))),
00999 && (m_isPeriodic[3]
01000 || ( p[3] >= m_domainBox.smallEnd(3)
01001 && p[3] <= m_domainBox.bigEnd (3))),
01002 && (m_isPeriodic[4]
01003 || ( p[4] >= m_domainBox.smallEnd(4)
01004 && p[4] <= m_domainBox.bigEnd (4))),
01005 && (m_isPeriodic[5]
01006 || ( p[5] >= m_domainBox.smallEnd(5)
01007 && p[5] <= m_domainBox.bigEnd (5))))));
01008 }
01009
01010 inline
01011 bool
01012 ProblemDomain::image(IntVect& p) const
01013 {
01014 if (m_domainBox.contains(p)) return true;
01015 if (!contains(p)) return false;
01016
01017 D_TERM6(
01018 if (m_isPeriodic[0])
01019 {
01020 if (p[0]<m_domainBox.smallEnd(0)) p[0]+= m_domainBox.size(0);
01021 else if (p[0]>m_domainBox.bigEnd(0)) p[0]-= m_domainBox.size(0);
01022 },
01023 if (m_isPeriodic[1])
01024 {
01025 if (p[1]<m_domainBox.smallEnd(1)) p[1]+= m_domainBox.size(1);
01026 else if (p[1]>m_domainBox.bigEnd(1)) p[1]-= m_domainBox.size(1);
01027 },
01028 if (m_isPeriodic[2])
01029 {
01030 if (p[2]<m_domainBox.smallEnd(2)) p[2]+= m_domainBox.size(2);
01031 else if (p[2]>m_domainBox.bigEnd(2)) p[2]-= m_domainBox.size(2);
01032 },
01033 if (m_isPeriodic[3])
01034 {
01035 if (p[3]<m_domainBox.smallEnd(3)) p[3]+= m_domainBox.size(3);
01036 else if (p[3]>m_domainBox.bigEnd(3)) p[3]-= m_domainBox.size(3);
01037 },
01038 if (m_isPeriodic[4])
01039 {
01040 if (p[4]<m_domainBox.smallEnd(4)) p[4]+= m_domainBox.size(4);
01041 else if (p[4]>m_domainBox.bigEnd(4)) p[4]-= m_domainBox.size(4);
01042 },
01043 if (m_isPeriodic[5])
01044 {
01045 if (p[5]<m_domainBox.smallEnd(5)) p[5]+= m_domainBox.size(5);
01046 else if (p[5]>m_domainBox.bigEnd(5)) p[5]-= m_domainBox.size(5);
01047 });
01048
01049 return true;
01050 }
01051
01052 inline
01053 bool
01054 ProblemDomain::contains (const Box& b) const
01055 {
01056
01057 if (b.type() == m_domainBox.type())
01058 {
01059 return ( !isEmpty()
01060 && (D_TERM6((m_isPeriodic[0]
01061 || ( b.smallEnd(0) >= m_domainBox.smallEnd(0)
01062 && b.bigEnd (0) <= m_domainBox.bigEnd (0))), &&
01063 (m_isPeriodic[1]
01064 || ( b.smallEnd(1) >= m_domainBox.smallEnd(1)
01065 && b.bigEnd (1) <= m_domainBox.bigEnd (1))), &&
01066 (m_isPeriodic[2]
01067 || ( b.smallEnd(2) >= m_domainBox.smallEnd(2)
01068 && b.bigEnd (2) <= m_domainBox.bigEnd (2))), &&
01069 (m_isPeriodic[3]
01070 || ( b.smallEnd(3) >= m_domainBox.smallEnd(3)
01071 && b.bigEnd (3) <= m_domainBox.bigEnd (3))), &&
01072 (m_isPeriodic[4]
01073 || ( b.smallEnd(4) >= m_domainBox.smallEnd(4)
01074 && b.bigEnd (4) <= m_domainBox.bigEnd (4))), &&
01075 (m_isPeriodic[5]
01076 || ( b.smallEnd(5) >= m_domainBox.smallEnd(5)
01077 && b.bigEnd (5) <= m_domainBox.bigEnd (5))))));
01078 }
01079 else
01080 {
01081 Box domainBox = m_domainBox;
01082
01083
01084 for (int dir = 0; dir < SpaceDim; dir++)
01085 {
01086 if (b.type(dir) != domainBox.type(dir))
01087 {
01088 if (b.type(dir) == IndexType::NODE)
01089 {
01090 domainBox.surroundingNodes(dir);
01091 }
01092 else
01093 {
01094 domainBox.enclosedCells(dir);
01095 }
01096 }
01097 }
01098
01099 return ( !isEmpty()
01100 && (D_TERM6((m_isPeriodic[0]
01101 || ( b.smallEnd(0) >= domainBox.smallEnd(0)
01102 && b.bigEnd (0) <= domainBox.bigEnd (0))), &&
01103 (m_isPeriodic[1]
01104 || ( b.smallEnd(1) >= domainBox.smallEnd(1)
01105 && b.bigEnd (1) <= domainBox.bigEnd (1))), &&
01106 (m_isPeriodic[2]
01107 || ( b.smallEnd(2) >= domainBox.smallEnd(2)
01108 && b.bigEnd (2) <= domainBox.bigEnd (2))), &&
01109 (m_isPeriodic[3]
01110 || ( b.smallEnd(3) >= domainBox.smallEnd(3)
01111 && b.bigEnd (3) <= domainBox.bigEnd (3))), &&
01112 (m_isPeriodic[4]
01113 || ( b.smallEnd(4) >= domainBox.smallEnd(4)
01114 && b.bigEnd (4) <= domainBox.bigEnd (4))), &&
01115 (m_isPeriodic[5]
01116 || ( b.smallEnd(5) >= domainBox.smallEnd(5)
01117 && b.bigEnd (5) <= domainBox.bigEnd (5))))));
01118 }
01119 }
01120
01121 #endif
01122
01123 #include "NamespaceFooter.H"
01124 #endif