DUNE PDELab (2.7)

directleaflocalordering.hh
1// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=8 sw=2 sts=2:
3
4#ifndef DUNE_PDELAB_ORDERING_DIRECTLEAFLOCALORDERING_HH
5#define DUNE_PDELAB_ORDERING_DIRECTLEAFLOCALORDERING_HH
6
7#include <dune/typetree/leafnode.hh>
8
9#include <dune/geometry/referenceelements.hh>
10#include <dune/localfunctions/common/interfaceswitch.hh>
11#include <dune/localfunctions/common/localkey.hh>
12#include <dune/pdelab/ordering/utility.hh>
13#include <dune/pdelab/gridfunctionspace/gridfunctionspacebase.hh>
14
15#include <vector>
16#include <numeric>
17
18namespace Dune {
19 namespace PDELab {
20
23
24 template<typename OrderingTag, typename FEM, typename ES, typename DI, typename CI>
25 class DirectLeafLocalOrdering
26 : public TypeTree::LeafNode
27 {
28
29 template<typename>
30 friend class LeafGridViewOrdering;
31
32 template<typename>
33 friend class LeafOrderingBase;
34
35 template<typename size_type>
36 friend struct ::Dune::PDELab::impl::update_ordering_data;
37
38 public:
39
40 typedef LocalOrderingTraits<ES,DI,CI> Traits;
41
42 private:
43
44 typedef impl::GridFunctionSpaceOrderingData<typename Traits::SizeType> GFSData;
45
46 public:
47
48 void map_local_index(const typename Traits::SizeType geometry_type_index,
49 const typename Traits::SizeType entity_index,
50 typename Traits::TreeIndexView mi,
51 typename Traits::ContainerIndex& ci) const
52 {
53 DUNE_THROW(NotImplemented,"not implemented");
54 }
55
56 template<typename ItIn, typename ItOut>
57 void map_lfs_indices(const ItIn begin, const ItIn end, ItOut out) const
58 {
59 // don't do anything - this is handled by the specialized GridViewOrdering
60 }
61
62 template<typename CIOutIterator, typename DIOutIterator = DummyDOFIndexIterator>
63 typename Traits::SizeType
64 extract_entity_indices(const typename Traits::DOFIndex::EntityIndex& ei,
65 typename Traits::SizeType child_index,
66 CIOutIterator ci_out, const CIOutIterator ci_end,
67 DIOutIterator di_out = DIOutIterator()) const
68 {
69 const typename Traits::SizeType s = size(ei);
70
71 // Handle DOF indices
72 for (typename Traits::SizeType i = 0; i < s; ++i, ++di_out)
73 di_out->treeIndex().push_back(i);
74
75 // only return the size, as the tree visitor expects that from all leaf nodes.
76 // The actual index processing is done by the specialized GridViewOrdering.
77 return s;
78 }
79
80 typename Traits::SizeType size(const typename Traits::DOFIndex::EntityIndex& index) const
81 {
82 return size(
83 Traits::DOFIndexAccessor::GeometryIndex::geometryType(index),
84 Traits::DOFIndexAccessor::GeometryIndex::entityIndex(index)
85 );
86 }
87
88 typename Traits::SizeType size(const typename Traits::SizeType geometry_type_index, const typename Traits::SizeType entity_index) const
89 {
90 typedef typename Traits::SizeType size_type;
91 if (_fixed_size)
92 return _gt_dof_sizes[geometry_type_index];
93 else if (_gt_used[geometry_type_index])
94 {
95 const size_type index = _gt_entity_offsets[geometry_type_index] + entity_index;
96 return _entity_dof_offsets[index+1] - _entity_dof_offsets[index];
97 }
98 else
99 return 0;
100 }
101
102 typename Traits::SizeType size(const typename Traits::SizeType geometry_type_index, const typename Traits::SizeType entity_index, const typename Traits::SizeType child_index) const
103 {
104 DUNE_THROW(NotImplemented,"not implemented");
105 }
106
107 typename Traits::SizeType offset(const typename Traits::SizeType geometry_type_index, const typename Traits::SizeType entity_index, const typename Traits::SizeType child_index) const
108 {
109 assert(child_index == 0);
110 return 0;
111 }
112
113 DirectLeafLocalOrdering(const std::shared_ptr<const FEM>& fem, const ES& es)
114 : _fem(fem)
115 , _es(es)
116 , _fixed_size(false)
117 , _container_blocked(false)
118 , _gfs_data(nullptr)
119 {}
120
121 const typename Traits::EntitySet& entitySet() const
122 {
123 return _es;
124 }
125
126 const FEM& finiteElementMap() const
127 {
128 return *_fem;
129 }
130
131 private:
132
133 static constexpr auto GT_UNUSED = ~std::size_t(0);
134
135 typedef FiniteElementInterfaceSwitch<
136 typename FEM::Traits::FiniteElement
137 > FESwitch;
138
139
140 void update_a_priori_fixed_size()
141 {
142 _fixed_size = _fem->fixedSize();
143 }
144
145 template<typename CodimMask>
146 void collect_used_codims(CodimMask& codims) const
147 {
148 for (typename ES::dim_type codim = 0; codim <= ES::dimension; ++codim)
149 if (_fem->hasDOFs(codim))
150 codims.set(codim);
151 }
152
153 template<typename It>
154 void update_fixed_size(It it, const It end)
155 {
156 assert(_fixed_size);
157
158 _max_local_size = _fem->maxLocalSize();
159
160 typedef typename Traits::SizeType size_type;
161 const size_type dim = Traits::GridView::dimension;
162 _codim_used.reset();
163 _gt_used.assign(GlobalGeometryTypeIndex::size(dim),false);
164 _gt_dof_sizes.assign(GlobalGeometryTypeIndex::size(dim),0);
165 for (; it != end; ++it)
166 {
167 size_type size = _fem->size(*it);
168 _gt_dof_sizes[GlobalGeometryTypeIndex::index(*it)] = size;
169 _gt_used[GlobalGeometryTypeIndex::index(*it)] = size > 0;
170 _codim_used[dim - it->dim()] = _codim_used[dim - it->dim()] || (size > 0);
171 }
172
173 _codim_fixed_size.set();
174 }
175
176
177 void pre_collect_used_geometry_types_from_cell()
178 {
179 typedef typename Traits::SizeType size_type;
180 const size_type dim = Traits::GridView::dimension;
181
182 _codim_used.reset();
183 _gt_used.assign(GlobalGeometryTypeIndex::size(dim),false);
184 _gt_dof_sizes.assign(GlobalGeometryTypeIndex::size(dim),GT_UNUSED);
185 _local_gt_dof_sizes.resize(GlobalGeometryTypeIndex::size(dim));
186 _max_local_size = 0;
187 _fixed_size_possible = true;
188 }
189
190
191 void collect_used_geometry_types_from_cell(const typename Traits::GridView::template Codim<0>::Entity& cell)
192 {
193 FESwitch::setStore(_fe_store,_fem->find(cell));
194
195 const typename FESwitch::Coefficients& coeffs =
196 FESwitch::coefficients(*_fe_store);
197
198 _max_local_size = std::max(_max_local_size,coeffs.size());
199
200 auto ref_el = ReferenceElements<typename Traits::GridView::ctype,Traits::GridView::dimension>::general(cell.type());
201
202 for (std::size_t i = 0; i < coeffs.size(); ++i)
203 {
204 const LocalKey& key = coeffs.localKey(i);
205 GeometryType gt = ref_el.type(key.subEntity(),key.codim());
206 _gt_used[GlobalGeometryTypeIndex::index(gt)] = true;
207 _codim_used.set(key.codim());
208 }
209 }
210
211
212 template<typename It>
213 void allocate_entity_offset_vector(It it, const It end)
214 {
215 _gt_entity_offsets.assign(GlobalGeometryTypeIndex::size(ES::dimension) + 1,0);
216 for (; it != end; ++it)
217 {
218 if (_gt_used[GlobalGeometryTypeIndex::index(*it)])
219 _gt_entity_offsets[GlobalGeometryTypeIndex::index(*it) + 1] = _es.indexSet().size(*it);
220 }
221 std::partial_sum(_gt_entity_offsets.begin(),_gt_entity_offsets.end(),_gt_entity_offsets.begin());
222 _entity_dof_offsets.assign(_gt_entity_offsets.back() + 1,0);
223
224 // Don't claim fixed size for any codim for now
225 _codim_fixed_size.reset();
226 }
227
228
229 void extract_per_entity_sizes_from_cell(const typename Traits::GridView::template Codim<0>::Entity& cell)
230 {
231 if (this->_fixed_size_possible)
232 std::fill(_local_gt_dof_sizes.begin(),_local_gt_dof_sizes.end(),0);
233
234 FESwitch::setStore(_fe_store,_fem->find(cell));
235
236 const typename FESwitch::Coefficients& coeffs =
237 FESwitch::coefficients(*_fe_store);
238
239 typedef typename Traits::SizeType size_type;
240
241 auto ref_el = ReferenceElements<typename Traits::GridView::ctype,Traits::GridView::dimension>::general(cell.type());
242
243 for (std::size_t i = 0; i < coeffs.size(); ++i)
244 {
245 const LocalKey& key = coeffs.localKey(i);
246 GeometryType gt = ref_el.type(key.subEntity(),key.codim());
247 const size_type geometry_type_index = GlobalGeometryTypeIndex::index(gt);
248
249 const size_type entity_index = _es.indexSet().subIndex(cell,key.subEntity(),key.codim());
250 const size_type index = _gt_entity_offsets[geometry_type_index] + entity_index;
251 _local_gt_dof_sizes[geometry_type_index] = _entity_dof_offsets[index+1] = std::max(_entity_dof_offsets[index+1],static_cast<size_type>(key.index() + 1));
252 }
253
254 if (_fixed_size_possible)
255 {
256 for (size_type i = 0; i < _local_gt_dof_sizes.size(); ++i)
257 {
258 if (_gt_dof_sizes[i] == GT_UNUSED)
259 _gt_dof_sizes[i] = _local_gt_dof_sizes[i];
260 else if (_gt_dof_sizes[i] != _local_gt_dof_sizes[i])
261 {
262 _fixed_size_possible = false;
263 break;
264 }
265 }
266 }
267 }
268
269
270 void finalize_non_fixed_size_update()
271 {
272 if (_fixed_size_possible)
273 {
274 // set size of unused geometry types to 0
275 for (auto& size : _gt_dof_sizes)
276 if (size == GT_UNUSED)
277 size = 0;
278 // free per-entity offsets
279 _entity_dof_offsets = std::vector<typename Traits::SizeType>();
280 _fixed_size = true;
281 _codim_fixed_size.set();
282 }
283 else
284 {
285 // convert per-entity sizes to offsets
286 std::partial_sum(_entity_dof_offsets.begin(),_entity_dof_offsets.end(),_entity_dof_offsets.begin());
287 _fixed_size = false;
288 _codim_fixed_size.reset();
289 }
290 }
291
292
293 typename Traits::SizeType maxLocalSize() const
294 {
295 return _max_local_size;
296 }
297
298 private:
299
300 bool update_gfs_data_size(typename Traits::SizeType& size, typename Traits::SizeType& block_count) const
301 {
302 return false;
303 }
304
305 protected:
306
307 std::shared_ptr<const FEM> _fem;
308 typename FESwitch::Store _fe_store;
309
310 ES _es;
311 bool _fixed_size;
312 bool _fixed_size_possible;
313 typename Traits::SizeType _max_local_size;
314 const bool _container_blocked;
315
316 typename Traits::CodimFlag _codim_used;
317 typename Traits::CodimFlag _codim_fixed_size;
318 std::vector<bool> _gt_used;
319
320 std::vector<typename Traits::SizeType> _gt_entity_offsets;
321 std::vector<typename Traits::SizeType> _gt_dof_sizes;
322 std::vector<typename Traits::SizeType> _entity_dof_offsets;
323 std::vector<typename Traits::SizeType> _local_gt_dof_sizes;
324
325 // This is only here to make the visitor happy that traverses all
326 // Orderings to manipulate the contained GFSData
327 GFSData* _gfs_data;
328
329 };
330
332
333 } // namespace PDELab
334} // namespace Dune
335
336#endif // DUNE_PDELAB_ORDERING_DIRECTLEAFLOCALORDERING_HH
Dune::GlobalGeometryTypeIndex::index
static constexpr std::size_t index(const GeometryType &gt)
Compute the index for the given geometry type over all dimensions.
Definition: typeindex.hh:133
Dune::GlobalGeometryTypeIndex::size
static constexpr std::size_t size(std::size_t maxdim)
Compute total number of geometry types up to and including the given dimension.
Definition: typeindex.hh:120
Dune::VTK::GeometryType
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:180
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::FloatCmp::gt
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
Dune::Simd::Overloads::max
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune
Dune namespace.
Definition: alignedallocator.hh:14
std
STL namespace.
Dune::FiniteElementInterfaceSwitch::coefficients
static const Coefficients & coefficients(const FiniteElement &fe)
access coefficients
Definition: interfaceswitch.hh:42
Dune::FiniteElementInterfaceSwitch::Coefficients
FiniteElement::Traits::Coefficients Coefficients
export the type of the coefficients
Definition: interfaceswitch.hh:33
Dune::FiniteElementInterfaceSwitch::setStore
static void setStore(Store &store, const FiniteElement &fe)
Store a finite element in the store.
Definition: interfaceswitch.hh:75
Dune::FiniteElementInterfaceSwitch::Store
std::shared_ptr< const FiniteElement > Store
Type for storing finite elements.
Definition: interfaceswitch.hh:68
Dune::Geo::ReferenceElements::general
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)