Changeset f57182 for src


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
Files:
22 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}

  • src/comm/comm.cpp

    rd13e27 rf57182  
    4747vmg_float Comm::ComputeResidualNorm(Multigrid& sol_mg, Multigrid& rhs_mg)
    4848{
    49 #ifdef DEBUG_MATRIX_CHECKS
    50   sol().IsCompatible(rhs());
    51   sol().IsConsistent();
    52   rhs().IsConsistent();
    53 #endif
    54 
    5549  Stencil::iterator stencil_iter;
    5650  vmg_float norm = 0.0;
    … …  
    6660  const Grid& sol = sol_mg();
    6761  const Grid& rhs = rhs_mg();
     62
     63#ifdef DEBUG_MATRIX_CHECKS
     64  sol.IsCompatible(rhs);
     65  sol.IsConsistent();
     66  rhs.IsConsistent();
     67#endif
    6868
    6969  for (int i=rhs.Local().Begin().X(); i<rhs.Local().End().X(); ++i)

  • src/comm/comm_mpi.cpp

    rd13e27 rf57182  
    487487          << "    GLOBAL_FINER_END:   " << grid.Global().GlobalFinerEnd() << std::endl
    488488          << "    GLOBAL_FINER_SIZE:  " << grid.Global().GlobalFinerSize() << std::endl
    489           << "    FINEST_ABS_BEGIN:   " << grid.Global().FinestAbsBegin() << std::endl
    490           << "    FINEST_ABS_END:     " << grid.Global().FinestAbsEnd() << std::endl
    491           << "    FINEST_ABS_SIZE:    " << grid.Global().FinestAbsSize() << std::endl
    492489          << "    GLOBAL_SIZE:        " << grid.Global().GlobalSize() << std::endl
    493490          << "  EXTENT:" << std::endl
    … …  
    577574  if (settings.CommunicatorLocal(grid) != MPI_COMM_NULL) {
    578575
    579     Index end, end_global;
    580 
    581     for (int i=0; i<3; ++i) {
    582       end[i] = grid.Local().End()[i];
    583       end_global[i] = grid.Global().LocalEnd()[i];
    584     }
    585 
    586576    vtkSmartPointer<vtkImageData> image = vtkSmartPointer<vtkImageData>::New();
    587     image->SetExtent(grid.Global().LocalBegin().X(), end_global.X()-1,
    588                      grid.Global().LocalBegin().Y(), end_global.Y()-1,
    589                      grid.Global().LocalBegin().Z(), end_global.Z()-1);
     577    image->SetExtent(grid.Global().LocalBegin().X(), grid.Global().LocalEnd().X()-1,
     578                     grid.Global().LocalBegin().Y(), grid.Global().LocalEnd().Y()-1,
     579                     grid.Global().LocalBegin().Z(), grid.Global().LocalEnd().Z()-1);
    590580    image->SetSpacing(grid.Extent().MeshWidth().vec());
    591581    image->SetOrigin(grid.Extent().Begin().vec());
    … …  
    595585    image->GetPointData()->GetScalars()->SetName(information);
    596586
    597     Index i;
    598     for (i.X()=grid.Local().Begin().X(); i.X()<end.X(); ++i.X())
    599       for (i.Y()=grid.Local().Begin().Y(); i.Y()<end.Y(); ++i.Y())
    600         for (i.Z()=grid.Local().Begin().Z(); i.Z()<end.Z(); ++i.Z())
    601           image->SetScalarComponentFromDouble(i.X() - grid.Local().Begin().X() + grid.Global().LocalBegin().X(),
    602                                               i.Y() - grid.Local().Begin().Y() + grid.Global().LocalBegin().Y(),
    603                                               i.Z() - grid.Local().Begin().Z() + grid.Global().LocalBegin().Z(),
     587    Index begin, end, i;
     588    for (int j=0; j<3; ++j) {
     589      begin[j] = (grid.Local().BoundarySize1()[j] > 0 ? grid.Local().BoundaryBegin1()[j] : grid.Local().Begin()[j]);
     590      end[j] = (grid.Local().BoundarySize2()[j] > 0 ? grid.Local().BoundaryEnd2()[j] : grid.Local().End()[j]);
     591    }
     592
     593    for (i.X() = begin.X(); i.X() < end.X(); ++i.X())
     594      for (i.Y() = begin.Y(); i.Y() < end.Y(); ++i.Y())
     595        for (i.Z() = begin.Z(); i.Z() < end.Z(); ++i.Z())
     596          image->SetScalarComponentFromDouble(i.X() - begin.X() + grid.Global().LocalBegin().X(),
     597                                              i.Y() - begin.Y() + grid.Global().LocalBegin().Y(),
     598                                              i.Z() - begin.Z() + grid.Global().LocalBegin().Z(),
    604599                                              0, grid.GetVal(i));
    605600

  • src/comm/comm_serial.cpp

    rd13e27 rf57182  
    333333        << "  LOCAL_FINER_END:    " << (*mg)(i).Global().LocalFinerEnd() << std::endl
    334334        << "  LOCAL_FINER_SIZE:   " << (*mg)(i).Global().LocalFinerSize() << std::endl
    335         << "  FINEST_ABS_BEGIN:   " << (*mg)(i).Global().FinestAbsBegin() << std::endl
    336         << "  FINEST_ABS_END:     " << (*mg)(i).Global().FinestAbsEnd() << std::endl
    337         << "  FINEST_ABS_SIZE:    " << (*mg)(i).Global().FinestAbsSize() << std::endl
    338335        << "  GLOBAL_SIZE:        " << (*mg)(i).Global().GlobalSize() << std::endl
    339336        << "  BOUNDARY_TYPE:      " << (*mg)(i).Global().BoundaryType() << std::endl

  • src/comm/domain_decomposition_mpi.cpp

    rd13e27 rf57182  
    5555    global_l.GlobalFinerEnd() = interface->Global()[i].GlobalFinerEnd();
    5656    global_l.GlobalFinerSize() = interface->Global()[i].GlobalFinerSize();
    57     global_l.FinestAbsBegin() = interface->Global()[i].FinestAbsBegin();
    58     global_l.FinestAbsEnd() = interface->Global()[i].FinestAbsEnd();
    59     global_l.FinestAbsSize() = interface->Global()[i].FinestAbsSize();
    6057    global_l.GlobalSize() = interface->Global()[i].GlobalSize();
    6158    global_l.BoundaryType() = interface->Global()[i].BoundaryType();
    … …  
    8077        global_l.LocalEnd() = global_l.LocalBegin() + global_l.LocalSize();
    8178
    82         global_l.LocalFinerBegin() = 0;
    83         global_l.LocalFinerEnd() = 0;
    84         global_l.LocalFinerSize() = 0;
     79        global_l.LocalFinerBegin() = global_l.LocalBegin().Clamp(global_l.GlobalFinerBegin(), global_l.GlobalFinerEnd());
     80        global_l.LocalFinerEnd() = global_l.LocalEnd().Clamp(global_l.GlobalFinerBegin(), global_l.GlobalFinerEnd());
     81        global_l.LocalFinerSize() = global_l.LocalFinerEnd() - global_l.LocalFinerBegin();
    8582
    8683      }else {

  • src/commands/com_check_relative_residual.cpp

    rd13e27 rf57182  
    5454    const vmg_float& init_res = factory.GetObjectStorageVal<vmg_float>(arguments[1]);
    5555    const vmg_float& precision = factory.GetObjectStorageVal<vmg_float>("PRECISION");
    56     const vmg_float rel_res = std::fabs(res / init_res);
     56    const vmg_float rel_res = std::abs(res / init_res);
    5757
    5858    MG::GetComm()->PrintOnce(Info, "Relative residual: %e", rel_res);

  • src/grid/grid.cpp

    rd13e27 rf57182  
    270270  return consistent;
    271271}
     272
     273Vector Grid::GetSpatialPos(const Index& index_local) const
     274{
     275  if (Global().BoundaryType() == LocallyRefined)
     276    return Extent().Begin() + Extent().MeshWidth() * (index_local + Global().LocalBegin() - 1);
     277  else
     278    return Extent().Begin() + Extent().MeshWidth() * (index_local - Local().HaloSize1() + Global().LocalBegin());
     279}
     280
     281
     282Index Grid::GetGlobalIndex(const Vector& pos) const
     283{
     284  const Index index = (pos - Extent().Begin()) / Extent().MeshWidth() + 0.5;
     285  return index + (Global().BoundaryType() == LocallyRefined ? 1 : 0);
     286}
     287

  • src/grid/grid.hpp

    rd13e27 rf57182  
    107107  vmg_float& operator()(const Index& index);
    108108
     109  const vmg_float& operator()(int x, int y, int z) const;  ///< Returns a reference to the requested gridpoint.
     110  const vmg_float& operator()(const Index& index) const;
     111
    109112  const vmg_float& GetVal(int x, int y, int z) const; ///< Returns the value of a requested gridpoint.
    110113  const vmg_float& GetVal(const Index& index) const;
    … …  
    138141  bool IsActive() const {return Local().Size().Product() > 0;}
    139142
     143  Vector GetSpatialPos(const Index& index_local) const;
     144  Index GetGlobalIndex(const Vector& pos) const;
     145
    140146private:
    141147  void InitGrid();
    … …  
    167173}
    168174
     175inline const vmg_float& Grid::operator()(int x, int y, int z) const
     176{
     177  return grid[z + local.SizeTotal().Z() * (y + local.SizeTotal().Y() * x)];
     178}
     179
     180inline const vmg_float& Grid::operator()(const Index& index) const
     181{
     182  return this->operator()(index.X(), index.Y(), index.Z());
     183}
     184
    169185inline const vmg_float& Grid::GetVal(int x, int y, int z) const
    170186{

  • src/grid/grid_iterator_suite.cpp

    rd13e27 rf57182  
    3838{
    3939  local.SetBounds(local_.Begin(), local_.End());
     40  local_finer.SetBounds(local_.FinerBegin(), local_.FinerEnd());
    4041  complete_grid.SetBounds(0, local_.SizeTotal());
    4142  inner_local_grid.SetBounds(local_.Begin()+local_.HaloSize1(), local_.End()-local_.HaloSize2());

  • src/grid/grid_iterator_suite.hpp

    rd13e27 rf57182  
    5252  GridIteratorSuite(const GridIteratorSuite& other) :
    5353    local(other.local),
     54    local_finer(other.local_finer),
    5455    complete_grid(other.complete_grid),
    5556    inner_local_grid(other.inner_local_grid),
    … …  
    6869
    6970  const GridIteratorSet& Local() const {return local;}
     71  const GridIteratorSet& LocalFiner() const {return local_finer;}
    7072  const GridIteratorSet& CompleteGrid() const {return complete_grid;}
    7173  const GridIteratorSet& InnerLocalGrid() const {return inner_local_grid;}
    … …  
    8284private:
    8385  GridIteratorSet local;
     86  GridIteratorSet local_finer;
    8487  GridIteratorSet complete_grid;
    8588  GridIteratorSet inner_local_grid;

  • src/grid/grid_properties.hpp

    rd13e27 rf57182  
    4343    global_finer_begin(0), global_finer_end(0), global_finer_size(0),
    4444    local_finer_begin(0), local_finer_end(0), local_finer_size(0),
    45     finest_abs_begin(0), finest_abs_end(0), finest_abs_size(0),
    4645    global_size(0),
    4746    boundary(EmptyGrid)
    … …  
    5251    global_finer_begin(other.global_finer_begin), global_finer_end(other.global_finer_end), global_finer_size(other.global_finer_size),
    5352    local_finer_begin(other.local_finer_begin), local_finer_end(other.local_finer_end), local_finer_size(other.local_finer_size),
    54     finest_abs_begin(other.finest_abs_begin), finest_abs_end(other.finest_abs_end), finest_abs_size(other.finest_abs_size),
    5553    global_size(other.global_size),
    5654    boundary(other.boundary)
    … …  
    8078  const Index& LocalFinerEnd() const {return local_finer_end;}
    8179  const Index& LocalFinerSize() const {return local_finer_size;}
    82 
    83   Index& FinestAbsBegin() {return finest_abs_begin;}
    84   Index& FinestAbsEnd() {return finest_abs_end;}
    85   Index& FinestAbsSize() {return finest_abs_size;}
    86 
    87   const Index& FinestAbsBegin() const {return finest_abs_begin;}
    88   const Index& FinestAbsEnd() const {return finest_abs_end;}
    89   const Index& FinestAbsSize() const {return finest_abs_size;}
    9080
    9181  Index& GlobalSize() {return global_size;}
    … …  
    9989  Index global_finer_begin, global_finer_end, global_finer_size;
    10090  Index local_finer_begin, local_finer_end, local_finer_size;
    101   Index finest_abs_begin, finest_abs_end, finest_abs_size;
    10291  Index global_size;
    10392  BT boundary;

  • src/grid/multigrid.cpp

    rd13e27 rf57182  
    8585            comm->BoundaryConditions()[j] == Open) {
    8686
    87           local_l.SizeTotal()[j] = global_separated[i].LocalSize()[j] +
    88             (global_separated[i].BoundaryType() == LocallyRefined ? 2 : 0);
     87          local_l.SizeTotal()[j] = global_separated[i].LocalSize()[j];
    8988
    9089          /*
    … …  
    141140      local_l.BoundarySize2() = local_l.BoundaryEnd2() - local_l.BoundaryBegin2();
    142141
    143       local_l.FinerSize() = global_separated[i].LocalFinerSize() - local_l.BoundarySize1() - local_l.BoundarySize2();
    144       local_l.FinerBegin() = global_separated[i].LocalFinerBegin() - global_separated[i].LocalBegin() + local_l.Begin();
     142      local_l.FinerSize() = global_separated[i].LocalFinerSize();
     143      local_l.FinerBegin() = global_separated[i].LocalFinerBegin() - global_separated[i].LocalBegin() + local_l.HaloSize1();
    145144      local_l.FinerEnd() = local_l.FinerBegin() + local_l.FinerSize();
    146145

  • src/grid/tempgrid.cpp

    rd13e27 rf57182  
    7171  global.LocalFinerEnd() = 0;
    7272  global.LocalFinerSize() = 0;
    73   global.FinestAbsBegin() = 0;
    74   global.FinestAbsEnd() = size;
    75   global.FinestAbsSize() = size;
    7673  global.GlobalSize() = size;
    7774  global.BoundaryType() = BTUndefined;
    … …  
    113110  const Index off = GridIndexTranslations::EndOffset(boundary);
    114111
     112  /*
     113   * Set global grid attributes
     114   */
     115
    115116  level = grid.Level() + 1;
     117
     118  global.BoundaryType() = grid_finer.Global().BoundaryType();
    116119
    117120  global.GlobalFinerBegin() = 0;
    … …  
    123126  global.LocalFinerSize() = 0;
    124127
    125   global.FinestAbsBegin() = grid_finer.Global().FinestAbsBegin();
    126   global.FinestAbsEnd() = grid_finer.Global().FinestAbsEnd();
    127   global.FinestAbsSize() = grid_finer.Global().FinestAbsSize();
     128  global.LocalBegin() = 2 * (grid.Global().LocalFinerBegin() - grid.Global().GlobalFinerBegin());
     129  global.LocalEnd() = 2 * (grid.Global().LocalFinerEnd() - grid.Global().GlobalFinerBegin() - off) + off;
     130  global.LocalSize() = global.LocalEnd() - global.LocalBegin();
    128131
    129132  global.GlobalSize() = 2 * (grid.Global().GlobalFinerSize() - off) + off;
    130   global.BoundaryType() = grid_finer.Global().BoundaryType();
    131 
    132   global.LocalBegin() = 2 * grid.Global().LocalBegin().Clamp(grid.Global().GlobalFinerBegin(), grid.Global().GlobalFinerEnd());
    133   global.LocalEnd() = 2 * grid.Global().LocalEnd().Clamp(grid.Global().GlobalFinerBegin(), grid.Global().GlobalFinerEnd()) - off;
     133
     134  for (int j=0; j<3; ++j) {
     135
     136    if (boundary[j] == Dirichlet ||
     137        boundary[j] == Open) {
     138
     139      if (grid.Global().BoundaryType() == LocallyRefined) {
     140
     141        global.GlobalSize()[j] += 2;
     142        global.LocalBegin()[j] += 1;
     143        global.LocalEnd()[j] +=1;
     144
     145      } else if (grid.Global().BoundaryType() == GlobalMax) {
     146
     147        global.GlobalSize()[j] += 2;
     148        global.LocalBegin()[j] += 1;
     149        global.LocalEnd()[j] +=1;
     150
     151      }
     152
     153    } else {
     154
     155    }
     156
     157  }
    134158
    135159  global.LocalSize() = global.LocalEnd() - global.LocalBegin();
    … …  
    168192    }
    169193
    170     if (grid.Local().BoundarySize1()[i] > 0) {
     194    if (grid.Local().BoundarySize1()[i] > 0 && global.BoundaryType() != LocallyRefined) {
     195      local.Size()[i] -= grid.Local().BoundarySize1()[i];
    171196      local.Begin()[i] += grid.Local().BoundarySize1()[i];
    172197      local.BoundaryEnd1()[i] = grid.Local().BoundarySize1()[i];
    … …  
    178203    }
    179204
    180     if (grid.Local().BoundarySize2()[i] > 0) {
     205    if (grid.Local().BoundarySize2()[i] > 0 && global.BoundaryType() != LocallyRefined) {
     206      local.Size()[i] -= grid.Local().BoundarySize2()[i];
    181207      local.End()[i] -= grid.Local().BoundarySize2()[i];
    182208      local.BoundaryBegin2()[i] = local.End()[i];
    183       local.BoundaryEnd2()[i] = local.End()[i]+grid.Local().BoundarySize2()[i];
     209      local.BoundaryEnd2()[i] = local.End()[i] + grid.Local().BoundarySize2()[i];
    184210    }
    185211
    … …  
    195221    local.BoundarySize1() + local.BoundarySize2();
    196222
    197   extent.Size() = grid.Extent().Size();
    198   extent.Begin() = grid.Extent().Begin();
    199   extent.End() = grid.Extent().End();
    200   extent.MeshWidth() = 0.5 * grid.Extent().MeshWidth();
     223  extent = grid_finer.Extent();
    201224
    202225  iterators.SetSubgrids(local);
    … …  
    218241  global.GlobalFinerEnd() = grid_coarser.Global().GlobalFinerEnd();
    219242  global.GlobalFinerSize() = grid_coarser.Global().GlobalFinerSize();
    220 
    221   global.FinestAbsBegin() = grid.Global().FinestAbsBegin();
    222   global.FinestAbsEnd() = grid.Global().FinestAbsEnd();
    223   global.FinestAbsSize() = grid.Global().FinestAbsSize();
    224243
    225244  global.GlobalSize() = (grid.Global().GlobalSize() - off)/2 + off;
    … …  
    326345  for (iter=Iterators().Local().Begin(); iter!=Iterators().Local().End(); ++iter)
    327346    (*this)(*iter) = rhs.GetVal(*iter) - prefactor * A.Apply(sol, *iter);
     347
     348  this->ClearBoundary();
    328349}
    329350

  • src/level/level_operator_cs.cpp

    rd13e27 rf57182  
    7575  const GridIteratorSet bounds_c(begin_c, end_c);
    7676
     77  assert(bounds_c.Begin().GetEnd()-bounds_c.Begin().GetBegin() == ((bounds_f.Begin().GetEnd()-bounds_f.Begin().GetBegin())-1)/2+1);
     78
    7779  // Compute residual
    7880  TempGrid *temp = MG::GetTempGrid();
    … …  
    119121  const GridIteratorSet bounds_c(sol_c.Local().FinerBegin(), sol_c.Local().FinerEnd());
    120122
     123  assert(bounds_c.Begin().GetEnd()-bounds_c.Begin().GetBegin() == ((bounds_f.Begin().GetEnd()-bounds_f.Begin().GetBegin())-1)/2+1);
     124
    121125  sol_f_undist.ClearHalo();
    122126

  • src/level/level_operator_fas.cpp

    rd13e27 rf57182  
    3030#endif
    3131
     32#include <limits>
     33
    3234#include "base/discretization.hpp"
     35#include "base/helper.hpp"
    3336#include "base/index.hpp"
    3437#include "comm/comm.hpp"
    … …  
    5659  Grid& rhs_c_undist = comm.GetCoarserGrid(rhs);
    5760
     61  Index begin_c = rhs_c_undist.Local().FinerBegin();
     62  Index end_c = rhs_c_undist.Local().FinerEnd();
     63
    5864  Index begin_f = rhs_f.Local().Begin();
    5965  Index end_f = rhs_f.Local().End();
    6066
    6167  if (rhs_c_undist.Global().BoundaryType() == GlobalCoarsened) {
     68    begin_c += rhs_c_undist.Local().BoundarySize1();
     69    end_c -= rhs_c_undist.Local().BoundarySize2();
    6270    begin_f += rhs_f.Local().BoundarySize1();
    6371    end_f -= rhs_f.Local().BoundarySize2();
    6472  }
    6573
     74  Helper::AssertVectorsEqual(rhs_c_undist.GetSpatialPos(begin_c), rhs_f.GetSpatialPos(begin_f));
     75
    6676  const GridIteratorSet bounds_f(begin_f, end_f);
    67   const GridIteratorSet bounds_c(rhs_c_undist.Local().FinerBegin(), rhs_c_undist.Local().FinerEnd());
     77  const GridIteratorSet bounds_c(begin_c, end_c);
     78
     79  assert(bounds_c.Begin().GetEnd()-bounds_c.Begin().GetBegin() == ((bounds_f.Begin().GetEnd()-bounds_f.Begin().GetBegin())-1)/2+1);
    6880
    6981  const Stencil& op_res = OperatorRestrict();
    … …  
    7587  comm.CommToGhosts(sol_f);
    7688  comm.CommSubgrid(sol_c_dist, sol_c_undist, 0);
     89  comm.CommSubgrid(rhs_c_dist, rhs_c_undist, 0);
    7790
    7891  MG::GetDiscretization()->SetInnerBoundary(sol_f, rhs_f, sol_c_undist);
    … …  
    8396  res_grid->ImportFromResidual(sol_f, rhs_f);
    8497
    85   for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c)
     98  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
     99    Helper::AssertVectorsEqual(sol_c_undist.GetSpatialPos(*iter_c), sol_f.GetSpatialPos(*iter_f));
    86100    sol_c_undist(*iter_c) = op_sol.Apply(sol_f, *iter_f);
     101  }
    87102
    88   for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c)
    89     rhs_c_undist(*iter_c) = op_res.Apply(*res_grid, *iter_f);
     103  comm.CommToGhosts(sol_c_undist);
     104
     105  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
     106    Helper::AssertVectorsEqual(rhs_c_undist.GetSpatialPos(*iter_c), (*res_grid).GetSpatialPos(*iter_f));
     107    rhs_c_undist(*iter_c) = op_res.Apply(*res_grid, *iter_f) + prefactor * op_pde.Apply(sol_c_undist, *iter_c);
     108  }
     109
     110  /*
     111  if (rhs_c_undist.Global().BoundaryType() == GlobalCoarsened)
     112    rhs_c_undist.SetBoundary(sol_c_undist);
     113  */
    90114
    91115  comm.CommSubgrid(sol_c_undist, sol_c_dist, 1);
    92   comm.CommToGhosts(sol_c_dist);
    93 
    94116  comm.CommSubgrid(rhs_c_undist, rhs_c_dist, 1);
    95117
    … …  
    98120  sol_old->SetGrid(sol_c_dist);
    99121
    100   for (iter_c=rhs_c_dist.Iterators().Local().Begin(); iter_c!=rhs_c_dist.Iterators().Local().End(); ++iter_c)
    101     rhs_c_dist(*iter_c) += prefactor * op_pde.Apply(sol_c_dist, *iter_c);
    102 
    103   comm.CommToGhosts(rhs_c_dist);
     122  //  comm.CommToGhosts(rhs_c_dist);
    104123
    105124  sol.ToCoarserLevel();
    … …  
    125144  TempGrid *sol_old = this->GetTempGrid(sol_c.Level());
    126145
     146  Index begin_c = sol_c.Local().FinerBegin();
     147  Index end_c = sol_c.Local().FinerEnd();
     148
    127149  Index begin_f = sol_f_undist.Local().Begin();
    128150  Index end_f = sol_f_undist.Local().End();
    129151
    130152  if (sol_c.Global().BoundaryType() == GlobalCoarsened) {
    131     begin_f += rhs_f_undist.Local().BoundarySize1();
    132     end_f -= rhs_f_undist.Local().BoundarySize2();
     153    begin_c += sol_c.Local().BoundarySize1();
     154    end_c -= sol_c.Local().BoundarySize2();
     155    begin_f += sol_f_undist.Local().BoundarySize1();
     156    end_f -= sol_f_undist.Local().BoundarySize2();
    133157  }
    134158
     159  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(begin_c), sol_f_undist.GetSpatialPos(begin_f));
     160
    135161  const GridIteratorSet bounds_f(begin_f, end_f);
    136   const GridIteratorSet bounds_c(sol_c.Local().FinerBegin(), sol_c.Local().FinerEnd());
     162  const GridIteratorSet bounds_c(begin_c, end_c);
    137163
     164  assert(bounds_c.Begin().GetEnd()-bounds_c.Begin().GetBegin() == ((bounds_f.Begin().GetEnd()-bounds_f.Begin().GetBegin())-1)/2+1);
     165
     166
     167  comm.CommSubgrid(sol_f_dist, sol_f_undist, 0);
    138168  sol_f_undist.ClearHalo();
    139169
    140170  for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_f!=bounds_f.End(); iter_f+=2, ++iter_c) {
     171    Helper::AssertVectorsEqual(sol_c.GetSpatialPos(*iter_c), sol_f_undist.GetSpatialPos(*iter_f));
    141172    val = sol_c.GetVal(*iter_c) - sol_old->GetVal(*iter_c);
    142173    sol_f_undist(*iter_f) += op.GetDiag() * val;

  • src/samples/discretization_poisson_fd.cpp

    rd13e27 rf57182  
    1818
    1919/**
    20  * @file   discretization_poisson_fd_collatz.cpp
     20 * @file   discretization_poisson_fd.cpp
    2121 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
    2222 * @date   Mon Apr 18 13:03:47 2011
    2323 *
    2424 * @brief  Discretization of the poisson equation
    25  *         using the Collatz Mehrstellen Ansatz.
    26  *         Discretization error: O(h^4)
     25 *         using the standard 7-point stencil or
     26 *         the Collatz Mehrstellen Ansatz.
     27 *         Discretization error: O(h^2) or O(h^4)
    2728 *
    2829 */
    … …  
    3839using namespace VMG;
    3940
    40 DiscretizationPoissonFD::DiscretizationPoissonFD(const int& order) :
    41   Discretization(order)
     41void DiscretizationPoissonFD::InitDiscretizationPoissonFD()
    4242{
    4343  switch (order)
    … …  
    8383  if (order == 4) {
    8484
    85     Grid& rhs = MG::GetRhsMaxLevel();
    86 
    8785    Stencil stencil(6.0/12.0);
    8886    stencil.push_back(-1,  0,  0, 1.0/12.0);
    … …  
    9391    stencil.push_back( 0,  0,  1, 1.0/12.0);
    9492
    95     stencil.Apply(rhs);
     93    stencil.Apply(MG::GetRhsMaxLevel());
    9694
    9795  }

  • src/samples/discretization_poisson_fd.hpp

    rd13e27 rf57182  
    3939{
    4040public:
    41   DiscretizationPoissonFD(const int& order);
     41  DiscretizationPoissonFD(const int& order) :
     42    Discretization(order)
     43  {
     44    InitDiscretizationPoissonFD();
     45  }
    4246
    4347  vmg_float OperatorPrefactor(const Grid& grid) const
    … …  
    4751
    4852  void ModifyRightHandSide();
     53
     54private:
     55  void InitDiscretizationPoissonFD();
    4956};
    5057

  • src/samples/discretization_poisson_fv.cpp

    rd13e27 rf57182  
    3737using namespace VMG;
    3838
     39void DiscretizationPoissonFV::InitDiscretizationPoissonFV()
     40{
     41  switch (order)
     42    {
     43    case 2:
     44      stencil.SetDiag(6.0);
     45      stencil.push_back(-1,  0,  0, -1.0);
     46      stencil.push_back( 1,  0,  0, -1.0);
     47      stencil.push_back( 0, -1,  0, -1.0);
     48      stencil.push_back( 0,  1,  0, -1.0);
     49      stencil.push_back( 0,  0, -1, -1.0);
     50      stencil.push_back( 0,  0,  1, -1.0);
     51      break;
     52    case 4:
     53      stencil.SetDiag(24.0/6.0);
     54      stencil.push_back(-1,  0,  0, -2.0/6.0);
     55      stencil.push_back( 1,  0,  0, -2.0/6.0);
     56      stencil.push_back( 0, -1,  0, -2.0/6.0);
     57      stencil.push_back( 0,  1,  0, -2.0/6.0);
     58      stencil.push_back( 0,  0, -1, -2.0/6.0);
     59      stencil.push_back( 0,  0,  1, -2.0/6.0);
     60      stencil.push_back(-1, -1,  0, -1.0/6.0);
     61      stencil.push_back(-1,  1,  0, -1.0/6.0);
     62      stencil.push_back( 1, -1,  0, -1.0/6.0);
     63      stencil.push_back( 1,  1,  0, -1.0/6.0);
     64      stencil.push_back(-1,  0, -1, -1.0/6.0);
     65      stencil.push_back(-1,  0,  1, -1.0/6.0);
     66      stencil.push_back( 1,  0, -1, -1.0/6.0);
     67      stencil.push_back( 1,  0,  1, -1.0/6.0);
     68      stencil.push_back( 0, -1, -1, -1.0/6.0);
     69      stencil.push_back( 0, -1,  1, -1.0/6.0);
     70      stencil.push_back( 0,  1, -1, -1.0/6.0);
     71      stencil.push_back( 0,  1,  1, -1.0/6.0);
     72      break;
     73    default:
     74      assert(0 != "vmg choose discretization order 2 or 4");
     75      break;
     76    }
     77}
     78
     79void DiscretizationPoissonFV::ModifyRightHandSide()
     80{
     81  if (order == 4) {
     82
     83    Grid& rhs = MG::GetRhsMaxLevel();
     84
     85    Stencil stencil(6.0/12.0);
     86    stencil.push_back(-1,  0,  0, 1.0/12.0);
     87    stencil.push_back( 1,  0,  0, 1.0/12.0);
     88    stencil.push_back( 0, -1,  0, 1.0/12.0);
     89    stencil.push_back( 0,  1,  0, 1.0/12.0);
     90    stencil.push_back( 0,  0, -1, 1.0/12.0);
     91    stencil.push_back( 0,  0,  1, 1.0/12.0);
     92
     93    stencil.Apply(rhs);
     94
     95  }
     96}
     97
    3998void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
    4099{
    … …  
    48107  const Index b2_f = rhs_f.Local().SizeTotal() - 1;
    49108
    50   const Index begin_f = sol_f.Local().Begin();
    51   const Index end_f = sol_f.Local().End();
    52   const Index begin_c = sol_c.Local().FinerBegin();
     109  const Index& begin_f = sol_f.Local().Begin();
     110  const Index& end_f = sol_f.Local().End();
     111  const Index& begin_c = sol_c.Local().FinerBegin();
    53112
    54113  const vmg_float c_1_3 = 1.0 / 3.0;
    … …  
    56115  const vmg_float c_4_3 = 4.0 / 3.0;
    57116
     117  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerBegin()), sol_f.GetSpatialPos(sol_f.Local().Begin()));
     118  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerEnd()-1), sol_f.GetSpatialPos(sol_f.Local().End()-1));
     119
    58120  //
    59121  // X-direction
    60122  //
     123  i_f = begin_f; i_c = begin_c;
    61124  for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
    62125    for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z()) {
     126      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
    63127      rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b1_c.X()-1,i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z())) * h2_inv.X();
    64128      rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b2_c.X()+1,i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z())) * h2_inv.X();
    … …  
    106170  // Y-direction
    107171  //
    108 
     172  i_f = begin_f; i_c = begin_c;
    109173  for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
    110174    for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z()) {
     175      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
    111176      rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b1_c.Y()-1,i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z())) * h2_inv.Y();
    112177      rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b2_c.Y()+1,i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z())) * h2_inv.Y();
    … …  
    154219  // Z-direction
    155220  //
    156 
     221  i_f = begin_f; i_c = begin_c;
    157222  for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
    158223    for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y()) {
     224      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
    159225      rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b1_c.Z()-1) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1)) * h2_inv.Z();
    160226      rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b2_c.Z()+1) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1)) * h2_inv.Z();

  • src/samples/discretization_poisson_fv.hpp

    rd13e27 rf57182  
    4444    Discretization(order)
    4545  {
    46     stencil.SetDiag(6.0);
    47     stencil.push_back(-1,  0,  0, -1.0);
    48     stencil.push_back( 1,  0,  0, -1.0);
    49     stencil.push_back( 0, -1,  0, -1.0);
    50     stencil.push_back( 0,  1,  0, -1.0);
    51     stencil.push_back( 0,  0, -1, -1.0);
    52     stencil.push_back( 0,  0,  1, -1.0);
     46    InitDiscretizationPoissonFV();
    5347  }
    5448
    … …  
    5852  }
    5953
     54  void ModifyRightHandSide();
     55
    6056private:
     57  void InitDiscretizationPoissonFV();
    6158  void SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const;
    6259};
Note: See TracChangeset for help on using the changeset viewer.