File: SubGridPart.hpp

package info (click to toggle)
opm-grid 2024.10%2Bds-3
  • links: PTS, VCS
  • area: main
  • in suites: trixie
  • size: 3,060 kB
  • sloc: cpp: 33,052; ansic: 3,210; sh: 67; makefile: 18
file content (381 lines) | stat: -rw-r--r-- 12,894 bytes parent folder | download | duplicates (2)
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
308
309
310
311
312
313
314
315
316
317
318
319
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
369
370
371
372
373
374
375
376
377
378
379
380
381
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=2 sts=2:

// Note: this file is based on defaultgridview.hh from Dune, and is therefore
// licensed under the Dune license (GPLv2 + runtime exception),
// see https://dune-project.org/about/license/
// rather than the OPM license (GPLv3+)
// Copyright 2021 Dune contributors.
// Copyright 2021 SINTEF Digital, Mathematics and Cybernetics.

#ifndef OPM_SUBGRIDPART_HEADER
#define OPM_SUBGRIDPART_HEADER

#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>

#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/gridview.hh>

#include <algorithm>
#include <cassert>
#include <cstddef>
#include <map>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <vector>

namespace Dune
{

template <class GridImp>
class SubGridPart;

template <class GridImp>
struct SubGridPartTraits {
    using GridPartImp = SubGridPart<GridImp>;

    /** \brief type of the grid */
    using Grid = typename std::remove_const<GridImp>::type;

    /** \brief type of the index set */
    using IndexSet = typename Grid ::Traits ::LeafIndexSet;

    /** \brief type of the intersection */
    using Intersection = typename Grid ::Traits ::LeafIntersection;

    /** \brief type of the intersection iterator */
    using IntersectionIterator = typename Grid ::Traits ::LeafIntersectionIterator;

    /** \brief type of the collective communication */
    using CollectiveCommunication = typename Grid ::Traits ::Communication;


    template <class BaseEntityType>
    class SubEntity : public BaseEntityType
    {
    public:
        SubEntity()
            : BaseEntityType()
            , is_owned_(false)
        {
        }
        SubEntity(const BaseEntityType& base, const bool owned)
            : BaseEntityType(base)
            , is_owned_(owned)
        {
        }
        auto partitionType() const
        {
            if (is_owned_) {
                return Dune::InteriorEntity;
            } else {
                return Dune::OverlapEntity;
            }
        }
    private:
        bool is_owned_;
    };

    template <int cd>
    struct Codim {
        using BaseEntity = typename Grid ::Traits ::template Codim<cd>::Entity;
        using Entity = SubEntity<BaseEntity>;
        using EntitySeed = typename Grid ::Traits ::template Codim<cd>::EntitySeed;

        using Geometry = typename Grid ::template Codim<cd>::Geometry;
        using LocalGeometry = typename Grid ::template Codim<cd>::LocalGeometry;

        /** \brief Define types needed to iterate over entities of a given partition type */
        template <PartitionIteratorType pit>
        struct Partition {
            /** \brief iterator over a given codim and partition type */
            using BaseIterator = typename Grid ::template Codim<cd>::template Partition<pit>::LeafIterator;
        };
    };

    enum { conforming = Capabilities ::isLeafwiseConforming<Grid>::v };
};


/// \brief A class to represent a part of a grid, similar to a GridView.
///
/// The differences from a GridView are:
///   - The SubGridPart consists of a set of elements (codim 0 entities),
///     considered to be Interior, and their neighbours, considered to be Overlap.
///   - When iterating over intersections on the elements, and accessing the outside()
///     elements, this can give you access to grid entities that are not in the
///     SubGridPart itself. This can only happen for intersections in the Overlap part.
///     For intersections of elements in the Interior part, the outside() element will
///     be either Interior or Overlap.
template <class GridImp>
class SubGridPart
{
    using ThisType = SubGridPart<GridImp>;

public:
    using Traits = SubGridPartTraits<GridImp>;

    /** \brief type of the grid */
    using Grid = typename Traits::Grid;

