AdaptiveGrid.hpp 18.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#pragma once

#include <algorithm>
#include <list>
#include <memory>
#include <mutex>
#include <shared_mutex>
#include <type_traits>
#include <utility>

Praetorius, Simon's avatar
Praetorius, Simon committed
11
#include <dune/common/hybridutilities.hh>
12
#include <dune/common/timer.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
13
14
#include <dune/common/version.hh>
#include <dune/common/parallel/mpihelper.hh>
15
16

#include <dune/geometry/type.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
17
18
19
20
#include <dune/grid/common/backuprestore.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/grid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
21

Praetorius, Simon's avatar
Praetorius, Simon committed
22
23
#include <amdis/common/Concepts.hpp>
#include <amdis/common/DefaultGridView.hpp>
24
25
26
27
28
29
30
31
32
33
34
#include <amdis/common/SharedPtr.hpp>
#include <amdis/Observer.hpp>
#include <amdis/Output.hpp>

namespace AMDiS
{
  namespace event
  {
    /** Event generated from an AdaptiveGrid when calling preAdapt(). It contains the return value
     *  of preAdapt() as its 'mightCoarsen' member and is passed to registered observers after
     *  calling preAdapt on the underlying grid.
Praetorius, Simon's avatar
Praetorius, Simon committed
35
     **/
36
37
38
39
40
    struct preAdapt { bool mightCoarsen; };

    /** Event generated from an AdaptiveGrid when calling adapt(). Its 'adapted' member contains the
     *  value true if either preAdapt() or adapt() returned true. This event is passed to registered
     *  observers after calling adapt on the underlying grid.
Praetorius, Simon's avatar
Praetorius, Simon committed
41
     **/
42
43
44
45
    struct adapt { bool adapted; };

    /** Event generated from an AdaptiveGrid when calling postAdapt().This event is passed to
     *  registered observers after calling postAdapt on the underlying grid.
Praetorius, Simon's avatar
Praetorius, Simon committed
46
     **/
47
48
49
    struct postAdapt {};
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
50
51
52
53
54
55
56
57
58

  // forward declaration
  template <class HostGrid>
  class AdaptiveGridFamily;


  /// \brief Wrapper class for Dune-grids that allows automatic signalling of events
  /// during grid adaptation.
  /**
59
60
   *  Calls to grid.preAdapt(), grid.adapt() and grid.postAdapt() need to be replaced by calls to
   *  the corresponding methods of this class to use the automatic update functionality.
Praetorius, Simon's avatar
Praetorius, Simon committed
61
62
63
64
   *
   * \tparam HG  Host grid to be wrapped. Must implement the dune grid interface.
   **/
  template <class HG>
65
  class AdaptiveGrid
Praetorius, Simon's avatar
Praetorius, Simon committed
66
67
68
      : public Dune::GridDefaultImplementation<
          HG::dimension, HG::dimensionworld, typename HG::ctype, AdaptiveGridFamily<HG> >
      , public Signals<event::preAdapt, event::adapt, event::postAdapt>
69
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
70
    using Self = AdaptiveGrid<HG>;
71
72
73

  public:

Praetorius, Simon's avatar
Praetorius, Simon committed
74
75
76
77
78
    using HostGrid   = HG;
    using GridFamily = AdaptiveGridFamily<HG>;
    using Traits     = typename GridFamily::Traits;
    using Element    = typename Traits::template Codim<0>::Entity;

79

Praetorius, Simon's avatar
Praetorius, Simon committed
80
81
82
83
84
85
  public:

    template <class HostGrid_,
      Dune::disableCopyMove<Self, HostGrid_> = 0>
    explicit AdaptiveGrid(HostGrid_&& hostGrid)
      : hostGrid_(wrap_or_share(FWD(hostGrid)))
86
87
    {}

Praetorius, Simon's avatar
Praetorius, Simon committed
88
89
    /// Return the underlying grid
    std::shared_ptr<HostGrid> const& hostGrid() const { return hostGrid_; }
90

Praetorius, Simon's avatar
Praetorius, Simon committed
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

  public:

    /// Grid iterators
    /// @{

    /// Iterator to first entity of given codim on level
    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto lbegin(int level) const { return hostGrid_->levelGridView(level).template begin<codim,pt>(); }

    /// one past the end on this level
    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto lend(int level) const { return hostGrid_->levelGridView(level).template end<codim,pt>(); }


    /// Iterator to first leaf entity of given codim
    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto leafbegin() const { return hostGrid_->leafGridView().template begin<codim,pt>(); }

