| [de061d] | 1 | /*
 | 
|---|
 | 2 |  *    vmg - a versatile multigrid solver
 | 
|---|
 | 3 |  *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
 | 
|---|
 | 4 |  *
 | 
|---|
 | 5 |  *  vmg is free software: you can redistribute it and/or modify
 | 
|---|
 | 6 |  *  it under the terms of the GNU General Public License as published by
 | 
|---|
 | 7 |  *  the Free Software Foundation, either version 3 of the License, or
 | 
|---|
 | 8 |  *  (at your option) any later version.
 | 
|---|
 | 9 |  *
 | 
|---|
 | 10 |  *  vmg is distributed in the hope that it will be useful,
 | 
|---|
 | 11 |  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 | 
|---|
 | 12 |  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 | 
|---|
 | 13 |  *  GNU General Public License for more details.
 | 
|---|
 | 14 |  *
 | 
|---|
 | 15 |  *  You should have received a copy of the GNU General Public License
 | 
|---|
 | 16 |  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 | 
|---|
 | 17 |  */
 | 
|---|
 | 18 | 
 | 
|---|
 | 19 | /**
 | 
|---|
 | 20 |  * @file   tempgrid.cpp
 | 
|---|
 | 21 |  * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 | 
|---|
 | 22 |  * @date   Mon Apr 18 12:55:05 2011
 | 
|---|
 | 23 |  *
 | 
|---|
 | 24 |  * @brief  VMG::TempGrid
 | 
|---|
 | 25 |  *
 | 
|---|
 | 26 |  */
 | 
|---|
 | 27 | 
 | 
|---|
 | 28 | #ifdef HAVE_CONFIG_H
 | 
|---|
 | 29 | #include <libvmg_config.h>
 | 
|---|
 | 30 | #endif
 | 
|---|
 | 31 | 
 | 
|---|
 | 32 | #include "base/discretization.hpp"
 | 
|---|
 | 33 | #include "base/interface.hpp"
 | 
|---|
 | 34 | #include "base/stencil.hpp"
 | 
|---|
 | 35 | #include "comm/comm.hpp"
 | 
|---|
 | 36 | #include "grid/grid_index_translations.hpp"
 | 
|---|
 | 37 | #include "grid/tempgrid.hpp"
 | 
|---|
 | 38 | #include "mg.hpp"
 | 
|---|
 | 39 | 
 | 
|---|
 | 40 | using namespace VMG;
 | 
|---|
 | 41 | 
 | 
|---|
 | 42 | void TempGrid::SetProperties(const Grid& rhs)
 | 
|---|
 | 43 | {
 | 
|---|
 | 44 |   local = rhs.Local();
 | 
|---|
 | 45 |   global = rhs.Global();
 | 
|---|
 | 46 |   extent = rhs.Extent();
 | 
|---|
 | 47 |   iterators.SetSubgrids(rhs.Local());
 | 
|---|
 | 48 |   level = rhs.Level();
 | 
|---|
 | 49 |   Allocate();
 | 
|---|
 | 50 | }
 | 
|---|
 | 51 | 
 | 
|---|
 | 52 | void TempGrid::SetProperties(const GlobalIndices& global_, const LocalIndices& local_, const SpatialExtent& extent_)
 | 
|---|
 | 53 | {
 | 
|---|
 | 54 |   local = local_;
 | 
|---|
 | 55 |   global = global_;
 | 
|---|
 | 56 |   extent = extent_;
 | 
|---|
 | 57 |   iterators.SetSubgrids(local_);
 | 
|---|
 | 58 |   Allocate();
 | 
|---|
 | 59 | }
 | 
|---|
 | 60 | 
 | 
|---|
 | 61 | void TempGrid::SetProperties(const Index& size, const Index& halo_size,
 | 
|---|
 | 62 |                    const Vector& spatial_begin, const Vector& spatial_end)
 | 
