AdaptiveGrid.hpp 18.4 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
#include <amdis/common/SharedPtr.hpp>
#include <amdis/Observer.hpp>
#include <amdis/Output.hpp>

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
30
31
32
33
34
35
36
37
  // forward declaration
  template <class HostGrid>
  class AdaptiveGridFamily;


  /// \brief Wrapper class for Dune-grids that allows automatic signalling of events
  /// during grid adaptation.
  /**
38
39
   *  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
40
41
42
43
   *
   * \tparam HG  Host grid to be wrapped. Must implement the dune grid interface.
   **/
  template <class HG>
44
  class AdaptiveGrid
Praetorius, Simon's avatar
Praetorius, Simon committed
45
46
      : public Dune::GridDefaultImplementation<
          HG::dimension, HG::dimensionworld, typename HG::ctype, AdaptiveGridFamily<HG> >
Praetorius, Simon's avatar
Praetorius, Simon committed
47
      , public Notifier<event::preAdapt, event::adapt, event::postAdapt>
48
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
49
    using Self = AdaptiveGrid<HG>;
50
51
52

  public:

Praetorius, Simon's avatar
Praetorius, Simon committed
53
54
55
56
57
    using HostGrid   = HG;
    using GridFamily = AdaptiveGridFamily<HG>;
    using Traits     = typename GridFamily::Traits;
    using Element    = typename Traits::template Codim<0>::Entity;

58

Praetorius, Simon's avatar
Praetorius, Simon committed
59
60
61
62
63
64
  public:

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

  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.
113
114
    int maxLevel() const { return hostGrid_->maxLevel(); }

Praetorius, Simon's avatar
Praetorius, Simon committed
115
    /// Number of grid entities per level and codim
116
117
118
119
120
121
122
123
124
125
126
    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
127
128
    /// Returns the number of boundary segments within the macro grid
    std::size_t numBoundarySegments() const { return hostGrid_->numBoundarySegments(); }
129

Praetorius, Simon's avatar
Praetorius, Simon committed
130
    /// @}
131
132


Praetorius, Simon's avatar
Praetorius, Simon committed
133
134
    /// Access to index and id sets
    /// @{
135

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
142
143
144
145
146
147
148
149
150
151
152
    /// 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
    /// @{
153
154
155
156
157
158
159
160
161
162
163
164
165
166

    /// 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
167
    /// Marks an entity to be refined/coarsened in a subsequent adapt
168
169
170
171
172
173
174
175
176
177
178
179
    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
180
      info(2,"AdaptiveGrid::preAdapt needed {} seconds", t.elapsed());
181
182
183
184
185
186
187
188
189
190
      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
191
      info(2,"AdaptiveGrid::adapt needed {} seconds", t.elapsed());
192
193
194
195
196
197
198
199
200
201
      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
202
      info(2,"AdaptiveGrid::postAdapt needed {} seconds", t.elapsed());
203
204
    }

205
206
207
208
209
210
    /// Update all registered dependent objects if grid is changed manually
    void update()
    {
      this->notify(event::adapt{true});
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
211
212
    /// Returns the grid change index, see \ref changeIndex.
    unsigned long changeIndex() const
213
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
214
      return changeIndex_;
215
216
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
217
218
219
220
221
222
223
224
225
226
    // @}


    /// 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
227
    template <class DataHandle>
Praetorius, Simon's avatar
Praetorius, Simon committed
228
229
    void communicate(DataHandle& handle, Dune::InterfaceType iftype,
                     Dune::CommunicationDirection dir, int level) const
230
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
231
      hostGrid_->levelGridView(level).communicate(handle,iftype,dir);
232
233
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
234
235
236
237
    /// Communicate data of leaf gridView
    template <class DataHandle>
    void communicate(DataHandle& handle, Dune::InterfaceType iftype,
                     Dune::CommunicationDirection dir) const
238
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
239
      hostGrid_->leafGridView().communicate(handle,iftype,dir);
240
241
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
242
243
244
245
246
247
248
    /// \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()
249
    {
250
251
252
      if constexpr (std::is_convertible_v<decltype(std::declval<HG>().loadBalance()), bool>)
        return hostGrid_->loadBalance();
      else {
253
254
        hostGrid_->loadBalance();
        return true;
255
      }
256
257
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
258
259
260
261
262
263
264
265
266
267
    /// \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)
268
    {
269
270
271
      if constexpr (std::is_convertible_v<decltype(std::declval<HG>().loadBalance(handle)), bool>)
        return hostGrid_->loadBalance(handle);
      else {
272
        hostGrid_->loadBalance(handle);
273
274
        return true;
      }
275
276
277
    }


Praetorius, Simon's avatar
Praetorius, Simon committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
    /// 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); }

297
298
299
300
301
302

  private:

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

Praetorius, Simon's avatar
Praetorius, Simon committed
303
304
    /// Flag set during \ref preAdapt(), indicating whether any element might be
    /// coarsened in \ref adapt()
305
306
307
308
309
    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
310
311
    /// This index is incremented every time the grid is changed, e.g. by refinement
    /// or coarsening.
312
313
314
    unsigned long changeIndex_ = 0;
  };

315
316
317
318
319
  // deduction guide
  template <class HostGrid>
  AdaptiveGrid(HostGrid const&)
    -> AdaptiveGrid<HostGrid>;

Praetorius, Simon's avatar
Praetorius, Simon committed
320
321
322
323
324
325
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

  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;

369
370

} // end namespace AMDiS
Praetorius, Simon's avatar
Praetorius, Simon committed
371
372
373
374
375
376
377
378
379
380
381


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;
382
    using Super = GridFactoryInterface<AMDiS::AdaptiveGrid<HostGrid> >;
Praetorius, Simon's avatar
Praetorius, Simon committed
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
    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);
    }

415
#if DUNE_VERSION_LT(DUNE_GRID,2,8)
Praetorius, Simon's avatar
Praetorius, Simon committed
416
417
418
419
420
421
    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
422
423
424
425
426
427
428
429
430
#else
    using ElementParametrizationType
      = std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dim>)>;

    /// Insert a parametrized element into the coarse grid
    void insertElement(GeometryType const& type,
                       std::vector<unsigned int> const& vertices,
                       ElementParametrizationType elementParametrization) override
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
431
432
433
434
    {
      using A0 = GeometryType;
      using A1 = std::vector<unsigned int>;
      using A2 = ElementParametrizationType;
435
436
437
438
      if constexpr (Std::is_detected<HasInsertElement, HostGridFactory, A0,A1,A2>::value)
        hostFactory_.insertElement(type, vertices, elementParametrization);
      else
        AMDiS::error_exit("insertElement() not implemented for HostGrid type.");
Praetorius, Simon's avatar
Praetorius, Simon committed
439
    }
440
    using Super::insertElement;
Praetorius, Simon's avatar
Praetorius, Simon committed
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458

    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;
459
460
461
462
      if constexpr (Std::is_detected<HasInsertBoundarySegment, HostGridFactory, A0,A1>::value)
        hostFactory_.insertBoundarySegment(vertices, boundarySegment);
      else
        AMDiS::error_exit("insertBoundarySegment() not implemented for HostGrid type.");
Praetorius, Simon's avatar
Praetorius, Simon committed
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
560
561
562
563
564
565
566
567
    }

    /// 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