    /// One past the end of the sequence of leaf entities
    template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
    auto leafend() const { return hostGrid_->leafGridView().template end<codim,pt>(); }


    /// Obtain begin intersection iterator with respect to the level GridView
    auto ilevelbegin(Element const& e) const { return hostGrid_->levelGridView(e.level()).ibegin(e); }

    /// Obtain end intersection iterator with respect to the level GridView
    auto ilevelend(Element const& e) const { return hostGrid_->levelGridView(e.level()).iend(e); }

    /// Obtain begin intersection iterator with respect to the leaf GridView
    auto ileafbegin(Element const& e) const { return hostGrid_->leafGridView().ibegin(e); }

    /// Obtain end intersection iterator with respect to the leaf GridView
    auto ileafend(Element const& e) const { return hostGrid_->leafGridView().iend(e); }

    /// @}


    /// Size methods
    /// @{

    /// Return maximum level defined in this grid.
134
135
    int maxLevel() const { return hostGrid_->maxLevel(); }

Praetorius, Simon's avatar
Praetorius, Simon committed
136
    /// Number of grid entities per level and codim
137
138
139
140
141
142
143
144
145
146
147
    int size(int level, int codim) const { return hostGrid_->size(level, codim); }

    /// Return number of leaf entities of a given codim in this process
    int size(int codim) const { return hostGrid_->size(codim); }

    /// Return number of entities per level and geometry type in this process
    int size(int level, Dune::GeometryType type) const { return hostGrid_->size(level, type); }

    /// Return number of leaf entities per geometry type in this process
    int size(Dune::GeometryType type) const { return hostGrid_->size(type); }

Praetorius, Simon's avatar
Praetorius, Simon committed
148
149
    /// Returns the number of boundary segments within the macro grid
    std::size_t numBoundarySegments() const { return hostGrid_->numBoundarySegments(); }
150

Praetorius, Simon's avatar
Praetorius, Simon committed
151
    /// @}
152
153


Praetorius, Simon's avatar
Praetorius, Simon committed
154
155
    /// Access to index and id sets
    /// @{
156

Praetorius, Simon's avatar
Praetorius, Simon committed
157
158
    /// Return const reference to the grids global id set
    auto const& globalIdSet() const { return hostGrid_->globalIdSet(); }
159

Praetorius, Simon's avatar
Praetorius, Simon committed
160
161
    /// Return const reference to the host grids local id set
    auto const& localIdSet() const { return hostGrid_->localIdSet(); }
162

Praetorius, Simon's avatar
Praetorius, Simon committed
163
164
165
166
167
168
169
170
171
172
173
    /// Return const reference to the host grids level index set for level level
    auto const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); }

    /// Return const reference to the host grids leaf index set
    auto const& leafIndexSet() const { return hostGrid_->leafIndexSet(); }

    /// @}


    /// Adaptivity and grid refinement
    /// @{
