DUNE PDELab (git)

borderdofexchanger.hh
1// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=8 sw=2 sts=2:
3#ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_BORDERDOFEXCHANGER_HH
4#define DUNE_PDELAB_GRIDOPERATOR_COMMON_BORDERDOFEXCHANGER_HH
5
6#include <algorithm>
7#include <cstddef>
8#include <tuple>
9#include <unordered_map>
10#include <unordered_set>
11#include <vector>
12
13#include <dune/common/deprecated.hh>
14#include <dune/common/parallel/mpihelper.hh>
15
16#include <dune/geometry/typeindex.hh>
17
18#include <dune/grid/common/datahandleif.hh>
19#include <dune/grid/common/gridenums.hh>
20
21#include <dune/pdelab/common/borderindexidcache.hh>
22#include <dune/pdelab/common/globaldofindex.hh>
23#include <dune/pdelab/gridfunctionspace/entityindexcache.hh>
24
25namespace Dune {
26 namespace PDELab {
27
28
32
34
66 template<typename GridOperator>
67 class NonOverlappingBorderDOFExchanger
68 {
69 typedef typename GridOperator::Traits::Jacobian M;
70 typedef M Matrix;
71 typedef typename GridOperator::Traits GridOperatorTraits;
72 typedef typename GridOperatorTraits::JacobianField Scalar;
73 typedef typename GridOperatorTraits::TrialGridFunctionSpace GFSU;
74 typedef typename GridOperatorTraits::TestGridFunctionSpace GFSV;
75 using EntitySet = typename GFSV::Traits::EntitySet;
76 static const int dim = EntitySet::dimension;
77 using Grid = typename EntitySet::Traits::GridView::Traits::Grid;
78 typedef typename Matrix::block_type BlockType;
79 typedef typename Grid::Traits::GlobalIdSet IdSet;
80 typedef typename IdSet::IdType IdType;
81
83 typedef Dune::PDELab::GlobalDOFIndex<
84 typename GFSV::Ordering::Traits::DOFIndex::value_type,
85 GFSV::Ordering::Traits::DOFIndex::max_depth,
86 typename IdSet::IdType
87 > GlobalDOFIndex;
88
89 public:
91 typedef std::unordered_map<
92 typename GFSV::Ordering::Traits::DOFIndex,
93 std::unordered_set<GlobalDOFIndex>
94 > BorderPattern;
95
96 private:
97 typedef typename GFSV::Ordering::Traits::DOFIndex RowDOFIndex;
98 typedef typename GFSU::Ordering::Traits::DOFIndex ColDOFIndex;
99
100 typedef std::pair<
101 typename RowDOFIndex::TreeIndex,
102 typename BorderPattern::mapped_type::value_type
103 > PatternMPIData;
104
105 typedef std::tuple<
106 typename RowDOFIndex::TreeIndex,
107 typename BorderPattern::mapped_type::value_type,
108 typename M::field_type
109 > ValueMPIData;
110
111 public:
118 NonOverlappingBorderDOFExchanger(const GridOperator& grid_operator)
119 : _communication_cache(std::make_shared<CommunicationCache>(grid_operator))
120 , _entity_set(grid_operator.testGridFunctionSpace().entitySet())
121 {}
122
123 void update(const GridOperator& grid_operator)
124 {
125 _communication_cache = std::make_shared<CommunicationCache>(grid_operator);
126 }
127
128 class CommunicationCache
129 : public BorderIndexIdCache<GFSV>
130 {
131
132 friend class NonOverlappingBorderDOFExchanger;
133 typedef BorderIndexIdCache<GFSV> BaseT;
134
135 public:
136
137 CommunicationCache(const GridOperator& go)
138 : BaseT(go.testGridFunctionSpace())
139 , _gfsu(go.trialGridFunctionSpace())
140 , _initialized(false)
141 , _entity_cache(go.testGridFunctionSpace())
142 {}
143
144 typedef IdType EntityID;
145 typedef typename GFSU::Ordering::Traits::DOFIndex::TreeIndex ColumnTreeIndex;
146 typedef std::size_t size_type;
147
148 bool initialized() const
149 {
150 return _initialized;
151 }
152
153 void finishInitialization()
154 {
155 _initialized = true;
156 }
157
158 void update()
159 {
160 BaseT::update();
161 _border_pattern.clear();
162 _initialized = false;
163 }
164
165
166 const BorderPattern& pattern() const
167 {
168 assert(initialized());
169 return _border_pattern;
170 }
171
172 template<typename LFSVCache, typename LFSUCache, typename LocalPattern>
173 void addEntries(const LFSVCache& lfsv_cache, const LFSUCache& lfsu_cache, const LocalPattern& pattern)
174 {
175 assert(!initialized());
176
177 for (typename LocalPattern::const_iterator it = pattern.begin(),
178 end_it = pattern.end();
179 it != end_it;
180 ++it)
181 {
182 // skip constrained entries for now. TODO: Is this correct??
183 if (lfsv_cache.isConstrained(it->i()) || lfsu_cache.isConstrained(it->j()))
184 continue;
185
186 const typename LFSVCache::DOFIndex& di = lfsv_cache.dofIndex(it->i());
187 const typename LFSUCache::DOFIndex& dj = lfsu_cache.dofIndex(it->j());
188
189 size_type row_gt_index = GFSV::Ordering::Traits::DOFIndexAccessor::geometryType(di);
190 size_type row_entity_index = GFSV::Ordering::Traits::DOFIndexAccessor::entityIndex(di);
191
192 size_type col_gt_index = GFSU::Ordering::Traits::DOFIndexAccessor::geometryType(dj);
193 size_type col_entity_index = GFSU::Ordering::Traits::DOFIndexAccessor::entityIndex(dj);
194
195 // We are only interested in connections between two border entities.
196 if (!this->isBorderEntity(row_gt_index,row_entity_index) ||
197 !this->isBorderEntity(col_gt_index,col_entity_index))
198 continue;
199
200 _border_pattern[di].insert(GlobalDOFIndex(this->id(col_gt_index,col_entity_index),dj.treeIndex()));
201 }
202 }
203
204
205 template<typename Entity>
206 size_type size(const Entity& e) const
207 {
208 if (!_gfsu.entitySet().contains(e))
209 return 0;
210 _entity_cache.update(e);
211 size_type n = 0;
212 for (size_type i = 0; i < _entity_cache.size(); ++i)
213 {
214 typename BorderPattern::const_iterator it = _border_pattern.find(_entity_cache.dofIndex(i));
215 if (!transfer_dof(i,it))
216 continue;
217 n += it->second.size();
218 }
219
220 return n;
221 }
222
223 template<typename Buffer, typename Entity>
224 void gather_pattern(Buffer& buf, const Entity& e) const
225 {
226 if (!_gfsu.entitySet().contains(e))
227 return;
228 _entity_cache.update(e);
229 for (size_type i = 0; i < _entity_cache.size(); ++i)
230 {
231 typename BorderPattern::const_iterator it = _border_pattern.find(_entity_cache.dofIndex(i));
232 if (!transfer_dof(i,it))
233 continue;
234 for (typename BorderPattern::mapped_type::const_iterator col_it = it->second.begin(),
235 col_end = it->second.end();
236 col_it != col_end;
237 ++col_it)
238 buf.write(std::make_pair(_entity_cache.dofIndex(i).treeIndex(),*col_it));
239 }
240 }
241
242 template<typename Buffer, typename Entity>
243 void gather_data(Buffer& buf, const Entity& e, const M& matrix) const
244 {
245 if (!_gfsu.entitySet().contains(e))
246 return;
247 _entity_cache.update(e);
248 for (size_type i = 0; i < _entity_cache.size(); ++i)
249 {
250 typename BorderPattern::const_iterator it = _border_pattern.find(_entity_cache.dofIndex(i));
251 if (!transfer_dof(i,it))
252 continue;
253 for (typename BorderPattern::mapped_type::const_iterator col_it = it->second.begin(),
254 col_end = it->second.end();
255 col_it != col_end;
256 ++col_it)
257 {
258 typename BaseT::EntityIndex col_entity = this->index(col_it->entityID());
259
260 ColDOFIndex dj;
261 GFSU::Ordering::Traits::DOFIndexAccessor::store(dj,col_entity.geometryTypeIndex(),col_entity.entityIndex(),col_it->treeIndex());
262 buf.write(std::make_tuple(_entity_cache.dofIndex(i).treeIndex(),*col_it,matrix(_entity_cache.containerIndex(i),_gfsu.ordering().mapIndex(dj))));
263 }
264 }
265 }
266
267 private:
268
269 bool transfer_dof(size_type i, typename BorderPattern::const_iterator it) const
270 {
271 // not a border DOF
272 if (it == _border_pattern.end())
273 return false;
274 else
275 return true;
276
277 /* Constraints check moved to addEntry()
278 // check for constraint
279 typename GridOperator::Traits::TrialGridFunctionSpaceConstraints::const_iterator cit = _constraints->find(_entity_cache.containerIndex(i));
280
281 // transfer if DOF is not constrained
282 // TODO: What about non-Dirichlet constraints??
283 return cit == _constraints->end();
284 */
285 }
286
287 const GFSU& _gfsu;
288 BorderPattern _border_pattern;
289 bool _initialized;
290 mutable EntityIndexCache<GFSV,true> _entity_cache;
291
292 };
293
295 template<typename Pattern>
296 class PatternExtender
297 : public CommDataHandleIF<PatternExtender<Pattern>, PatternMPIData>
298 {
299
300 typedef std::size_t size_type;
301
302 public:
304 typedef PatternMPIData DataType;
305
306 bool contains (int dim, int codim) const
307 {
308 return
309 codim > 0 &&
310 (_gfsu.dataHandleContains(codim) ||
311 _gfsv.dataHandleContains(codim));
312 }
313
314 bool fixedSize (int dim, int codim) const
315 {
316 return false;
317 }
318
321 template<typename Entity>
322 size_type size (Entity& e) const
323 {
324 if (Entity::codimension == 0)
325 return 0;
326
327 return _communication_cache.size(e);
328 }
329
332 template<typename MessageBuffer, typename Entity>
333 void gather (MessageBuffer& buff, const Entity& e) const
334 {
335 if (Entity::codimension == 0)
336 return;
337
338 _communication_cache.gather_pattern(buff,e);
339 }
340
343 template<typename MessageBuffer, typename Entity>
344 void scatter (MessageBuffer& buff, const Entity& e, size_t n)
345 {
346 if (Entity::codimension == 0)
347 return;
348
349 for (size_type i = 0; i < n; ++i)
350 {
351 DataType data;
352 buff.read(data);
353
354 std::pair<bool,typename CommunicationCache::EntityIndex> col_index = _communication_cache.findIndex(data.second.entityID());
355 if (!col_index.first)
356 continue;
357
358 RowDOFIndex di;
359 GFSV::Ordering::Traits::DOFIndexAccessor::store(di,
360 e.type(),
361 _entity_set.indexSet().index(e),
362 data.first);
363
364 ColDOFIndex dj;
365 GFSU::Ordering::Traits::DOFIndexAccessor::store(dj,
366 col_index.second.geometryTypeIndex(),
367 col_index.second.entityIndex(),
368 data.second.treeIndex());
369
370 _pattern.add_link(_gfsv.ordering().mapIndex(di),_gfsu.ordering().mapIndex(dj));
371 }
372 }
373
374 PatternExtender(const NonOverlappingBorderDOFExchanger& dof_exchanger,
375 const GFSU& gfsu,
376 const GFSV& gfsv,
377 Pattern& pattern)
378 : _communication_cache(dof_exchanger.communicationCache())
379 , _entity_set(dof_exchanger.entitySet())
380 , _gfsu(gfsu)
381 , _gfsv(gfsv)
382 , _pattern(pattern)
383 {}
384
385 private:
386
387 const CommunicationCache& _communication_cache;
388 const EntitySet _entity_set;
389 const GFSU& _gfsu;
390 const GFSV& _gfsv;
391 Pattern& _pattern;
392
393 };
394
396 class EntryAccumulator
397 : public CommDataHandleIF<EntryAccumulator,ValueMPIData>
398 {
399
400 typedef std::size_t size_type;
401
402 public:
404 typedef ValueMPIData DataType;
405
406 bool contains(int dim, int codim) const
407 {
408 return
409 codim > 0 &&
410 (_gfsu.dataHandleContains(codim) ||
411 _gfsv.dataHandleContains(codim));
412 }
413
414 bool fixedSize(int dim, int codim) const
415 {
416 return false;
417 }
418
419 template<typename Entity>
420 size_type size(Entity& e) const
421 {
422 if (Entity::codimension == 0)
423 return 0;
424
425 return _communication_cache.size(e);
426 }
427
428 template<typename MessageBuffer, typename Entity>
429 void gather(MessageBuffer& buff, const Entity& e) const
430 {
431 if (Entity::codimension == 0)
432 return;
433
434 _communication_cache.gather_data(buff,e,_matrix);
435 }
436
439 template<typename MessageBuffer, typename Entity>
440 void scatter(MessageBuffer& buff, const Entity& e, size_type n)
441 {
442 if (Entity::codimension == 0)
443 return;
444
445 for (size_type i = 0; i < n; ++i)
446 {
447 DataType data;
448 buff.read(data);
449
450 std::pair<bool,typename CommunicationCache::EntityIndex> col_index = _communication_cache.findIndex(std::get<1>(data).entityID());
451 if (!col_index.first)
452 continue;
453
454 RowDOFIndex di;
455 GFSV::Ordering::Traits::DOFIndexAccessor::store(di,
456 e.type(),
457 _entity_set.indexSet().index(e),
458 std::get<0>(data));
459
460 ColDOFIndex dj;
461 GFSU::Ordering::Traits::DOFIndexAccessor::store(dj,
462 col_index.second.geometryTypeIndex(),
463 col_index.second.entityIndex(),
464 std::get<1>(data).treeIndex());
465
466 _matrix(_gfsv.ordering().mapIndex(di),_gfsu.ordering().mapIndex(dj)) += std::get<2>(data);
467 }
468 }
469
470
471 EntryAccumulator(const NonOverlappingBorderDOFExchanger& dof_exchanger,
472 const GFSU& gfsu,
473 const GFSV& gfsv,
474 Matrix& matrix)
475 : _communication_cache(dof_exchanger.communicationCache())
476 , _entity_set(dof_exchanger.entitySet())
477 , _gfsu(gfsu)
478 , _gfsv(gfsv)
479 , _matrix(matrix)
480 {}
481
482 private:
483
484 const CommunicationCache& _communication_cache;
485 EntitySet _entity_set;
486 const GFSU& _gfsu;
487 const GFSV& _gfsv;
488 Matrix& _matrix;
489
490 };
491
498 void accumulateBorderEntries(const GridOperator& grid_operator, Matrix& matrix)
499 {
500 if (_entity_set.gridView().comm().size() > 1)
501 {
502 EntryAccumulator data_handle(*this,
503 grid_operator.testGridFunctionSpace(),
504 grid_operator.trialGridFunctionSpace(),
505 matrix);
506 _entity_set.gridView().communicate(data_handle,
507 InteriorBorder_InteriorBorder_Interface,
508 ForwardCommunication);
509 }
510 }
511
512 CommunicationCache& communicationCache()
513 {
514 return *_communication_cache;
515 }
516
517 const CommunicationCache& communicationCache() const
518 {
519 return *_communication_cache;
520 }
521
522 std::shared_ptr<CommunicationCache> communicationCacheStorage()
523 {
524 return _communication_cache;
525 }
526
527 const EntitySet& entitySet() const
528 {
529 return _entity_set;
530 }
531
532 private:
533
534 std::shared_ptr<CommunicationCache> _communication_cache;
535 EntitySet _entity_set;
536
537 };
538
539
540 template<typename GridOperator>
541 class NoDataBorderDOFExchanger
542 {
543
544 public:
545
546 typedef NoDataBorderDOFExchanger CommunicationCache;
547
549 typedef Empty BorderPattern;
550
551 NoDataBorderDOFExchanger()
552 {}
553
554 NoDataBorderDOFExchanger(const GridOperator& grid_operator)
555 {}
556
557 void accumulateBorderEntries(const GridOperator& grid_operator, typename GridOperator::Traits::Jacobian& matrix)
558 {}
559
560 CommunicationCache& communicationCache()
561 {
562 return *this;
563 }
564
565 const CommunicationCache& communicationCache() const
566 {
567 return *this;
568 }
569
570 void update(const GridOperator& grid_operator)
571 {}
572
573 };
574
575
576 template<typename GridOperator>
577 class OverlappingBorderDOFExchanger :
578 public NoDataBorderDOFExchanger<GridOperator>
579 {
580
581 public:
582
583 OverlappingBorderDOFExchanger()
584 {}
585
586 OverlappingBorderDOFExchanger(const GridOperator& grid_operator)
587 {}
588
589 };
590
591
592 } // namespace PDELab
593} // namespace Dune
594
595#endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_BORDERDOFEXCHANGER_HH
Dune::CommDataHandleIF
CommDataHandleIF describes the features of a data handle for communication in parallel runs using the...
Definition: datahandleif.hh:78
Dune::Entity
Wrapper class for entities.
Definition: entity.hh:66
Dune::Entity::codimension
static constexpr int codimension
Know your own codimension.
Definition: entity.hh:106
Dune::IdSet::IdType
IdTypeImp IdType
Type used to represent an id.
Definition: indexidset.hh:453
Dune::Matrix
A generic dynamic dense matrix.
Definition: matrix.hh:561
Dune::Matrix::block_type
T block_type
Export the type representing the components.
Definition: matrix.hh:568
Dune::PDELab::GridOperator
Standard grid operator implementation.
Definition: gridoperator.hh:36
Dune::PDELab::GridOperator::trialGridFunctionSpace
const GFSU & trialGridFunctionSpace() const
Get the trial grid function space.
Definition: gridoperator.hh:92
Dune::PDELab::GridOperator::testGridFunctionSpace
const GFSV & testGridFunctionSpace() const
Get the test grid function space.
Definition: gridoperator.hh:98
Dune::PDELab::NonOverlappingBorderDOFExchanger::EntryAccumulator
A DataHandle class to exchange matrix entries.
Definition: borderdofexchanger.hh:398
Dune::PDELab::NonOverlappingBorderDOFExchanger::PatternExtender
A DataHandle class to exchange matrix sparsity patterns.
Definition: borderdofexchanger.hh:298
Dune::PDELab::NonOverlappingBorderDOFExchanger
Helper class for adding up matrix entries on border.
Definition: borderdofexchanger.hh:68
datahandleif.hh
Describes the parallel communication interface class for MessageBuffers and DataHandles.
deprecated.hh
Definition of the DUNE_NO_DEPRECATED_* macros.
Dune::ForwardCommunication
@ ForwardCommunication
communicate as given in InterfaceType
Definition: gridenums.hh:171
Dune::InteriorBorder_InteriorBorder_Interface
@ InteriorBorder_InteriorBorder_Interface
send/receive interior and border entities
Definition: gridenums.hh:87
Dune::PDELab::NonOverlappingBorderDOFExchanger::accumulateBorderEntries
void accumulateBorderEntries(const GridOperator &grid_operator, Matrix &matrix)
Sums up the entries corresponding to border vertices.
Definition: borderdofexchanger.hh:498
Dune::PDELab::NonOverlappingBorderDOFExchanger::PatternExtender::size
size_type size(Entity &e) const
How many objects of type DataType have to be sent for a given entity.
Definition: borderdofexchanger.hh:322
Dune::PDELab::NonOverlappingBorderDOFExchanger::NonOverlappingBorderDOFExchanger
NonOverlappingBorderDOFExchanger(const GridOperator &grid_operator)
Constructor. Sets up the local to global relations.
Definition: borderdofexchanger.hh:118
Dune::PDELab::NonOverlappingBorderDOFExchanger::BorderPattern
std::unordered_map< typename GFSV::Ordering::Traits::DOFIndex, std::unordered_set< GlobalDOFIndex > > BorderPattern
Data structure for storing border-border matrix pattern entries in a communication-optimized form.
Definition: borderdofexchanger.hh:94
Dune::PDELab::NonOverlappingBorderDOFExchanger::PatternExtender::DataType
PatternMPIData DataType
Export type of data for message buffer.
Definition: borderdofexchanger.hh:304
Dune::PDELab::NonOverlappingBorderDOFExchanger::PatternExtender::gather
void gather(MessageBuffer &buff, const Entity &e) const
Pack data from user to message buffer.
Definition: borderdofexchanger.hh:333
Dune::PDELab::NonOverlappingBorderDOFExchanger::PatternExtender::scatter
void scatter(MessageBuffer &buff, const Entity &e, size_t n)
Unpack data from message buffer to user.
Definition: borderdofexchanger.hh:344
Dune::PDELab::NonOverlappingBorderDOFExchanger::EntryAccumulator::DataType
ValueMPIData DataType
Export type of data for message buffer.
Definition: borderdofexchanger.hh:404
Dune::PDELab::NonOverlappingBorderDOFExchanger::EntryAccumulator::scatter
void scatter(MessageBuffer &buff, const Entity &e, size_type n)
Unpack data from message buffer to user.
Definition: borderdofexchanger.hh:440
Dune::PDELab::NoDataBorderDOFExchanger::BorderPattern
Empty BorderPattern
Data structure for storing border-border matrix pattern entries in a communication-optimized form.
Definition: borderdofexchanger.hh:549
mpihelper.hh
Helpers for dealing with MPI.
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
std
STL namespace.
Dune::Empty
Just an empty class.
Definition: typetraits.hh:55
Dune::PDELab::GridOperatorTraits
Traits class for the grid operator.
Definition: gridoperatorutilities.hh:34
Dune::PDELab::GridOperatorTraits::JacobianField
JF JacobianField
The field type of the jacobian.
Definition: gridoperatorutilities.hh:69
Dune::PDELab::GridOperatorTraits::TrialGridFunctionSpace
GFSU TrialGridFunctionSpace
The trial grid function space.
Definition: gridoperatorutilities.hh:37
Dune::PDELab::GridOperatorTraits::TestGridFunctionSpace
GFSV TestGridFunctionSpace
The test grid function space.
Definition: gridoperatorutilities.hh:40
Dune::PDELab::GridOperatorTraits::Jacobian
Dune::PDELab::Backend::Matrix< MB, Domain, Range, JF > Jacobian
The type of the jacobian.
Definition: gridoperatorutilities.hh:72
typeindex.hh
Helper classes to provide indices for geometrytypes for use in a vector.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)