|---|
 | 63 | {
 | 
|---|
 | 64 |   global.LocalBegin() = 0;
 | 
|---|
 | 65 |   global.LocalEnd() = size;
 | 
|---|
 | 66 |   global.LocalSize() = size;
 | 
|---|
 | 67 |   global.GlobalBegin() = 0;
 | 
|---|
 | 68 |   global.GlobalEnd() = size;
 | 
|---|
 | 69 |   global.GlobalSize() = size;
 | 
|---|
 | 70 |   global.GlobalBeginFinest() = 0;
 | 
|---|
 | 71 |   global.GlobalEndFinest() = size;
 | 
|---|
 | 72 |   global.GlobalSizeFinest() = size;
 | 
|---|
 | 73 |   global.BoundaryType() = BTUndefined;
 | 
|---|
 | 74 | 
 | 
|---|
 | 75 |   local.Begin() = halo_size;
 | 
|---|
 | 76 |   local.End() = this->local.Begin() + size;
 | 
|---|
 | 77 |   local.Size() = size;
 | 
|---|
 | 78 |   local.SizeTotal() = size + 2 * halo_size;
 | 
|---|
 | 79 |   local.HaloBegin1() = 0;
 | 
|---|
 | 80 |   local.HaloEnd1() = halo_size;
 | 
|---|
 | 81 |   local.HaloSize1() = halo_size;
 | 
|---|
 | 82 |   local.HaloBegin2() = this->local.End();
 | 
|---|
 | 83 |   local.HaloEnd2() = this->local.HaloBegin2() = halo_size;
 | 
|---|
 | 84 |   local.HaloSize2() = halo_size;
 | 
|---|
 | 85 |   local.BoundaryBegin1() = 0;
 | 
|---|
 | 86 |   local.BoundaryEnd1() = 0;
 | 
|---|
 | 87 |   local.BoundarySize1() = 0;
 | 
|---|
 | 88 |   local.BoundaryBegin2() = 0;
 | 
|---|
 | 89 |   local.BoundaryEnd2() = 0;
 | 
|---|
 | 90 |   local.BoundarySize2() = 0;
 | 
|---|
 | 91 | 
 | 
|---|
 | 92 |   extent.Begin() = spatial_begin;
 | 
|---|
 | 93 |   extent.End() = spatial_end;
 | 
|---|
 | 94 |   extent.Size() = spatial_end - spatial_begin;
 | 
|---|
 | 95 |   extent.MeshWidth() = this->extent.Size() / static_cast<Vector>(size-1);
 | 
|---|
 | 96 | 
 | 
|---|
 | 97 |   Allocate();
 | 
|---|
 | 98 | }
 | 
|---|
 | 99 | 
 | 
|---|
 | 100 | void TempGrid::SetPropertiesToFiner(const Grid& grid, const Boundary& boundary)
 | 