    /** \brief type of the index set */
    using IndexSet = typename Traits ::IndexSet;

    /** \brief type of the intersection */
    using Intersection = typename Traits ::Intersection;

    /** \brief type of the intersection iterator */
    using IntersectionIterator = typename Traits ::IntersectionIterator;

    /** \brief type of the collective communication */
    using CollectiveCommunication = typename Traits ::CollectiveCommunication;


    /** \brief Codim Structure */
    template <int cd>
    struct Codim : public Traits ::template Codim<cd> {
        using Entity = typename Traits::template Codim<cd>::Entity;
        class SubIterator
        {
        public:
            SubIterator(const SubGridPart& view, std::size_t index)
                : view_(&view)
                , index_(index)
            {
            }
            const Entity& operator*() const
            {
                entity_ = this->view_->get(index_);
                return entity_;
            }
            const Entity* operator->() const
            {
                entity_ = this->view_->get(index_);
                return &entity_;
            }
            SubIterator operator++()
            {
                ++index_;
                return *this;
            }
            SubIterator operator++(int)
            {
                Iterator copy(*this);
                ++index_;
                return copy;
            }
            bool operator==(const SubIterator& other) const
            {
                assert(view_ == other.view_);
                return index_ == other.index_;
            }
            bool operator!=(const SubIterator& other) const
            {
                assert(view_ == other.view_);
                return index_ != other.index_;
            }
        private:
            const SubGridPart* view_;
            std::size_t index_;
            mutable Entity entity_; // This may be low-performing for grids with large Entity objects.
        };
        using Iterator = SubIterator;

        /** \brief Define types needed to iterate over entities of a given partition type */
        template <PartitionIteratorType pit>
        struct Partition {
            /** \brief iterator over a given codim and partition type */
            using Iterator = SubIterator;
        };

    };

    enum { conforming = Traits::conforming };
    enum { dimension = GridImp::dimension };

public:

    /// Construct a view of the codim 0 entities that can be constructed from the seeds input.
    ///
    /// The seeds input is moved from and will be in a valid but indeterminate state after the call.
    SubGridPart(const Grid& grid,
                std::vector<typename Codim<0>::Entity::EntitySeed>&& seeds,
                const bool overlap = true)
        : grid_(&grid)
        , subset_(std::move(seeds))
        , num_owned_(subset_.size())
    {
        // Nothing more to do if we do not want to have overlap entities.
        if (!overlap) {
            return;
        }

        // Add neighbouring not-owned entities to subset_
        using Seed = typename Codim<0>::Entity::EntitySeed;
        std::unordered_set<int> owned;
        std::unordered_map<int, Seed> neighbors;
        const auto& iset = grid_->leafIndexSet();
        const auto& leaf_view = grid_->leafGridView();
        for (const auto& seed : subset_) {
            // Add this entity to the set of owned indices.
            const auto& entity = grid_->entity(seed);
            owned.insert(iset.index(entity));
            // Iterating over all intersections, ...
            const auto end = leaf_view.iend(entity);
            for (auto it = leaf_view.ibegin(entity); it != end; ++it) {
                if (it->boundary()) {
                    continue;
                }
                if (it->neighbor()) {
                    const auto outside_entity = it->outside();
                    // ...for all neighbour entities, add to neighbors.
                    neighbors.try_emplace(iset.index(outside_entity), outside_entity.seed());
                }
            }
        }
        // Now that owned is complete, we can eliminate any owned entries.
        std::map<int, Seed> unowned_neighbors;
        for (const auto& nb : neighbors) {
            if (owned.count(nb.first) == 0) {
                unowned_neighbors.insert(nb);
            }
        }
        subset_.resize(subset_.size() + unowned_neighbors.size());
        std::size_t count = num_owned_;
        for (const auto& neighbor : unowned_neighbors) {
            subset_[count] = neighbor.second;
            ++count;
        }
        assert(count == subset_.size());
    }

    /** \brief obtain a const reference to the underlying hierarchic grid */
    const Grid& grid() const
    {
        assert(grid_);
        return *grid_;
    }