174
175
176
177
178
179
180
181
182
183
184
185
186
187

    /// Refines all grid elements refCount times.
    void globalRefine(int refCount)
    {
      for (int i = 0; i < refCount; ++i) {
        // mark all entities for grid refinement
        for (const auto& element : elements(hostGrid_->leafGridView()))
          hostGrid_->mark(1, element);
        preAdapt();
        adapt();
        postAdapt();
      }
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
188
    /// Marks an entity to be refined/coarsened in a subsequent adapt
189
190
191
192
193
194
195
196
197
198
199
200
    bool mark(int refCount, Element const& e) { return hostGrid_->mark(refCount, e); }

    /// Return refinement mark for entity
    int getMark(Element const& e) const { return hostGrid_->getMark(e); }

    /// Prepare the grid for adaptation and notify observers of the preAdapt event
    bool preAdapt()
    {
      Dune::Timer t;
      mightCoarsen_ = hostGrid_->preAdapt();
      Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1);
      this->notify(event::preAdapt{mightCoarsen_});
Praetorius, Simon's avatar
Praetorius, Simon committed
201
      info(2,"AdaptiveGrid::preAdapt needed {} seconds", t.elapsed());
202
203
204
205
206
207
208
209
210
211
      return mightCoarsen_;
    }

    /// Adapt the grid and notify observers of the adapt event
    bool adapt()
    {
      Dune::Timer t;
      refined_ = hostGrid_->adapt();
      Dune::MPIHelper::getCollectiveCommunication().max(&refined_, 1);
      this->notify(event::adapt{mightCoarsen_ || refined_});
Praetorius, Simon's avatar
Praetorius, Simon committed
212
      info(2,"AdaptiveGrid::adapt needed {} seconds", t.elapsed());
213
214
215
216
217
218
219
220
221
222
      return refined_;
    }

    /// Perform cleanup after grid adaptation and notify observers of the postAdapt event
    void postAdapt()
    {
      Dune::Timer t;
      hostGrid_->postAdapt();
      this->notify(event::postAdapt{});
      changeIndex_++;
Praetorius, Simon's avatar
Praetorius, Simon committed
223
      info(2,"AdaptiveGrid::postAdapt needed {} seconds", t.elapsed());
224
225
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
226
227
    /// Returns the grid change index, see \ref changeIndex.
    unsigned long changeIndex() const
228
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
229
      return changeIndex_;
230
231
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
232
233
234
235
236
237
238
239
240
241
    // @}


    /// Parallel data distribution and communication
    /// @{

    /// Return const reference to a collective communication object.
    auto const& comm() const { return hostGrid_->comm(); }

    /// Communicate data of level gridView
242
    template <class DataHandle>
Praetorius, Simon's avatar
Praetorius, Simon committed
243
244
    void communicate(DataHandle& handle, Dune::InterfaceType iftype,
                     Dune::CommunicationDirection dir, int level) const
245
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
246
      hostGrid_->levelGridView(level).communicate(handle,iftype,dir);
247
248
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
249
250
251
252
    /// Communicate data of leaf gridView
    template <class DataHandle>
    void communicate(DataHandle& handle, Dune::InterfaceType iftype,
                     Dune::CommunicationDirection dir) const
253
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
254
      hostGrid_->leafGridView().communicate(handle,iftype,dir);
255
256
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
257
258
259
260
261
262
263
    /// \brief Calls loadBalance on the underlying grid.
    /**
     *  Re-balances the load each process has to handle for a parallel grid.
     *
     *  \return true if the grid has changed.
     **/
    bool loadBalance()
264
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
265
266
267
268
      return Dune::Hybrid::ifElse(
        /* if */ std::is_convertible<decltype(std::declval<HG>().loadBalance()), bool>{},
        [&](auto id) { return id(hostGrid_)->loadBalance(); },
        [&](auto id) { id(hostGrid_)->loadBalance(); return true; });
269
270
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
271
272
273
274
275
276
277
278
279
280
    /// \brief Calls loadBalance(handle) on the underlying grid.
    /**
     *  Re-balances the load each process has to handle for a parallel grid and moves the data.
     *
     *  \param data  A data handle telling the method which data should be communicated
     *               and how. Has to adhere to the interface described by Dune::CommDataHandleIf.
     *  \return true if the grid has changed.
     **/
    template <class DataHandle>
    bool loadBalance(DataHandle& handle)
281
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
282
283
284
285
      return Dune::Hybrid::ifElse(
        /* if */ std::is_convertible<decltype(std::declval<HG>().loadBalance(handle)), bool>{},
        [&](auto id) { return id(hostGrid_)->loadBalance(handle); },
        [&](auto id) { id(hostGrid_)->loadBalance(handle); return true; });
286
287
288
    }


Praetorius, Simon's avatar
Praetorius, Simon committed
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
    /// Return size of the overlap region for a given codim on the level grid view.
    int overlapSize(int level, int codim) const { return hostGrid_->levelGridView(level).overlapSize(codim); }

    /// Return size of the overlap region for a given codim on the leaf grid view.
    int overlapSize(int codim) const { return hostGrid_->leafGridView().overlapSize(codim); }

    /// Return size of the ghost region for a given codim on the level grid view.
    int ghostSize(int level, int codim) const { return hostGrid_->levelGridView(level).ghostSize(codim); }

    /// Return size of the ghost region for a given codim on the leaf grid view.
    int ghostSize(int codim) const { return hostGrid_->leafGridView().ghostSize(codim); }

    /// @}


    /// Obtain Entity from EntitySeed of the HostGrid.
    template <class EntitySeed>
    auto entity(EntitySeed const& seed) const { return hostGrid_->entity(seed); }

308
309
310
311
312
313

  private:

    /// The underlying grid implementation
    std::shared_ptr<HostGrid> hostGrid_;

Praetorius, Simon's avatar
Praetorius, Simon committed
314
315
    /// Flag set during \ref preAdapt(), indicating whether any element might be
    /// coarsened in \ref adapt()