|---|
 | 101 | {
 | 
|---|
 | 102 |   assert(grid.Father() != NULL);
 | 
|---|
 | 103 |   assert(grid.Level() < grid.Father()->MaxLevel());
 | 
|---|
 | 104 | 
 | 
|---|
 | 105 |   const Grid& grid_finer = (*grid.Father())(grid.Level()+1);
 | 
|---|
 | 106 | 
 | 
|---|
 | 107 |   /*
 | 
|---|
 | 108 |    * Set global grid attributes
 | 
|---|
 | 109 |    */
 | 
|---|
 | 110 | 
 | 
|---|
 | 111 |   level = grid.Level() + 1;
 | 
|---|
 | 112 | 
 | 
|---|
 | 113 |   global.GlobalBegin() = grid_finer.Global().GlobalBegin();
 | 
|---|
 | 114 |   global.GlobalEnd() = grid_finer.Global().GlobalEnd();
 | 
|---|
 | 115 |   global.GlobalSize() = grid_finer.Global().GlobalSize();
 | 
|---|
 | 116 | 
 | 
|---|
 | 117 |   global.GlobalBeginFinest() = grid_finer.Global().GlobalBeginFinest();
 | 
|---|
 | 118 |   global.GlobalEndFinest() = grid_finer.Global().GlobalEndFinest();
 | 
|---|
 | 119 |   global.GlobalSizeFinest() = grid_finer.Global().GlobalSizeFinest();
 | 
|---|
 | 120 | 
 | 
|---|
 | 121 |   global.BoundaryType() = grid_finer.Global().BoundaryType();
 | 
|---|
 | 122 | 
 | 
|---|
 | 123 |   global.LocalBegin() = (2*grid.Global().LocalBegin()).Clamp(global.GlobalBegin(), global.GlobalEnd());
 | 
|---|
 | 124 |   global.LocalEnd() = (2*grid.Global().LocalEnd()).Clamp(global.GlobalBegin(), global.GlobalEnd());
 | 
|---|
 | 125 |   global.LocalSize() = global.LocalEnd() - global.LocalBegin();
 | 
|---|
 | 126 | 
 | 
|---|
 | 127 |   if (global.LocalSize().Product() == 0) {
 | 
|---|
 | 128 |     global.LocalBegin() = 0;
 | 
|---|
 | 129 |     global.LocalEnd() = 0;
 | 
|---|
 | 130 |     global.LocalSize() = 0;
 | 
|---|
 | 131 |     global.BoundaryType() = EmptyGrid;
 | 
|---|
 | 132 |   }
 | 
|---|
 | 133 | 
 | 
|---|
 | 134 |   local.Begin() = 0;
 | 
|---|
 | 135 |   local.End() = global.LocalSize();
 | 
|---|
 | 136 |   local.Size() = global.LocalSize();
 | 
|---|
 | 137 |   local.SizeTotal() = global.LocalSize();
 | 
|---|
 | 138 | 
 | 
|---|
 | 139 |   local.HaloBegin1() = 0;
 | 
|---|
 | 140 |   local.HaloEnd1() = 0;
 | 
|---|
 | 141 |   local.HaloBegin2() = 0;
 | 
|---|
 | 142 |   local.HaloEnd2() = 0;
 | 
|---|
 | 143 | 
 | 
|---|
 | 144 |   local.BoundaryBegin1() = 0;
 | 
|---|
 | 145 |   local.BoundaryEnd1() = 0;
 | 
|---|
 | 146 |   local.BoundaryBegin2() = 0;
 | 
|---|
 | 147 |   local.BoundaryEnd2() = 0;
 | 
|---|
 | 148 | 
 | 
|---|
 | 149 |   for (int i=0; i<3; ++i) {
 | 
|---|
 | 150 | 
 | 
|---|
 | 151 |     if (grid.Local().HaloSize1()[i] > 0) {
 | 
|---|
 | 152 |       local.Begin()[i] += grid.Local().HaloSize1()[i];
 | 
|---|
 | 153 |       local.End()[i] += grid.Local().HaloSize1()[i];
 | 
|---|
 | 154 |       local.HaloEnd1()[i] = grid.Local().HaloSize1()[i];
 | 
|---|
 | 155 |     }
 | 
|---|
 | 156 | 
 | 
|---|
 | 157 |     if (grid.Local().BoundarySize1()[i] > 0 && global.BoundaryType() != LocallyRefined) {
 | 
|---|
 | 158 |       local.Size()[i] -= grid.Local().BoundarySize1()[i];
 | 
|---|
 | 159 |       local.Begin()[i] += grid.Local().BoundarySize1()[i];
 | 
|---|
 | 160 |       local.BoundaryEnd1()[i] = grid.Local().BoundarySize1()[i];
 | 
|---|
 | 161 |     }
 | 
|---|
 | 162 | 
 | 
|---|
 | 163 |     if (grid.Local().HaloSize2()[i] > 0) {
 | 
|---|
 | 164 |       local.HaloBegin2()[i] = local.End()[i];
 | 
|---|
 | 165 |       local.HaloEnd2()[i] = local.End()[i] + grid.Local().HaloSize2()[i];
 | 
|---|
 | 166 |     }
 | 
|---|
 | 167 | 
 | 
|---|
 | 168 |     if (grid.Local().BoundarySize2()[i] > 0 && global.BoundaryType() != LocallyRefined) {
 | 
|---|
 | 169 |       local.Size()[i] -= grid.Local().BoundarySize2()[i];
 | 
|---|
 | 170 |       local.End()[i] -= grid.Local().BoundarySize2()[i];
 | 
|---|
 | 171 |       local.BoundaryBegin2()[i] = local.End()[i];
 | 
|---|
 | 172 |       local.BoundaryEnd2()[i] = local.End()[i] + grid.Local().BoundarySize2()[i];
 | 
|---|
 | 173 |     }
 | 
|---|
 | 174 | 
 | 
|---|
 | 175 |   }
 | 
|---|
 | 176 | 
 | 
|---|
 | 177 |   local.HaloSize1() = local.HaloEnd1() - local.HaloBegin1();
 | 
|---|
 | 178 |   local.HaloSize2() = local.HaloEnd2() - local.HaloBegin2();
 | 
|---|
 | 179 |   local.BoundarySize1() = local.BoundaryEnd1() - local.BoundaryBegin1();
 | 
|---|
 | 180 |   local.BoundarySize2() = local.BoundaryEnd2() - local.BoundaryBegin2();
 | 
|---|
 | 181 |   local.Size() = local.End() - local.Begin();
 | 
|---|
 | 182 |   local.SizeTotal() = local.Size() +
 | 
|---|
 | 183 |     local.HaloSize1() + local.HaloSize2() +
 | 
|---|
 | 184 |     local.BoundarySize1() + local.BoundarySize2();
 | 
|---|
 | 185 | 
 | 
|---|
 | 186 |   extent = grid_finer.Extent();
 | 
|---|
 | 187 | 
 | 
|---|
 | 188 |   iterators.SetSubgrids(local);
 | 
|---|
 | 189 | 
 | 
|---|
 | 190 |   Allocate();
 | 
|---|
 | 191 | }
 | 
