AdaptiveGrid.hpp 11 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
#pragma once

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

#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>

#include <dune/geometry/type.hh>

#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/SharedPtr.hpp>
#include <amdis/Observer.hpp>
#include <amdis/Output.hpp>

namespace AMDiS
{
  namespace Impl
  {
    template <class Grid, class... Args,
    REQUIRES(std::is_convertible<decltype(std::declval<Grid>().loadBalance()), bool>::value)>
    bool loadBalanceImpl(Grid& grid, Args&&... args)
    {
      return grid.loadBalance(FWD(args)...);
    }

    // Workaround for grids not adhering to the correct interface
    template <class Grid, class... Args,
    REQUIRES(not std::is_convertible<decltype(std::declval<Grid>().loadBalance()), bool>::value)>
    bool loadBalanceImpl(Grid& grid, Args&&... args)
    {
      grid.loadBalance(FWD(args)...);
      return true;
    }
  }

  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.
     */
    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.
     */
    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.
     */
    struct postAdapt {};
  }

  /** Wrapper class for Dune-grids that allows automatic signalling of events during grid adaptation
   *  This class needs to be created after construction of the associated grid by calling the
   *  instance method.
   *  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.
   */
  template <class Grid>
  class AdaptiveGrid
      : public Signals<event::preAdapt, event::adapt, event::postAdapt>
  {
    using Self = AdaptiveGrid<Grid>;
    using Element = typename Grid::template Codim<0>::Entity;
    struct HiddenStruct {};

  public:
    using HostGrid = Grid;
    enum { dimension = HostGrid::dimension };
    enum { dimensionworld = HostGrid::dimensionworld };
    using LeafGridView = typename HostGrid::LeafGridView;
    using LevelGridView = typename HostGrid::LevelGridView;

    template <int cd>
    struct Codim
    {
      using Geometry = typename HostGrid::template Codim<cd>::Geometry;
      using Entity = typename HostGrid::template Codim<cd>::Entity;
    };

    using GlobalIdSet = typename HostGrid::GlobalIdSet;
    using LocalIdSet = typename HostGrid::LocalIdSet;
    using LevelIndexSet = typename HostGrid::LevelIndexSet;
    using LeafIndexSet = typename HostGrid::LeafIndexSet;

    /// Constructor that may only be called indirectly via the instance method
    explicit AdaptiveGrid(std::shared_ptr<HostGrid> grid, HiddenStruct)
      : hostGrid_(std::move(grid))
    {}

    /// Unreachable constructor required to compile the unreachable branch in instance method
    explicit AdaptiveGrid(std::shared_ptr<HostGrid const> grid, HiddenStruct)
    {
      error_exit("Cannot create AdaptiveGrid from const Grid.");
    }

    int maxLevel() const { return hostGrid_->maxLevel(); }

    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); }

    auto numBoundarySegments () const { return hostGrid_->numBoundarySegments(); }

    /// View for a grid level for All_Partition
    LevelGridView levelGridView(int l) const { return hostGrid_->levelGridView(l); }

    /// View for the leaf grid for All_Partition
    LeafGridView leafGridView() const { return hostGrid_->leafGridView(); }

    GlobalIdSet const& globalIdSet() const { return hostGrid_->globalIdSet(); }

    /// return const reference to the host grids local id set
    LocalIdSet const& localIdSet() const { return hostGrid_->localIdSet(); }

    /// return const reference to the host grids level index set for level level
    LevelIndexSet const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); }

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

    /// 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();
      }
    }

    /// Mark entity for refinement
    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_});
      info(2,"preAdapt needed {} seconds", t.elapsed());
      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_});
      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_++;
      info(2,"postAdapt needed {} seconds", t.elapsed());
    }

    /** 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()
    {
      return Impl::loadBalanceImpl(*hostGrid_);
    }

    /*  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)
    {
      return Impl::loadBalanceImpl(*hostGrid_, handle);
    }

    /// Returns the grid change index, see \ref changeIndex.
    unsigned long changeIndex() const
    {
      return changeIndex_;
    }

  private:
    template <class Grid_>
    static std::shared_ptr<Self>
    instanceImpl(std::shared_ptr<Grid_> grid)
    {
      using mutex_type = std::shared_timed_mutex;
      static mutex_type access_mutex;

      // 1. Cleanup. Since the iterators may be invalidated only one accessor is allowed. Erases all
      //    expired weak pointers (all pointers that no longer have valid instances attached).
      std::unique_lock<mutex_type> write_lock(access_mutex);
      for (auto it = adaptiveGrids_.begin(); it != adaptiveGrids_.end(); ++it)
        if (it->expired())
          it = adaptiveGrids_.erase(it);
      write_lock.unlock();

      std::shared_lock<mutex_type> read_lock(access_mutex);
      // 2. Find matching observed grid object. We obtain a lock here to avoid complications with
      //    insertions made by other threads in step 3.
      auto it = find_if(adaptiveGrids_.begin(), adaptiveGrids_.end(),
        [&](std::weak_ptr<Self>& wPtr) {
          auto ag = wPtr.lock();
          return ag->hostGrid_ == grid;
        });

      // 3. Register new object or return existing. In case of inserting a new instance we obtain a
      //    write lock.
      if (it == adaptiveGrids_.end())
      {
        test_exit(!std::is_const<Grid_>::value,
          "No existing AdaptiveGrid found and no mutable grid passed to create a new instance");
        auto ptr = std::make_shared<Self>(std::move(grid), HiddenStruct{});
        read_lock.unlock();
        write_lock.lock();
        adaptiveGrids_.emplace_back(ptr);
        return std::move(ptr);
      }
      else
      {
        return it->lock();
      }
    }

    // Do-nothing specialization if argument is already an AdaptiveGrid
    static std::shared_ptr<Self> instanceImpl(std::shared_ptr<Self> grid) { return grid; }

  public:
    /** Returns the AdaptiveGrid associated to the grid passed.
     *  If no AdaptiveGrid exists, this method creates a new one if the passed grid is mutable,
     *  otherwise the call fails with an error.
     */
    template <class Grid_>
    static std::shared_ptr<Self> instance(Grid_&& grid)
    {
      return instanceImpl(wrap_or_share(FWD(grid)));
    }

    // Test for equality of the grid pointers
    bool operator==(AdaptiveGrid<HostGrid> const& that) const
    {
      return (hostGrid_ == that.hostGrid_);
    }

    /// Return the underlying grid
    std::shared_ptr<HostGrid> hostGrid() const { return hostGrid_; }

  private:
    /// List of previously created instances, indexed by address of the HostGrid.
    /**
     *  For each grid type Grid we maintain a list of instances handed out by the instance method.
     *  We use weak pointers here that are valid as long as there is at least one other place where
     *  the shared pointer to the instance is used. When the instance is no longer used the weak
     *  pointers here do not hinder the object's destruction.
     */
    static std::list<std::weak_ptr<Self>> adaptiveGrids_;

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

    /// Flag set during \ref preAdapt(), indicating whether any element might be coarsened in \ref adapt()
    bool mightCoarsen_ = false;

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

    /// This index is incremented every time the grid is changed, e.g. by refinement or coarsening.
    unsigned long changeIndex_ = 0;
  };

  template <class Grid>
  std::list<std::weak_ptr<AdaptiveGrid<Grid>>> AdaptiveGrid<Grid>::adaptiveGrids_;

} // end namespace AMDiS