316
317
318
319
320
    bool mightCoarsen_ = false;

    /// Flag set during \ref adapt() indicating that at least one entity was refined
    bool refined_ = false;

Praetorius, Simon's avatar
Praetorius, Simon committed
321
322
    /// This index is incremented every time the grid is changed, e.g. by refinement
    /// or coarsening.
323
324
325
    unsigned long changeIndex_ = 0;
  };

Praetorius, Simon's avatar
Praetorius, Simon committed
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374

  template <class HostGrid>
  class AdaptiveGridFamily
  {
  public:
    struct Traits : HostGrid::Traits
    {
      using Grid = AdaptiveGrid<HostGrid>;
      using LeafGridView = Dune::GridView< AMDiS::DefaultLeafGridViewTraits<const Grid> >;
      using LevelGridView = Dune::GridView< AMDiS::DefaultLevelGridViewTraits<const Grid> >;
    };
  };


  /// Return a change index of a grid that is not an \ref AdaptiveGrid
  template <class HostGrid>
  unsigned long changeIndex(HostGrid const& /*hostGrid*/)
  {
    return 0;
  }

  /// Return a change index of an \ref AdaptiveGrid
  template <class HostGrid>
  unsigned long changeIndex(AdaptiveGrid<HostGrid> const& grid)
  {
    return grid.changeIndex();
  }


  namespace Impl
  {
    template <class HostGrid>
    struct AdaptiveGridImpl
    {
      using type = AdaptiveGrid<HostGrid>;
    };

    template <class HostGrid>
    struct AdaptiveGridImpl<AdaptiveGrid<HostGrid>>
    {
      using type = AdaptiveGrid<HostGrid>;
    };
  }

  /// Always returning an \ref AdaptiveGrid. Returns the grid itself if it is
  /// already an AdaptiveGrid.
  template <class HostGrid>
  using AdaptiveGrid_t = typename Impl::AdaptiveGridImpl<HostGrid>::type;

375
376

} // end namespace AMDiS
Praetorius, Simon's avatar
Praetorius, Simon committed
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559


namespace Dune
{
  /// Specialization of a \ref GridFactory to \ref AdaptiveGrid.
  /// Provide a generic factory class for unstructured grids.
  template <class HostGrid>
  class GridFactory<AMDiS::AdaptiveGrid<HostGrid> >
      : public GridFactoryInterface<AMDiS::AdaptiveGrid<HostGrid> >
  {
    using Self = GridFactory;
    using GridType = AMDiS::AdaptiveGrid<HostGrid>;
    using HostGridFactory = GridFactory<HostGrid>;

  public:

    using ctype = typename GridType::ctype;

    enum { dim = GridType::dimension };
    enum { dimworld = GridType::dimensionworld };

    template <class... Args,
      Dune::disableCopyMove<Self, Args...> = 0>
    GridFactory(Args&&... args)
      : hostFactory_(FWD(args)...)
    {}

    /// Insert a vertex into the coarse grid
    void insertVertex(FieldVector<ctype,dimworld> const& pos) override
    {
      hostFactory_.insertVertex(pos);
    }

    template <class F, class... Args>
    using HasInsertElement = decltype(std::declval<F>().insertElement(std::declval<Args>()...));

    /// Insert an element into the coarse grid
    void insertElement(GeometryType const& type,
                       std::vector<unsigned int> const& vertices) override
    {
      hostFactory_.insertElement(type, vertices);
    }

    using ElementParametrizationType = std::shared_ptr<VirtualFunction<FieldVector<ctype,dim>,FieldVector<ctype,dimworld> > >;

    /// Insert a parametrized element into the coarse grid
    void insertElement(GeometryType const& type,
                       std::vector<unsigned int> const& vertices,
                       ElementParametrizationType const& elementParametrization) override
    {
      using A0 = GeometryType;
      using A1 = std::vector<unsigned int>;
      using A2 = ElementParametrizationType;
      Hybrid::ifElse(Std::is_detected<HasInsertElement, HostGridFactory, A0,A1,A2>{},
        [&](auto id) { id(hostFactory_).insertElement(type, vertices, elementParametrization); },
        [&](auto id) { AMDiS::error_exit("insertElement() not implemented for HostGrid type."); });
    }

    template <class F, class... Args>
    using HasInsertBoundarySegment = decltype(std::declval<F>().insertBoundarySegment(std::declval<Args>()...));

    //// Insert a boundary segment
    void insertBoundarySegment(std::vector<unsigned int> const& vertices) override
    {
      hostFactory_.insertBoundarySegment(vertices);
    }

    using BoundarySegmentType = std::shared_ptr<BoundarySegment<dim,dimworld> >;

    /// Insert an arbitrarily shaped boundary segment
    void insertBoundarySegment(std::vector<unsigned int> const& vertices,
                               BoundarySegmentType const& boundarySegment) override
    {
      using A0 = std::vector<unsigned int>;
      using A1 = BoundarySegmentType;
      Hybrid::ifElse(Std::is_detected<HasInsertBoundarySegment, HostGridFactory, A0,A1>{},
        [&](auto id) { id(hostFactory_).insertBoundarySegment(vertices, boundarySegment); },
        [&](auto id) { AMDiS::error_exit("insertBoundarySegment() not implemented for HostGrid type."); });
    }

    /// Finalize grid creation and hand over the grid
#if DUNE_VERSION_LT(DUNE_GRID,2,7)
    GridType* createGrid() override
    {
      std::unique_ptr<HostGrid> hostGrid(hostFactory_.createGrid());
      return new GridType(std::move(hostGrid));
    }
#else
    ToUniquePtr<GridType> createGrid() override
    {
      std::unique_ptr<HostGrid> hostGrid(hostFactory_.createGrid());
      return makeToUnique<GridType>(std::move(hostGrid));
    }
#endif

  private:
    HostGridFactory hostFactory_;
  };