|---|
 | 192 | 
 | 
|---|
 | 193 | void TempGrid::SetPropertiesToCoarser(const Grid& grid, const Boundary& boundary)
 | 
|---|
 | 194 | {
 | 
|---|
 | 195 |   assert(grid.Father() != NULL);
 | 
|---|
 | 196 |   assert(grid.Level() > grid.Father()->MinLevel());
 | 
|---|
 | 197 | 
 | 
|---|
 | 198 |   const Grid& grid_coarser = (*grid.Father())(grid.Level()-1);
 | 
|---|
 | 199 | 
 | 
|---|
 | 200 |   level = grid.Level() - 1;
 | 
|---|
 | 201 | 
 | 
|---|
 | 202 |   global.GlobalBegin() = grid.Global().GlobalBegin();
 | 
|---|
 | 203 |   global.GlobalEnd() = grid.Global().GlobalEnd();
 | 
|---|
 | 204 |   global.GlobalSize() = grid.Global().GlobalSize();
 | 
|---|
 | 205 | 
 | 
|---|
 | 206 |   global.GlobalBeginFinest() = grid.Global().GlobalBeginFinest();
 | 
|---|
 | 207 |   global.GlobalEndFinest() = grid.Global().GlobalEndFinest();
 | 
|---|
 | 208 |   global.GlobalSizeFinest() = grid.Global().GlobalSizeFinest();
 | 
|---|
 | 209 | 
 | 
|---|
 | 210 |   global.BoundaryType() = grid_coarser.Global().BoundaryType();
 | 
|---|
 | 211 | 
 | 
|---|
 | 212 |   global.LocalBegin() = grid.Global().LocalBegin();
 | 
|---|
 | 213 |   global.LocalEnd() = grid.Global().LocalEnd();
 | 
|---|
 | 214 |   GridIndexTranslations::GlobalFineToCoarse(global.LocalBegin(), global.LocalEnd());
 | 
|---|
 | 215 |   global.LocalSize() = global.LocalEnd() - global.LocalBegin();
 | 
|---|
 | 216 | 
 | 
|---|
 | 217 |   if (global.LocalSize().Product() == 0) {
 | 
|---|
 | 218 |     global.LocalBegin() = 0;
 | 
|---|
 | 219 |     global.LocalEnd() = 0;
 | 
|---|
 | 220 |     global.LocalSize() = 0;
 | 
|---|
 | 221 |     global.BoundaryType() = EmptyGrid;
 | 
|---|
 | 222 |   }
 | 
|---|
 | 223 | 
 | 
|---|
 | 224 |   local.SizeTotal() = global.LocalSize();
 | 
|---|
 | 225 |   local.Size() = global.LocalSize();
 | 
|---|
 | 226 |   local.Begin() = 0;
 | 
|---|
 | 227 |   local.End() = global.LocalSize();
 | 
|---|
 | 228 | 
 | 
|---|
 | 229 |   for (int i=0; i<3; ++i) {
 | 
|---|
 | 230 | 
 | 
|---|
 | 231 |     if (grid.Local().HaloSize1()[i] > 0) {
 | 
|---|
 | 232 |       local.SizeTotal()[i] += grid.Local().HaloSize1()[i];
 | 
|---|
 | 233 |       local.Begin()[i] += grid.Local().HaloSize1()[i];
 | 
|---|
 | 234 |       local.End()[i] += grid.Local().HaloSize1()[i];
 | 
|---|
 | 235 |       local.HaloBegin1()[i] = 0;
 | 
|---|
 | 236 |       local.HaloEnd1()[i] = grid.Local().HaloSize1()[i];
 | 
|---|
 | 237 |     }else {
 | 
|---|
 | 238 |       local.HaloBegin1()[i] = 0;
 | 
|---|
 | 239 |       local.HaloEnd1()[i] = 0;
 | 
|---|
 | 240 |     }
 | 
|---|
 | 241 | 
 | 
|---|
 | 242 |     if (grid.Local().BoundarySize1()[i]> 0) {
 | 
|---|
 | 243 |       local.Size()[i] -= grid.Local().BoundarySize1()[i];
 | 
|---|
 | 244 |       local.Begin()[i] += grid.Local().BoundarySize1()[i];
 | 
|---|
 | 245 |       local.BoundaryBegin1()[i] = 0;
 | 
|---|
 | 246 |       local.BoundaryEnd1()[i] = grid.Local().BoundarySize1()[i];
 | 
|---|
 | 247 |     }else {
 | 
|---|
 | 248 |       local.BoundaryBegin1()[i] = 0;
 | 
|---|
 | 249 |       local.BoundaryEnd1()[i] = 0;
 | 
|---|
 | 250 |     }
 | 
|---|
 | 251 | 
 | 
|---|
 | 252 |     if (grid.Local().HaloSize2()[i] > 0) {
 | 
|---|
 | 253 |       local.SizeTotal()[i] += grid.Local().HaloSize2()[i];
 | 
|---|
 | 254 |       local.HaloBegin2()[i] = local.End()[i];
 | 
|---|
 | 255 |       local.HaloEnd2()[i] = local.End()[i] + grid.Local().HaloSize2()[i];
 | 
|---|
 | 256 |     }else {
 | 
|---|
 | 257 |       local.HaloBegin2()[i] = 0;
 | 
|---|
 | 258 |       local.HaloEnd2()[i] = 0;
 | 
|---|
 | 259 |     }
 | 
|---|
 | 260 | 
 | 
|---|
 | 261 |     if (grid.Local().BoundarySize2()[i] > 0) {
 | 
|---|
 | 262 |       local.Size()[i] -= grid.Local().BoundarySize2()[i];
 | 
|---|
 | 263 |       local.End()[i] -= grid.Local().BoundarySize2()[i];
 | 
|---|
 | 264 |       local.BoundaryBegin2()[i] = local.End()[i];
 | 
|---|
 | 265 |       local.BoundaryEnd2()[i] = local.End()[i] + grid.Local().BoundarySize2()[i];
 | 
|---|
 | 266 |     }else {
 | 
|---|
 | 267 |       local.BoundaryBegin2()[i] = 0;
 | 
|---|
 | 268 |       local.BoundaryEnd2()[i] = 0;
 | 
|---|
 | 269 |     }
 | 
|---|
 | 270 | 
 | 
|---|
 | 271 |   }
 | 
|---|
 | 272 | 
 | 
|---|
 | 273 |   local.HaloSize1() = local.HaloEnd1() - local.HaloBegin1();
 | 
|---|
 | 274 |   local.HaloSize2() = local.HaloEnd2() - local.HaloBegin2();
 | 
|---|
 | 275 |   local.BoundarySize1() = local.BoundaryEnd1() - local.BoundaryBegin1();
 | 
|---|
 | 276 |   local.BoundarySize2() = local.BoundaryEnd2() - local.BoundaryBegin2();
 | 
|---|
 | 277 | 
 | 
|---|
 | 278 |   Extent().Size() = grid.Extent().Size();
 | 
|---|
 | 279 |   Extent().Begin() = grid.Extent().Begin();
 | 
|---|
 | 280 |   Extent().End() = grid.Extent().End();
 | 
|---|
 | 281 |   Extent().MeshWidth() = 2.0 * grid.Extent().MeshWidth();
 | 
|---|
 | 282 | 
 | 
|---|
 | 283 |   iterators.SetSubgrids(local);
 | 
|---|
 | 284 |   Allocate();
 | 
|---|
 | 285 | }
 | 
