libMesh
coupling_matrix.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 #ifndef LIBMESH_COUPLING_MATRIX_H
21 #define LIBMESH_COUPLING_MATRIX_H
22 
23 // Local Includes
24 #include "libmesh/libmesh_common.h"
25 
26 // C++ includes
27 #include <algorithm>
28 #include <limits>
29 #include <utility> // std::pair
30 #include <vector>
31 
32 namespace libMesh
33 {
34 
35 // Forward declarations
36 class ConstCouplingAccessor;
37 class CouplingAccessor;
38 
39 class ConstCouplingRow;
40 
41 class ConstCouplingRowConstIterator;
42 
55 {
56 public:
57 
61  explicit
62  CouplingMatrix (const unsigned int n=0);
63 
67  bool operator() (const unsigned int i,
68  const unsigned int j) const;
69 
74  CouplingAccessor operator() (const unsigned int i,
75  const unsigned int j);
76 
81  unsigned int size() const;
82 
87  void resize(const unsigned int n);
88 
92  void clear();
93 
97  bool empty() const;
98 
99  CouplingMatrix & operator&= (const CouplingMatrix & other);
100 
101 private:
102 
103  friend class ConstCouplingAccessor;
104  friend class CouplingAccessor;
105  friend class ConstCouplingRow;
107 
119  typedef std::pair<std::size_t, std::size_t> range_type;
120  typedef std::vector<range_type> rc_type;
121  rc_type _ranges;
122 
126  unsigned int _size;
127 };
128 
129 
130 
135 {
136 public:
137  ConstCouplingAccessor(std::size_t loc_in,
138  const CouplingMatrix & mat_in) :
139  _location(loc_in), _mat(mat_in)
140  {
141  libmesh_assert_less(_location, _mat.size() * _mat.size());
142  }
143 
144  operator bool() const {
145  const std::size_t max_size = std::numeric_limits<std::size_t>::max();
146 
147  // Find the range that might contain i,j
148  // lower_bound isn't *quite* what we want
149  CouplingMatrix::rc_type::const_iterator lb = std::upper_bound
150  (_mat._ranges.begin(), _mat._ranges.end(),
151  std::make_pair(_location, max_size));
152  if (lb!=_mat._ranges.begin())
153  --lb;
154  else
155  lb=_mat._ranges.end();
156 
157  // If no range might contain i,j then it's 0
158  if (lb == _mat._ranges.end())
159  return false;
160 
161  const std::size_t lastloc = lb->second;
162 
163 #ifdef DEBUG
164  const std::size_t firstloc = lb->first;
165  libmesh_assert_less_equal(firstloc, lastloc);
166  libmesh_assert_less_equal(firstloc, _location);
167 
168  CouplingMatrix::rc_type::const_iterator next = lb;
169  next++;
170  if (next != _mat._ranges.end())
171  {
172  // Ranges should be sorted and should not touch
173  libmesh_assert_greater(next->first, lastloc+1);
174  }
175 #endif
176 
177  return (lastloc >= _location);
178  }
179 
180 protected:
181 
182  std::size_t _location;
184 };
185 
186 
191 {
192 public:
193  CouplingAccessor(std::size_t loc_in,
194  CouplingMatrix & mat_in) :
195  ConstCouplingAccessor(loc_in, mat_in), _my_mat(mat_in) {}
196 
197  template <typename T>
199  {
200  // For backwards compatibility we take integer arguments,
201  // but coupling matrix entries are really all zero or one.
202  const bool as_bool = new_value;
203  libmesh_assert_equal_to(new_value, as_bool);
204 
205  *this = as_bool;
206  return *this;
207  }
208 
209  CouplingAccessor & operator = (bool new_value)
210  {
211  const std::size_t max_size = std::numeric_limits<std::size_t>::max();
212 
213  // Find the range that might contain i,j
214  // lower_bound isn't *quite* what we want
215  CouplingMatrix::rc_type::iterator lb =
216  std::upper_bound (_my_mat._ranges.begin(), _my_mat._ranges.end(),
217  std::make_pair(_location, max_size));
218  if (lb!=_my_mat._ranges.begin())
219  --lb;
220  else
221  lb=_my_mat._ranges.end();
222 
223  // If no range might contain i,j then we might need to make a new
224  // one.
225  if (lb == _my_mat._ranges.end())
226  {
227  if (new_value == true)
228  _my_mat._ranges.insert(_my_mat._ranges.begin(),
229  std::make_pair(_location, _location));
230  }
231  else
232  {
233  const std::size_t firstloc = lb->first;
234  const std::size_t lastloc = lb->second;
235  libmesh_assert_less_equal(firstloc, lastloc);
236  libmesh_assert_less_equal(firstloc, _location);
237 
238 #ifdef DEBUG
239  CouplingMatrix::rc_type::const_iterator next = lb;
240  next++;
241  if (next != _my_mat._ranges.end())
242  {
243  // Ranges should be sorted and should not touch
244  libmesh_assert_greater(next->first, lastloc+1);
245  }
246 #endif
247 
248  // If we're in this range, we might need to shorten or remove
249  // or split it
250  if (new_value == false)
251  {
252  if (_location == firstloc)
253  {
254  if (_location == lastloc)
255  {
256  _my_mat._ranges.erase(lb);
257  }
258  else
259  {
260  libmesh_assert_less (lb->first, lastloc);
261  lb->first++;
262  }
263  }
264  else if (_location == lastloc)
265  {
266  libmesh_assert_less (firstloc, lb->second);
267 
268  lb->second--;
269  }
270  else if (_location < lastloc)
271  {
272  libmesh_assert_less_equal(_location+1, lastloc);
273 
274  lb->first = _location+1;
275 
276  libmesh_assert_less_equal(firstloc, _location-1);
277 
278  _my_mat._ranges.insert
279  (lb, std::make_pair(firstloc, _location-1));
280  }
281  }
282 
283  // If we're not in this range, we might need to extend it or
284  // join it with its neighbor or add a new one.
285  else // new_value == true
286  {
287  CouplingMatrix::rc_type::iterator next = lb;
288  next++;
289  const std::size_t nextloc =
290  (next == _my_mat._ranges.end()) ?
292  next->first;
293 
294  // Ranges should be sorted and should not touch
295  libmesh_assert_greater(nextloc, lastloc+1);
296 
297  if (_location > lastloc)
298  {
299  if (_location == lastloc + 1)
300  {
301  if (_location == nextloc - 1)
302  {
303  next->first = firstloc;
304  _my_mat._ranges.erase(lb);
305  }
306  else
307  lb->second++;
308  }
309  else
310  {
311  if (_location == nextloc - 1)
312  next->first--;
313  else
314  _my_mat._ranges.insert
315  (next, std::make_pair(_location, _location));
316  }
317  }
318  }
319  }
320 
321  return *this;
322  }
323 
324 private:
325 
327 };
328 
329 
330 
336 {
337 public:
338  ConstCouplingRow(unsigned int row_in,
339  const CouplingMatrix & mat_in) :
340  _row_i(row_in), _mat(mat_in)
341  {
342  libmesh_assert_less(_row_i, _mat.size());
343 
344  // Location for i,N
345  _begin_location = _row_i*_mat.size()+_mat.size()-1;
346 
347  const std::size_t max_size = std::numeric_limits<std::size_t>::max();
348 
349  // Find the range that might contain i,N
350  // lower_bound isn't *quite* what we want
351  _begin_it = std::upper_bound
352  (_mat._ranges.begin(), _mat._ranges.end(),
353  std::make_pair(_begin_location, max_size));
354  if (_begin_it !=_mat._ranges.begin())
355  --_begin_it;
356  else
357  _begin_it=_mat._ranges.end();
358 
359  // If that range doesn't exist then we're an empty row
360  if (_begin_it == _mat._ranges.end())
361  _begin_location = max_size;
362  else
363  {
364  const std::size_t lastloc = _begin_it->second;
365 #ifdef DEBUG
366  const std::size_t firstloc = _begin_it->first;
367  libmesh_assert_less_equal(firstloc, lastloc);
368 #endif
369 
370  // If that range ends before i,0 then we're an empty row
371  std::size_t zero_location = _row_i*_mat.size();
372  if (zero_location > lastloc)
373  {
374  _begin_location = max_size;
375  _begin_it = _mat._ranges.end();
376  }
377  else
378  // We have *some* entry(s) in this row, we just need to find
379  // the earliest
380  {
381  while (_begin_it != _mat._ranges.begin())
382  {
383  CouplingMatrix::rc_type::const_iterator prev =
384  _begin_it;
385  --prev;
386 
387  if (prev->second < zero_location)
388  break;
389 
390  _begin_it = prev;
391  }
392  if (_begin_it->first < zero_location)
393  _begin_location = zero_location;
394  else
395  _begin_location = _begin_it->first;
396  }
397  }
398  }
399 
400  /*
401  * A forward iterator type for looping over indices in this row
402  */
404 
405  /*
406  * An iterator to the first index in this row, or to end() for an
407  * empty row
408  */
409  const_iterator begin() const;
410 
411  /*
412  * An iterator representing past-the-end of this row
413  */
414  const_iterator end() const;
415 
416  bool operator== (const ConstCouplingRow & other) const
417  {
418  // Thinking that rows from different matrix objects are equal is
419  // not even wrong
420  libmesh_assert(&_mat == &other._mat);
421 
422  return ((_begin_location == other._begin_location) &&
423  (_begin_it == other._begin_it));
424  }
425 
426  bool operator!= (const ConstCouplingRow & other) const
427  {
428  return !(*this == other);
429  }
430 protected:
431 
433 
434  unsigned int _row_i;
436 
437  // The location (i*size+j) corresponding to the first entry in this
438  // row, or numeric_limits<size_t>::max() for an empty row.
439  std::size_t _begin_location;
440 
441  // Iterator to the range containing the first row element, or
442  // _row._mat._ranges.end() if no CouplingMatrix values are true for
443  // this row
444  CouplingMatrix::rc_type::const_iterator _begin_it;
445 };
446 
447 
448 
450 {
451 public:
453  std::size_t loc_in,
454  CouplingMatrix::rc_type::const_iterator it_in) :
455  _location(loc_in),
456  _row(row_in),
457  _it(it_in)
458  {
459 #ifndef NDEBUG
460  if (_it != _row._mat._ranges.end())
461  {
462  libmesh_assert_less_equal(_it->first, _location);
463  libmesh_assert_less_equal(_location, _it->second);
464  }
465  else
466  {
467  libmesh_assert_equal_to
468  (_location, std::numeric_limits<size_t>::max());
469  }
470 #endif
471  }
472 
473  unsigned int operator* ()
474  {
475  libmesh_assert_not_equal_to
477  return _location % _row._mat.size();
478  }
479 
481  {
482  libmesh_assert_not_equal_to
484 
485  if (_location == _it->second)
486  {
487  ++_it;
488 
489  // Are we past the end of the matrix?
490  if (_it == _row._mat._ranges.end())
492  else
493  {
494  _location = _it->first;
495  // Are we past the end of the row?
496  if (_location >= _row._mat.size()*(_row._row_i+1))
497  {
499  _it = _row._mat._ranges.end();
500  }
501  }
502  }
503  else
504  ++_location;
505 
506  return *this;
507  }
508 
509  bool operator== (const ConstCouplingRowConstIterator & other) const
510  {
511  // Thinking that iterators from different row objects are equal
512  // is not even wrong
513  libmesh_assert(_row == other._row);
514 
515  return ((_location == other._location) &&
516  (_it == other._it));
517  }
518 
519  bool operator!= (const ConstCouplingRowConstIterator & other) const
520  {
521  return !(*this == other);
522  }
523 
524 private:
525  // The location (i*size+j) corresponding to this iterator, or
526  // numeric_limits<size_t>::max() to signify end()
527  std::size_t _location;
529  // The range containing this iterator location, or
530  // _row._mat._ranges.end() if no range contains this location
531  CouplingMatrix::rc_type::const_iterator _it;
532 };
533 
534 
535 
536 //--------------------------------------------------
537 // ConstCouplingRow inline methods
538 inline
540  return const_iterator (*this, _begin_location, _begin_it);
541 }
542 
543 inline
545  return const_iterator
547  _mat._ranges.end());
548 }
549 
550 
551 
552 //--------------------------------------------------
553 // CouplingMatrix inline methods
554 inline
555 CouplingMatrix::CouplingMatrix (const unsigned int n) :
556  _ranges(), _size(n)
557 {
558  this->resize(n);
559 }
560 
561 
562 
563 inline
564 bool CouplingMatrix::operator() (const unsigned int i,
565  const unsigned int j) const
566 {
567  libmesh_assert_less (i, _size);
568  libmesh_assert_less (j, _size);
569 
570  const std::size_t location = std::size_t(i)*_size + j;
571 
572  return bool(ConstCouplingAccessor(location, *this));
573 }
574 
575 
576 
577 
578 inline
580  const unsigned int j)
581 {
582  const std::size_t location = std::size_t(i)*_size + j;
583 
584  return CouplingAccessor(location, *this);
585 }
586 
587 
588 
589 inline
590 unsigned int CouplingMatrix::size() const
591 {
592  return _size;
593 }
594 
595 
596 
597 inline
598 void CouplingMatrix::resize(const unsigned int n)
599 {
600  _size = n;
601 
602  _ranges.clear();
603 }
604 
605 
606 
607 inline
609 {
610  this->resize(0);
611 }
612 
613 
614 
615 inline
617 {
618  return (_size == 0);
619 }
620 
621 
622 } // namespace libMesh
623 
624 
625 #endif // LIBMESH_COUPLING_MATRIX_H
void clear()
Clears the matrix.
ConstCouplingAccessor(std::size_t loc_in, const CouplingMatrix &mat_in)
boostcopy::enable_if_c< ScalarTraits< Scalar >::value, TypeTensor< typename CompareTypes< T, Scalar >::supertype > >::type operator*(const Scalar factor, const TypeTensor< T > &t)
Definition: type_tensor.h:929
bool operator()(const unsigned int i, const unsigned int j) const
static const unsigned int prev[3]
A lookup table for the decrement modulo 3 operation, for iterating through the three nodes per elemen...
This proxy class acts like a container of indices from a single coupling row.
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
const_iterator end() const
The libMesh namespace provides an interface to certain functionality in the library.
CouplingMatrix(const unsigned int n=0)
Constructor.
bool operator!=(const OrderWrapper &lhs, const OrderWrapper &rhs)
Definition: fe_type.h:95
long double max(long double a, double b)
libmesh_assert(j)
This accessor class allows simple setting of CouplingMatrix values.
CouplingMatrix & operator&=(const CouplingMatrix &other)
unsigned int _size
The size of the matrix.
ConstCouplingRowConstIterator(const ConstCouplingRow &row_in, std::size_t loc_in, CouplingMatrix::rc_type::const_iterator it_in)
bool operator==(const OrderWrapper &lhs, const OrderWrapper &rhs)
Overload comparison operators for OrderWrapper.
Definition: fe_type.h:94
friend class ConstCouplingAccessor
const_iterator begin() const
Iterator & operator++()
op++() forwards on to Iter::op++()
static const unsigned int next[3]
A lookup table for the increment modulo 3 operation, for iterating through the three nodes per elemen...
void resize(const unsigned int n)
Resizes the matrix and initializes all entries to be 0.
Iterator & operator=(const Iterator &rhs)
Assignment operator.
std::pair< std::size_t, std::size_t > range_type
Coupling matrices are typically either full or very sparse, and all values are only zero or one...
std::vector< range_type > rc_type
CouplingMatrix::rc_type::const_iterator _it
This accessor class allows simple access to CouplingMatrix values.
ConstCouplingRow(unsigned int row_in, const CouplingMatrix &mat_in)
unsigned int size() const
ConstCouplingRowConstIterator const_iterator
const CouplingMatrix & _mat
const CouplingMatrix & _mat
CouplingMatrix::rc_type::const_iterator _begin_it
CouplingAccessor(std::size_t loc_in, CouplingMatrix &mat_in)
This class defines a coupling matrix.