  namespace Impl
  {
    template <class HostGrid, bool = Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v>
    class BackupRestoreFacilityImpl {};

    template <class HostGrid>
    class BackupRestoreFacilityImpl<HostGrid,true>
    {
      using Grid = AMDiS::AdaptiveGrid<HostGrid>;
      using HostBackupRestoreFacility = BackupRestoreFacility<HostGrid>;

    public:

      /// Backup the grid to file or stream
      template <class Output>
      static void backup(Grid const& grid, Output const& filename_or_stream)
      {
        HostBackupRestoreFacility::backup(*grid.hostGrid(), filename_or_stream);
      }

      /// Restore the grid from file or stream
      template <class Input>
      static Grid* restore(Input const& filename_or_stream)
      {
        std::unique_ptr<HostGrid> hostGrid(HostBackupRestoreFacility::restore(filename_or_stream));
        return new Grid(std::move(hostGrid));
      }
    };

  } // end namespace Impl


  /// Specialization of \ref BackupRestoreFacility to \ref AdaptiveGrid.
  /// Facility for writing and reading grids.
  template <class HostGrid>
  class BackupRestoreFacility<AMDiS::AdaptiveGrid<HostGrid>>
      : public Impl::BackupRestoreFacilityImpl<HostGrid>
  {
  public:
    using Impl::BackupRestoreFacilityImpl<HostGrid>::BackupRestoreFacilityImpl;
  };


  namespace Capabilities
  {
    template <class HostGrid, int codim>
    struct hasEntity<AMDiS::AdaptiveGrid<HostGrid>, codim>
        : hasEntity<HostGrid,codim>{};

    template <class HostGrid, int codim>
    struct hasEntityIterator<AMDiS::AdaptiveGrid<HostGrid>, codim>
        : hasEntityIterator<HostGrid, codim> {};

    template <class HostGrid>
    struct isLevelwiseConforming<AMDiS::AdaptiveGrid<HostGrid> >
        : isLevelwiseConforming<HostGrid> {};

    template <class HostGrid>
    struct isLeafwiseConforming<AMDiS::AdaptiveGrid<HostGrid> >
        : isLeafwiseConforming<HostGrid> {};

    template <class HostGrid>
    struct hasSingleGeometryType<AMDiS::AdaptiveGrid<HostGrid> >
        : hasSingleGeometryType<HostGrid> {};

    template <class HostGrid, int codim >
    struct canCommunicate<AMDiS::AdaptiveGrid<HostGrid>, codim>
        : canCommunicate<HostGrid, codim> {};

    template <class HostGrid>
    struct hasBackupRestoreFacilities<AMDiS::AdaptiveGrid<HostGrid> >
        : hasBackupRestoreFacilities<HostGrid> {};

    template <class HostGrid>
    struct threadSafe<AMDiS::AdaptiveGrid<HostGrid> >
        : threadSafe<HostGrid> {};

    template <class HostGrid>
    struct viewThreadSafe<AMDiS::AdaptiveGrid<HostGrid> >
        : viewThreadSafe<HostGrid> {};

  } // end namespace Capabilities
} // end namespace Dune