|---|
 | 286 | 
 | 
|---|
 | 287 | void TempGrid::ImportFromResidual(Grid& sol, Grid& rhs)
 | 
|---|
 | 288 | {
 | 
|---|
 | 289 |   Grid::iterator iter;
 | 
|---|
 | 290 | 
 | 
|---|
 | 291 |   const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol);
 | 
|---|
 | 292 |   const Stencil& A = MG::GetDiscretization()->GetStencil();
 | 
|---|
 | 293 | 
 | 
|---|
 | 294 |   this->Clear();
 | 
|---|
 | 295 | 
 | 
|---|
 | 296 |   MG::GetComm()->CommToGhosts(sol);
 | 
|---|
 | 297 | 
 | 
|---|
 | 298 |   for (iter=Iterators().Local().Begin(); iter!=Iterators().Local().End(); ++iter)
 | 
|---|
 | 299 |     (*this)(*iter) = rhs.GetVal(*iter) - prefactor * A.Apply(sol, *iter);
 | 
|---|
 | 300 | 
 | 
|---|
 | 301 |   this->ClearBoundary();
 | 
|---|
 | 302 | }
 | 
|---|
 | 303 | 
 | 
|---|
 | 304 | void TempGrid::Allocate()
 | 
|---|
 | 305 | {
 | 
|---|
 | 306 |   const int size = local.SizeTotal().Product();
 | 
|---|
 | 307 | 
 | 
|---|
 | 308 |   if (size > size_max) {
 | 
|---|
 | 309 |     size_max = size;
 | 
|---|
 | 310 |     delete [] grid;
 | 
|---|
 | 311 |     grid = new vmg_float[size];
 | 
|---|
 | 312 |   }
 | 
|---|
 | 313 | }
 | 