    /** \brief obtain the index set */
    // const IndexSet& indexSet() const // Not implemented

    /** \brief obtain number of entities in a given codimension */
    int size(int codim) const
    {
        if (codim == 0) {
            return subset_.size();
        } else {
            return 0;
        }
    }

    /** \brief obtain number of entities with a given geometry type */
    // int size(const GeometryType& type) const // Not implemented

    /** \brief obtain begin iterator for this view */
    template <int cd>
    typename Codim<cd>::Iterator begin() const
    {
        static_assert(cd == 0, "Only codimension 0 iterators for SubGridPart.");
        using Iterator = typename Codim<cd>::Iterator;
        return Iterator(*this, 0);
    }

    /** \brief obtain end iterator for this view */
    template <int cd>
    typename Codim<cd>::Iterator end() const
    {
        static_assert(cd == 0, "Only codimension 0 iterators for SubGridPart.");
        using Iterator = typename Codim<cd>::Iterator;
        return Iterator(*this, subset_.size());
    }

    // We support iterating over Interior_Partition, Overlap_Partition and All_Partition

    /** \brief obtain begin iterator for this view */
    template <int cd, PartitionIteratorType pit>
    typename Codim<cd>::template Partition<pit>::Iterator begin() const
    {
        static_assert(cd == 0, "Only codimension 0 iterators for SubGridPart.");
        static_assert(pit == Interior_Partition || pit == Overlap_Partition || pit == All_Partition);
        if constexpr (pit == Interior_Partition || pit == All_Partition) {
            return begin<0>();
        } else {
            // Overlap partition starts at index num_owned_.
            // Note that it may be empty, i.e. begin() == end().
            return typename Codim<cd>::Iterator(*this, num_owned_);
        }
    }

    /** \brief obtain end iterator for this view */
    template <int cd, PartitionIteratorType pit>
    typename Codim<cd>::template Partition<pit>::Iterator end() const
    {
        static_assert(cd == 0, "Only codimension 0 iterators for SubGridPart.");
        static_assert(pit == Interior_Partition || pit == Overlap_Partition || pit == All_Partition);
        if constexpr (pit == Overlap_Partition || pit == All_Partition) {
            return end<0>();
        } else {
            // Interior partition ends before index num_owned_.
            return typename Codim<cd>::Iterator(*this, num_owned_);
        }
    }

    /** \brief obtain begin intersection iterator with respect to this view */
    IntersectionIterator ibegin(const typename Codim<0>::Entity& entity) const
    {
        return entity.impl().ileafbegin();
    }

    /** \brief obtain end intersection iterator with respect to this view */
    IntersectionIterator iend(const typename Codim<0>::Entity& entity) const
    {
        return entity.impl().ileafend();
    }

    /** \brief obtain collective communication object */
    const CollectiveCommunication& comm() const
    {
        return grid().comm();
    }

    /** \brief Return size of the overlap region for a given codim on the grid view.  */
    int overlapSize(int codim) const
    {
        if (codim == 0) {
            return subset_.size() - num_owned_;
        } else {
            return 0;
        }
    }

    /** \brief Return size of the ghost region for a given codim on the grid view.  */
    int ghostSize([[maybe_unused]] int codim) const
    {
        return 0;
    }

    /** communicate data on this view */
    // template <class DataHandleImp, class DataType>
    // void
    // communicate(CommDataHandleIF<DataHandleImp, DataType>& data, InterfaceType iftype, CommunicationDirection dir) const
    // Not implemented.

private:
    using Entity0 = typename Codim<0>::Entity;
    Entity0 get(std::size_t ii) const
    {
        const bool owned = ii < num_owned_;
        return Entity0(grid_->entity(subset_[ii]), owned);
    }
    const Grid* grid_;
    std::vector<typename Entity0::EntitySeed> subset_;
    const std::size_t num_owned_;
};

} // namespace Dune

#endif // OPM_SUBGRIDPART_HEADER