dune-pdelab  2.0.0
borderdofexchanger.hh
Go to the documentation of this file.
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 <cstddef>
7 #include <vector>
8 #include <algorithm>
9 
10 #include <dune/common/deprecated.hh>
11 #include <dune/common/parallel/mpihelper.hh>
12 #include <dune/common/tuples.hh>
13 
14 #include <dune/geometry/typeindex.hh>
15 #include <dune/grid/common/gridenums.hh>
16 #include <dune/grid/common/datahandleif.hh>
17 
23 
24 namespace Dune {
25  namespace PDELab {
26 
27 
31 
33 
65  template<typename GridOperator>
67  {
68  typedef typename GridOperator::Traits::Jacobian M;
69  typedef M Matrix;
71  typedef typename GridOperatorTraits::JacobianField Scalar;
72  typedef typename GridOperatorTraits::TrialGridFunctionSpace GFSU;
73  typedef typename GridOperatorTraits::TestGridFunctionSpace GFSV;
74  typedef typename GFSV::Traits::GridViewType GridView;
75  static const int dim = GridView::dimension;
76  typedef typename GridView::Traits::Grid Grid;
77  typedef typename Matrix::block_type BlockType;
78  typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
79  typedef typename Grid::Traits::GlobalIdSet IdSet;
80  typedef typename IdSet::IdType IdType;
81 
84  typename GFSV::Ordering::Traits::DOFIndex::value_type,
85  GFSV::Ordering::Traits::DOFIndex::max_depth,
86  typename GFSV::Traits::GridView::Grid::GlobalIdSet::IdType
88 
89  public:
91  typedef unordered_map<
92  typename GFSV::Ordering::Traits::DOFIndex,
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:
117  : _communication_cache(make_shared<CommunicationCache>(grid_operator))
118  , _grid_view(grid_operator.testGridFunctionSpace().gridView())
119  {}
120 
121  void update(const GridOperator& grid_operator)
122  {
123  _communication_cache = make_shared<CommunicationCache>(grid_operator);
124  }
125 
127  : public BorderIndexIdCache<GFSV>
128  {
129 
132 
133  public:
134 
136  : BaseT(go.testGridFunctionSpace())
137  , _gfsu(go.trialGridFunctionSpace())
138  , _initialized(false)
139  , _entity_cache(go.testGridFunctionSpace())
140  {}
141 
142  typedef IdType EntityID;
143  typedef typename GFSU::Ordering::Traits::DOFIndex::TreeIndex ColumnTreeIndex;
144  typedef std::size_t size_type;
145 
146  bool initialized() const
147  {
148  return _initialized;
149  }
150 
152  {
153  _initialized = true;
154  }
155 
156  void update()
157  {
158  BaseT::update();
159  _border_pattern.clear();
160  _initialized = false;
161  }
162 
163 
164  const BorderPattern& pattern() const
165  {
166  assert(initialized());
167  return _border_pattern;
168  }
169 
170  template<typename LFSVCache, typename LFSUCache, typename LocalPattern>
171  void addEntries(const LFSVCache& lfsv_cache, const LFSUCache& lfsu_cache, const LocalPattern& pattern)
172  {
173  assert(!initialized());
174 
175  for (typename LocalPattern::const_iterator it = pattern.begin(),
176  end_it = pattern.end();
177  it != end_it;
178  ++it)
179  {
180  // skip constrained entries for now. TODO: Is this correct??
181  if (lfsv_cache.isConstrained(it->i()) || lfsu_cache.isConstrained(it->j()))
182  continue;
183 
184  const typename LFSVCache::DOFIndex& di = lfsv_cache.dofIndex(it->i());
185  const typename LFSUCache::DOFIndex& dj = lfsu_cache.dofIndex(it->j());
186 
187  size_type row_gt_index = GFSV::Ordering::Traits::DOFIndexAccessor::geometryType(di);
188  size_type row_entity_index = GFSV::Ordering::Traits::DOFIndexAccessor::entityIndex(di);
189 
190  size_type col_gt_index = GFSU::Ordering::Traits::DOFIndexAccessor::geometryType(dj);
191  size_type col_entity_index = GFSU::Ordering::Traits::DOFIndexAccessor::entityIndex(dj);
192 
193  // We are only interested in connections between two border entities.
194  if (!this->isBorderEntity(row_gt_index,row_entity_index) ||
195  !this->isBorderEntity(col_gt_index,col_entity_index))
196  continue;
197 
198  _border_pattern[di].insert(GlobalDOFIndex(this->id(col_gt_index,col_entity_index),dj.treeIndex()));
199  }
200  }
201 
202 
203  template<typename Entity>
204  size_type size(const Entity& e) const
205  {
206  _entity_cache.update(e);
207  size_type n = 0;
208  for (size_type i = 0; i < _entity_cache.size(); ++i)
209  {
210  typename BorderPattern::const_iterator it = _border_pattern.find(_entity_cache.dofIndex(i));
211  if (!transfer_dof(i,it))
212  continue;
213  n += it->second.size();
214  }
215 
216  return n;
217  }
218 
219  template<typename Buffer, typename Entity>
220  void gather_pattern(Buffer& buf, const Entity& e) const
221  {
222  _entity_cache.update(e);
223  for (size_type i = 0; i < _entity_cache.size(); ++i)
224  {
225  typename BorderPattern::const_iterator it = _border_pattern.find(_entity_cache.dofIndex(i));
226  if (!transfer_dof(i,it))
227  continue;
228  for (typename BorderPattern::mapped_type::const_iterator col_it = it->second.begin(),
229  col_end = it->second.end();
230  col_it != col_end;
231  ++col_it)
232  buf.write(std::make_pair(_entity_cache.dofIndex(i).treeIndex(),*col_it));
233  }
234  }
235 
236  template<typename Buffer, typename Entity>
237  void gather_data(Buffer& buf, const Entity& e, const M& matrix) const
238  {
239  _entity_cache.update(e);
240  for (size_type i = 0; i < _entity_cache.size(); ++i)
241  {
242  typename BorderPattern::const_iterator it = _border_pattern.find(_entity_cache.dofIndex(i));
243  if (!transfer_dof(i,it))
244  continue;
245  for (typename BorderPattern::mapped_type::const_iterator col_it = it->second.begin(),
246  col_end = it->second.end();
247  col_it != col_end;
248  ++col_it)
249  {
250  typename BaseT::EntityIndex col_entity = this->index(col_it->entityID());
251 
252  ColDOFIndex dj;
253  GFSU::Ordering::Traits::DOFIndexAccessor::store(dj,col_entity.geometryTypeIndex(),col_entity.entityIndex(),col_it->treeIndex());
254  buf.write(make_tuple(_entity_cache.dofIndex(i).treeIndex(),*col_it,matrix(_entity_cache.containerIndex(i),_gfsu.ordering().mapIndex(dj))));
255  }
256  }
257  }
258 
259  private:
260 
261  bool transfer_dof(size_type i, typename BorderPattern::const_iterator it) const
262  {
263  // not a border DOF
264  if (it == _border_pattern.end())
265  return false;
266  else
267  return true;
268 
269  /* Constraints check moved to addEntry()
270  // check for constraint
271  typename GridOperator::Traits::TrialGridFunctionSpaceConstraints::const_iterator cit = _constraints->find(_entity_cache.containerIndex(i));
272 
273  // transfer if DOF is not constrained
274  // TODO: What about non-Dirichlet constraints??
275  return cit == _constraints->end();
276  */
277  }
278 
279  const GFSU& _gfsu;
280  BorderPattern _border_pattern;
281  bool _initialized;
282  mutable EntityIndexCache<GFSV,true> _entity_cache;
283 
284  };
285 
287  template<typename Pattern>
289  : public CommDataHandleIF<PatternExtender<Pattern>, PatternMPIData>
290  {
291 
292  typedef std::size_t size_type;
293 
294  public:
296  typedef PatternMPIData DataType;
297 
298  bool contains (int dim, int codim) const
299  {
300  return
301  codim > 0 &&
302  (_gfsu.dataHandleContains(codim) ||
303  _gfsv.dataHandleContains(codim));
304  }
305 
306  bool fixedsize (int dim, int codim) const
307  {
308  return false;
309  }
310 
313  template<typename Entity>
314  size_type size (Entity& e) const
315  {
316  if (Entity::codimension == 0)
317  return 0;
318 
319  return _communication_cache.size(e);
320  }
321 
324  template<typename MessageBuffer, typename Entity>
325  void gather (MessageBuffer& buff, const Entity& e) const
326  {
327  if (Entity::codimension == 0)
328  return;
329 
330  _communication_cache.gather_pattern(buff,e);
331  }
332 
335  template<typename MessageBuffer, typename Entity>
336  void scatter (MessageBuffer& buff, const Entity& e, size_t n)
337  {
338  if (Entity::codimension == 0)
339  return;
340 
341  for (size_type i = 0; i < n; ++i)
342  {
343  DataType data;
344  buff.read(data);
345 
346  std::pair<bool,typename CommunicationCache::EntityIndex> col_index = _communication_cache.findIndex(data.second.entityID());
347  if (!col_index.first)
348  continue;
349 
350  RowDOFIndex di;
351  GFSV::Ordering::Traits::DOFIndexAccessor::store(di,
352  e.type(),
353  _grid_view.indexSet().index(e),
354  data.first);
355 
356  ColDOFIndex dj;
357  GFSU::Ordering::Traits::DOFIndexAccessor::store(dj,
358  col_index.second.geometryTypeIndex(),
359  col_index.second.entityIndex(),
360  data.second.treeIndex());
361 
362  _pattern.add_link(_gfsv.ordering().mapIndex(di),_gfsu.ordering().mapIndex(dj));
363  }
364  }
365 
367  const GFSU& gfsu,
368  const GFSV& gfsv,
369  Pattern& pattern)
370  : _communication_cache(dof_exchanger.communicationCache())
371  , _grid_view(dof_exchanger.gridView())
372  , _gfsu(gfsu)
373  , _gfsv(gfsv)
374  , _pattern(pattern)
375  {}
376 
377  private:
378 
379  const CommunicationCache& _communication_cache;
380  GridView _grid_view;
381  const GFSU& _gfsu;
382  const GFSV& _gfsv;
383  Pattern& _pattern;
384 
385  };
386 
389  : public CommDataHandleIF<EntryAccumulator,ValueMPIData>
390  {
391 
392  typedef std::size_t size_type;
393 
394  public:
396  typedef ValueMPIData DataType;
397 
398  bool contains(int dim, int codim) const
399  {
400  return
401  codim > 0 &&
402  (_gfsu.dataHandleContains(codim) ||
403  _gfsv.dataHandleContains(codim));
404  }
405 
406  bool fixedsize(int dim, int codim) const
407  {
408  return false;
409  }
410 
411  template<typename Entity>
412  size_type size(Entity& e) const
413  {
414  if (Entity::codimension == 0)
415  return 0;
416 
417  return _communication_cache.size(e);
418  }
419 
420  template<typename MessageBuffer, typename Entity>
421  void gather(MessageBuffer& buff, const Entity& e) const
422  {
423  if (Entity::codimension == 0)
424  return;
425 
426  _communication_cache.gather_data(buff,e,_matrix);
427  }
428 
431  template<typename MessageBuffer, typename Entity>
432  void scatter(MessageBuffer& buff, const Entity& e, size_type n)
433  {
434  if (Entity::codimension == 0)
435  return;
436 
437  for (size_type i = 0; i < n; ++i)
438  {
439  DataType data;
440  buff.read(data);
441 
442  std::pair<bool,typename CommunicationCache::EntityIndex> col_index = _communication_cache.findIndex(get<1>(data).entityID());
443  if (!col_index.first)
444  continue;
445 
446  RowDOFIndex di;
447  GFSV::Ordering::Traits::DOFIndexAccessor::store(di,
448  e.type(),
449  _grid_view.indexSet().index(e),
450  get<0>(data));
451 
452  ColDOFIndex dj;
453  GFSU::Ordering::Traits::DOFIndexAccessor::store(dj,
454  col_index.second.geometryTypeIndex(),
455  col_index.second.entityIndex(),
456  get<1>(data).treeIndex());
457 
458  _matrix(_gfsv.ordering().mapIndex(di),_gfsu.ordering().mapIndex(dj)) += get<2>(data);
459  }
460  }
461 
462 
464  const GFSU& gfsu,
465  const GFSV& gfsv,
466  Matrix& matrix)
467  : _communication_cache(dof_exchanger.communicationCache())
468  , _grid_view(dof_exchanger.gridView())
469  , _gfsu(gfsu)
470  , _gfsv(gfsv)
471  , _matrix(matrix)
472  {}
473 
474  private:
475 
476  const CommunicationCache& _communication_cache;
477  GridView _grid_view;
478  const GFSU& _gfsu;
479  const GFSV& _gfsv;
480  Matrix& _matrix;
481 
482  };
483 
487  void accumulateBorderEntries(const GridOperator& grid_operator, Matrix& matrix)
488  {
489  if (_grid_view.comm().size() > 1)
490  {
491  EntryAccumulator data_handle(*this,
492  grid_operator.testGridFunctionSpace(),
493  grid_operator.trialGridFunctionSpace(),
494  matrix);
495  _grid_view.communicate(data_handle,
496  InteriorBorder_InteriorBorder_Interface,
497  ForwardCommunication);
498  }
499  }
500 
502  {
503  return *_communication_cache;
504  }
505 
507  {
508  return *_communication_cache;
509  }
510 
511  shared_ptr<CommunicationCache> communicationCacheStorage()
512  {
513  return _communication_cache;
514  }
515 
516  const GridView& gridView() const
517  {
518  return _grid_view;
519  }
520 
521  private:
522 
523  shared_ptr<CommunicationCache> _communication_cache;
524  GridView _grid_view;
525 
526  };
527 
528 
529  template<typename GridOperator>
531  {
532 
533  public:
534 
536 
538  typedef Empty BorderPattern;
539 
541  {}
542 
544  {}
545 
546  void accumulateBorderEntries(const GridOperator& grid_operator, typename GridOperator::Traits::Jacobian& matrix)
547  {}
548 
550  {
551  return *this;
552  }
553 
555  {
556  return *this;
557  }
558 
559  void update(const GridOperator& grid_operator)
560  {}
561 
562  };
563 
564 
565  template<typename GridOperator>
567  public NoDataBorderDOFExchanger<GridOperator>
568  {
569 
570  public:
571 
573  {}
574 
576  {}
577 
578  };
579 
580 
581  } // namespace PDELab
582 } // namespace Dune
583 
584 #endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_BORDERDOFEXCHANGER_HH
void update(const Entity &e)
Definition: entityindexcache.hh:49
const DI & dofIndex(size_type i) const
Definition: entityindexcache.hh:58
CommunicationCache & communicationCache()
Definition: borderdofexchanger.hh:501
JF JacobianField
The field type of the jacobian.
Definition: gridoperatorutilities.hh:69
void update(const GridOperator &grid_operator)
Definition: borderdofexchanger.hh:121
NoDataBorderDOFExchanger CommunicationCache
Definition: borderdofexchanger.hh:535
A DataHandle class to exchange matrix sparsity patterns.
Definition: borderdofexchanger.hh:288
Provide common name for std::unordered_map and std::unordered_multimap classes in Dune::PDELab namesp...
size_type size(Entity &e) const
How many objects of type DataType have to be sent for a given entity.
Definition: borderdofexchanger.hh:314
GFSU TrialGridFunctionSpace
The trial grid function space.
Definition: gridoperatorutilities.hh:37
void update()
Definition: borderdofexchanger.hh:156
size_type size(const Entity &e) const
Definition: borderdofexchanger.hh:204
Definition: borderdofexchanger.hh:566
Helper class for adding up matrix entries on border.
Definition: borderdofexchanger.hh:66
Traits class for the grid operator.
Definition: gridoperatorutilities.hh:33
const CI & containerIndex(size_type i) const
Definition: entityindexcache.hh:64
void gather_pattern(Buffer &buf, const Entity &e) const
Definition: borderdofexchanger.hh:220
bool isBorderEntity(std::size_t gt_index, std::size_t entity_index) const
Definition: borderindexidcache.hh:114
PatternExtender(const NonOverlappingBorderDOFExchanger &dof_exchanger, const GFSU &gfsu, const GFSV &gfsv, Pattern &pattern)
Definition: borderdofexchanger.hh:366
const GFSU & trialGridFunctionSpace() const
Get the trial grid function space.
Definition: gridoperator.hh:87
size_type size(Entity &e) const
Definition: borderdofexchanger.hh:412
Provide common name for std::unordered_set and std::unordered_multiset classes in Dune::PDELab namesp...
Standard grid operator implementation.
Definition: gridoperator.hh:34
PatternMPIData DataType
Export type of data for message buffer.
Definition: borderdofexchanger.hh:296
const BorderPattern & pattern() const
Definition: borderdofexchanger.hh:164
void update(const GridOperator &grid_operator)
Definition: borderdofexchanger.hh:559
bool initialized() const
Definition: borderdofexchanger.hh:146
NonOverlappingBorderDOFExchanger(const GridOperator &grid_operator)
Constructor. Sets up the local to global relations.
Definition: borderdofexchanger.hh:116
void accumulateBorderEntries(const GridOperator &grid_operator, typename GridOperator::Traits::Jacobian &matrix)
Definition: borderdofexchanger.hh:546
size_type size() const
Definition: entityindexcache.hh:74
const GridView & gridView() const
Definition: borderdofexchanger.hh:516
A DataHandle class to exchange matrix entries.
Definition: borderdofexchanger.hh:388
void gather(MessageBuffer &buff, const Entity &e) const
Pack data from user to message buffer.
Definition: borderdofexchanger.hh:325
std::pair< bool, EntityIndex > findIndex(id_type entity_id) const
Definition: borderindexidcache.hh:139
NoDataBorderDOFExchanger()
Definition: borderdofexchanger.hh:540
OverlappingBorderDOFExchanger(const GridOperator &grid_operator)
Definition: borderdofexchanger.hh:575
void finishInitialization()
Definition: borderdofexchanger.hh:151
IdType EntityID
Definition: borderdofexchanger.hh:142
void gather(MessageBuffer &buff, const Entity &e) const
Definition: borderdofexchanger.hh:421
shared_ptr< CommunicationCache > communicationCacheStorage()
Definition: borderdofexchanger.hh:511
void accumulateBorderEntries(const GridOperator &grid_operator, Matrix &matrix)
Sums up the entries corresponding to border vertices.
Definition: borderdofexchanger.hh:487
void gather_data(Buffer &buf, const Entity &e, const M &matrix) const
Definition: borderdofexchanger.hh:237
Definition: globaldofindex.hh:14
unordered_map< typename GFSV::Ordering::Traits::DOFIndex, unordered_set< GlobalDOFIndex > > BorderPattern
Data structure for storing border-border matrix pattern entries in a communication-optimized form...
Definition: borderdofexchanger.hh:94
void update()
Definition: borderindexidcache.hh:88
const CommunicationCache & communicationCache() const
Definition: borderdofexchanger.hh:554
Empty BorderPattern
Data structure for storing border-border matrix pattern entries in a communication-optimized form...
Definition: borderdofexchanger.hh:538
EntryAccumulator(const NonOverlappingBorderDOFExchanger &dof_exchanger, const GFSU &gfsu, const GFSV &gfsv, Matrix &matrix)
Definition: borderdofexchanger.hh:463
bool contains(int dim, int codim) const
Definition: borderdofexchanger.hh:398
Definition: borderindexidcache.hh:25
Definition: unordered_map.hh:46
const GFSV & testGridFunctionSpace() const
Get the test grid function space.
Definition: gridoperator.hh:93
CommunicationCache & communicationCache()
Definition: borderdofexchanger.hh:549
void scatter(MessageBuffer &buff, const Entity &e, size_type n)
Unpack data from message buffer to user.
Definition: borderdofexchanger.hh:432
Definition: borderdofexchanger.hh:530
GFSU::Ordering::Traits::DOFIndex::TreeIndex ColumnTreeIndex
Definition: borderdofexchanger.hh:143
NoDataBorderDOFExchanger(const GridOperator &grid_operator)
Definition: borderdofexchanger.hh:543
bool contains(int dim, int codim) const
Definition: borderdofexchanger.hh:298
GFSV TestGridFunctionSpace
The test grid function space.
Definition: gridoperatorutilities.hh:40
void addEntries(const LFSVCache &lfsv_cache, const LFSUCache &lfsu_cache, const LocalPattern &pattern)
Definition: borderdofexchanger.hh:171
OverlappingBorderDOFExchanger()
Definition: borderdofexchanger.hh:572
const CommunicationCache & communicationCache() const
Definition: borderdofexchanger.hh:506
void scatter(MessageBuffer &buff, const Entity &e, size_t n)
Unpack data from message buffer to user.
Definition: borderdofexchanger.hh:336
CommunicationCache(const GridOperator &go)
Definition: borderdofexchanger.hh:135
Dune::PDELab::BackendMatrixSelector< MB, Domain, Range, JF >::Type Jacobian
The type of the jacobian.
Definition: gridoperatorutilities.hh:72
ValueMPIData DataType
Export type of data for message buffer.
Definition: borderdofexchanger.hh:396
EntityIndex index(id_type entity_id) const
Definition: borderindexidcache.hh:129
const E & e
Definition: interpolate.hh:172
Definition: unordered_set.hh:45
bool fixedsize(int dim, int codim) const
Definition: borderdofexchanger.hh:406
bool fixedsize(int dim, int codim) const
Definition: borderdofexchanger.hh:306
std::size_t size_type
Definition: borderdofexchanger.hh:144