|---|
 | 314 | 
 | 
|---|
 | 315 | TempGrid::TempGrid() :
 | 
|---|
 | 316 |   size_max(0)
 | 
|---|
 | 317 | {
 | 
|---|
 | 318 | }
 | 
|---|
 | 319 | 
 | 
|---|
 | 320 | TempGrid::TempGrid(const Grid& rhs) :
 | 
|---|
 | 321 |   size_max(0)
 | 
|---|
 | 322 | {
 | 
|---|
 | 323 |   SetProperties(rhs);
 | 
|---|
 | 324 | }
 | 
|---|
 | 325 | 
 | 
|---|
 | 326 | TempGrid::TempGrid(const GlobalIndices& global, const LocalIndices& local, const SpatialExtent& extent) :
 | 
|---|
 | 327 |   size_max(0)
 | 
|---|
 | 328 | {
 | 
|---|
 | 329 |   SetProperties(global, local, extent);
 | 
|---|
 | 330 | }
 | 
|---|
 | 331 | 
 | 
|---|
 | 332 | TempGrid::TempGrid(const Index& size, const Index& halo_size,
 | 
|---|
 | 333 |                    const Vector& spatial_begin, const Vector& spatial_end) :
 | 
|---|
 | 334 |   size_max(0)
 | 
|---|
 | 335 | {
 | 
|---|
 | 336 |   SetProperties(size, halo_size, spatial_begin, spatial_end);
 | 
|---|
 | 337 | }
 | 
|---|
 | 338 | 
 | 
|---|
 | 339 | TempGrid::~TempGrid()
 | 
|---|
 | 340 | {
 | 
|---|
 | 341 | }
 | 
|---|