Changeset f57182 for src/base


Ignore:
Timestamp:
Mar 30, 2013, 2:44:52 PM (13 years ago)
Author:
Julian Iseringhausen <isering@…>
Children:
8180d8
Parents:
d13e27
git-author:
Julian Iseringhausen <isering@…> (06/11/12 14:02:16)
git-committer:
Julian Iseringhausen <isering@…> (03/30/13 14:44:52)
Message:

Open boundary conditions.

Conflicts:

lib/vmg/src/Makefile.am
lib/vmg/src/base/factory.cpp
lib/vmg/test/unit_test/library/dirichlet_fas_lr_mpi.cpp

Location:
src/base
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/base/factory.cpp

    rd13e27 rf57182  
    7171
    7272  MG::GetComm()->PrintOnce(Debug, "Error: Object %s is not registered", id.c_str());
    73   assert(0);
     73  assert(0 == "Mandatory object not registered.");
    7474
    7575  return NULL;

  • src/base/helper.cpp

    rd13e27 rf57182  
    9090  return (1.0 - coord.Z()) * interpolate_vals[0] + coord.Z() * interpolate_vals[1];
    9191}
     92
     93bool Helper::AssertVectorsEqual(const Vector& pos_1, const Vector& pos_2)
     94{
     95  bool equal = (pos_1 - pos_2).Abs().IsComponentwiseLessOrEqual(std::numeric_limits<vmg_float>::epsilon());
     96  assert(equal);
     97  return equal;
     98}

  • src/base/helper.hpp

    rd13e27 rf57182  
    8080   * Checks a number for validity, i.e. it is neither nan nor inf.
    8181   */
    82   inline bool CheckNumber(const vmg_float& number)
     82  template <class T>
     83  inline bool CheckNumber(const T& number)
    8384  {
    8485    bool valid = true;
    8586
    86     if (std::numeric_limits<vmg_float>::has_quiet_NaN) {
    87       valid &= number != std::numeric_limits<vmg_float>::quiet_NaN();
    88       assert(number != std::numeric_limits<vmg_float>::quiet_NaN());
    89     }
     87    // Check for nan
     88    valid &= number == number;
     89    assert(number == number);
    9090
    91     if (std::numeric_limits<vmg_float>::has_signaling_NaN) {
    92       valid &= number != std::numeric_limits<vmg_float>::signaling_NaN();
    93       assert(number != std::numeric_limits<vmg_float>::signaling_NaN());
    94     }
    95 
    96     if (std::numeric_limits<vmg_float>::has_infinity) {
    97       valid &= number != std::numeric_limits<vmg_float>::infinity();
    98       valid &= number != -1 * std::numeric_limits<vmg_float>::infinity();
    99       assert(number != std::numeric_limits<vmg_float>::infinity());
    100       assert(number != -1 * std::numeric_limits<vmg_float>::infinity());
     91    // Check for infinity
     92    if (!std::numeric_limits<T>::is_integer && std::numeric_limits<T>::has_denorm == std::denorm_present) {
     93      valid &= number >= -1 * std::numeric_limits<T>::max();
     94      valid &= number <= std::numeric_limits<T>::max();
     95      assert(number >= -1 * std::numeric_limits<T>::max());
     96      assert(number <= std::numeric_limits<T>::max());
     97    }else {
     98      valid &= number >= std::numeric_limits<T>::min();
     99      valid &= number <= std::numeric_limits<T>::max();
     100      assert(number >= std::numeric_limits<T>::min());
     101      assert(number <= std::numeric_limits<T>::max());
    101102    }
    102103
    … …  
    179180  }
    180181
     182  bool AssertVectorsEqual(const Vector& pos_1, const Vector& pos_2);
     183
    181184  vmg_float InterpolateTrilinear(const Vector& point, const Grid& grid);
    182185

  • src/base/interface.cpp

    rd13e27 rf57182  
    3838using namespace VMG;
    3939
     40static Index GetGlobalIndex(const Vector& pos, const SpatialExtent& extent, const BT& bt)
     41{
     42  const Index index = (pos - extent.Begin()) / extent.MeshWidth() + 0.5;
     43  return index + (bt == LocallyRefined ? 1 : 0);
     44}
     45
    4046void Interface::InitInterface(const Vector& box_offset, const vmg_float& box_size,
    4147                              const int& coarseningSteps, const vmg_float& alpha)
    4248{
    4349  int i;
    44   Index num_cells = Helper::intpow(2, levelMax);
    4550  Index size_factor;
    4651
    4752  const Vector box_center = box_offset + 0.5 * box_size;
    4853
     54  /*
     55   * Get Extents
     56   */
    4957  for (i=0; i<coarseningSteps; ++i) {
    50 
    51     global.push_back(GlobalIndices());
    5258    extent.push_back(SpatialExtent());
    53 
    5459    for (int j=0; j<3; ++j)
    5560      size_factor[j] = (bc[j] == Periodic ? 1 : Helper::intpow(2, static_cast<int>(log(pow(alpha, i+1)) / log(2.0) + 1.0)));
    … …  
    5863    extent.back().Begin() = box_center - 0.5 * extent.back().Size();
    5964    extent.back().End() = extent.back().Begin() + extent.back().Size();
    60     extent.back().MeshWidth() = pow(2.0, i-levelMax);
    61 
    62     num_cells = Helper::intpow(2,levelMax-i) * size_factor;
    63 
    64     global.back().GlobalSize() = num_cells + 1;
    65     global.back().LocalSize() = num_cells + 1;
    66     global.back().LocalBegin() = 0;
    67     global.back().LocalEnd() = num_cells + 1;
    68 
    69     global.back().FinestAbsSize() = Helper::intpow(2,i) * num_cells + 1;
    70     global.back().FinestAbsBegin() = ((global.back().FinestAbsSize()-1) * (1-size_factor)) / (2*size_factor);
    71     global.back().FinestAbsEnd() = global.back().FinestAbsBegin() + global.back().FinestAbsSize();
    72 
    73     if (i==0)
    74       global.back().GlobalFinerBegin() = (num_cells - Helper::intpow(2, levelMax-i))/2;
    75     else
    76       global.back().GlobalFinerBegin() = (global.back().FinestAbsSize() - (++global.rbegin())->FinestAbsSize()) / Helper::intpow(2,i+1);
    77 
    78     global.back().GlobalFinerEnd() = global.back().GlobalSize() - global.back().GlobalFinerBegin();
    79     global.back().GlobalFinerSize() = global.back().GlobalFinerEnd() - global.back().GlobalFinerBegin();
    80 
    81     global.back().LocalFinerBegin() = global.back().GlobalFinerBegin();
    82     global.back().LocalFinerEnd() = global.back().GlobalFinerEnd();
    83     global.back().LocalFinerSize() = global.back().GlobalFinerSize();
     65    extent.back().MeshWidth() = pow(2.0, i-levelMax) * extent.back().Size() / size_factor;
    8466  }
    8567
    86   while (global.size() == 0 || global.back().GlobalSize().Min() > Helper::intpow(2, levelMin)+1) {
     68  if (extent.size() == 0) {
     69    extent.push_back(SpatialExtent());
     70    extent.back().Size() = box_size;
     71    extent.back().Begin() = box_offset;
     72    extent.back().End() = box_offset + box_size;
     73    extent.back().MeshWidth() = box_size / Helper::intpow(2,levelMax);
     74  }
    8775
    88     if (global.size() > 0)
    89       num_cells /= 2;
     76  while ((extent.back().Size() / extent.back().MeshWidth()).Min() > Helper::intpow(2,levelMin) + 1) {
     77    extent.push_back(SpatialExtent(extent.back()));
     78    extent.back().MeshWidth() *= 2;
     79  }
    9080
     81  /*
     82   * Find GlobalMax
     83   */
     84  for (std::vector<SpatialExtent>::const_iterator iter=extent.begin(); iter!=extent.end(); ++iter) {
    9185    global.push_back(GlobalIndices());
    92     extent.push_back(SpatialExtent());
     86    global.back().BoundaryType() = GlobalCoarsened;
     87  }
    9388
    94     if (global.size() == 1) {
    95       extent.back().Size() = box_size;
    96       extent.back().Begin() = box_offset;
    97       extent.back().End() = box_offset + box_size;
    98       extent.back().MeshWidth() = box_size / static_cast<Vector>(num_cells);
    99     }else {
    100       extent.back().Size() = (++extent.rbegin())->Size();
    101       extent.back().Begin() = (++extent.rbegin())->Begin();
    102       extent.back().End() = (++extent.rbegin())->End();
    103       extent.back().MeshWidth() = 2.0 * (++extent.rbegin())->MeshWidth();
     89  for (i=global.size()-1; i>=0; --i) {
     90    if (i == 0 || extent[i-1].Size() != extent[i].Size()) {
     91      global[i].BoundaryType() = GlobalMax;
     92      break;
     93    }
     94  }
     95  for (--i; i>=0; --i)
     96    global[i].BoundaryType() = LocallyRefined;
     97
     98  /*
     99   * Compute global grid values
     100   */
     101  for (i=0; i<global.size(); ++i) {
     102
     103    const Index num_cells = extent[i].Size() / extent[i].MeshWidth() + 0.5;
     104
     105    for (int j=0; j<3; ++j) {
     106
     107      if (bc[j] == Dirichlet || bc[j] == Open) {
     108        /*
     109         * Dirichlet
     110         */
     111
     112        if (global[i].BoundaryType() == LocallyRefined) {
     113          /*
     114           * Locally Refined
     115           */
     116
     117          global[i].GlobalSize()[j] = num_cells[j] + 3;
     118
     119          if (i == 0) {
     120            global[i].GlobalFinerBegin()[j] = GetGlobalIndex(box_offset, extent[i], global[i].BoundaryType())[j];
     121            global[i].GlobalFinerEnd()[j] = GetGlobalIndex(box_offset + box_size, extent[i], global[i].BoundaryType())[j] + 1;
     122            global[i].GlobalFinerSize()[j] = global[i].GlobalFinerEnd()[j] - global[i].GlobalFinerBegin()[j];
     123          } else {
     124            global[i].GlobalFinerBegin()[j] = GetGlobalIndex(extent[i-1].Begin(), extent[i], global[i].BoundaryType())[j];
     125            global[i].GlobalFinerEnd()[j] = GetGlobalIndex(extent[i-1].End(), extent[i], global[i].BoundaryType())[j] + 1;
     126            global[i].GlobalFinerSize()[j] = global[i].GlobalFinerEnd()[j] - global[i].GlobalFinerBegin()[j];
     127          }
     128
     129        } else if (global[i].BoundaryType() == GlobalMax) {
     130          /*
     131           * Global Max
     132           */
     133
     134          global[i].GlobalSize()[j] = num_cells[j] + 1;
     135
     136          if (i == 0) {
     137            global[i].GlobalFinerBegin()[j] = 0;
     138            global[i].GlobalFinerEnd()[j] = num_cells[j] + 1;
     139            global[i].GlobalFinerSize()[j] = num_cells[j] + 1;
     140          } else {
     141            global[i].GlobalFinerBegin()[j] = GetGlobalIndex(extent[i-1].Begin(), extent[i], global[i].BoundaryType())[j];
     142            global[i].GlobalFinerEnd()[j] = GetGlobalIndex(extent[i-1].End(), extent[i], global[i].BoundaryType())[j] + 1;
     143            global[i].GlobalFinerSize()[j] = global[i].GlobalFinerEnd()[j] - global[i].GlobalFinerBegin()[j];
     144          }
     145
     146        } else {
     147          /*
     148           * Global Coarsened
     149           */
     150
     151          global[i].GlobalSize()[j] = num_cells[j] + 1;
     152
     153          global[i].GlobalFinerBegin()[j] = 0;
     154          global[i].GlobalFinerEnd()[j] = num_cells[j] + 1;
     155          global[i].GlobalFinerSize()[j] = num_cells[j] + 1;
     156
     157        }
     158
     159      } else if (bc[j] == Periodic) {
     160        /*
     161         * Periodic
     162         */
     163
     164        global[i].GlobalSize()[j] = num_cells[j];
     165
     166        global[i].GlobalFinerBegin()[j] = 0;
     167        global[i].GlobalFinerEnd()[j] = num_cells[j];
     168        global[i].GlobalFinerSize()[j] = num_cells[j];
     169
     170      }
    104171    }
    105172
    106     global.back().GlobalSize() = num_cells + 1;
    107     global.back().LocalSize() = num_cells + 1;
    108     global.back().LocalBegin() = 0;
    109     global.back().LocalEnd() = num_cells + 1;
     173    global[i].LocalBegin() = 0;
     174    global[i].LocalEnd() = global[i].GlobalSize();
     175    global[i].LocalSize() = global[i].GlobalSize();
    110176
    111     if (global.size() == 1) {
    112       global.back().FinestAbsBegin() = 0;
    113       global.back().FinestAbsEnd() = global.back().GlobalSize();
    114       global.back().FinestAbsSize() = global.back().GlobalSize();
    115     }else {
    116       global.back().FinestAbsBegin() = (++global.rbegin())->FinestAbsBegin();
    117       global.back().FinestAbsEnd() = (++global.rbegin())->FinestAbsEnd();
    118       global.back().FinestAbsSize() = (++global.rbegin())->FinestAbsSize();
    119     }
    120 
    121     global.back().GlobalFinerBegin() = 0;
    122     global.back().GlobalFinerEnd() = global.back().GlobalSize();
    123     global.back().GlobalFinerSize() = global.back().GlobalSize();
    124 
    125     global.back().LocalFinerBegin() = global.back().GlobalFinerBegin();
    126     global.back().LocalFinerEnd() = global.back().GlobalFinerEnd();
    127     global.back().LocalFinerSize() = global.back().GlobalFinerSize();
     177    global[i].LocalFinerBegin() = global[i].GlobalFinerBegin();
     178    global[i].LocalFinerEnd() = global[i].GlobalFinerEnd();
     179    global[i].LocalFinerSize() = global[i].GlobalFinerSize();
    128180
    129181  }
    130182
    131   for (i=0; i<3; ++i)
    132     if (bc[i] == Periodic)
    133       for (unsigned int j=0; j<global.size(); ++j) {
    134         global[j].GlobalSize()[i] -= 1;
    135         global[j].FinestAbsSize()[i] -= Helper::intpow(2, j);
    136         global[j].FinestAbsEnd()[i] -= Helper::intpow(2, j);
    137         global[j].LocalSize()[i] -= 1;
    138         global[j].LocalEnd()[i] -= 1;
    139         global[j].GlobalFinerSize()[i] -= 1;
    140         global[j].GlobalFinerEnd()[i] -= 1;
    141         global[j].LocalFinerSize()[i] -= 1;
    142         global[j].LocalFinerEnd()[i] -= 1;
    143       }
    144 
    145183  levelMin = levelMax - global.size() + 1;
    146184
    147   global.back().BoundaryType() = GlobalCoarsened;
    148 
    149   for (i=global.size()-2; i>=0; --i) {
    150     if (global[i].FinestAbsSize().Product() >= global[i+1].FinestAbsSize().Product()) {
    151       global[i].BoundaryType() = GlobalCoarsened;
    152     }else {
    153       global[i].BoundaryType() = LocallyRefined;
    154       global[i+1].BoundaryType() = GlobalMax;
    155       break;
    156     }
    157   }
    158 
    159   for (; i>=0; --i)
    160     global[i].BoundaryType() = LocallyRefined;
    161 
    162   if (global.front().BoundaryType() != LocallyRefined &&
    163       global.front().BoundaryType() != GlobalMax)
    164     global.front().BoundaryType() = GlobalMax;
    165 
    166185}
Note: See TracChangeset for help on using the changeset viewer.