Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
43 : /**
44 : * This class defines a coupling matrix. A coupling
45 : * matrix is simply a matrix of ones and zeros describing
46 : * how different components in a system couple with each
47 : * other. A coupling matrix is necessarily square but not
48 : * necessarily symmetric.
49 : *
50 : * \author Benjamin S. Kirk
51 : * \date 2002
52 : * \brief Defines the coupling between variables of a System.
53 : */
54 0 : class CouplingMatrix
55 : {
56 : public:
57 :
58 : /**
59 : * Constructor.
60 : */
61 : explicit
62 : CouplingMatrix (const std::size_t n=0);
63 :
64 : /**
65 : * \returns The (i,j) entry of the matrix.
66 : */
67 : bool operator() (const std::size_t i,
68 : const std::size_t j) const;
69 :
70 : /**
71 : * \returns The (i,j) entry of the matrix as
72 : * a smart-reference.
73 : */
74 : CouplingAccessor operator() (const std::size_t i,
75 : const std::size_t j);
76 :
77 : /**
78 : * \returns The size of the matrix, i.e. N for an
79 : * NxN matrix.
80 : */
81 : std::size_t size() const;
82 :
83 : /**
84 : * Resizes the matrix and initializes
85 : * all entries to be 0.
86 : */
87 : void resize(const std::size_t n);
88 :
89 : /**
90 : * Clears the matrix.
91 : */
92 : void clear();
93 :
94 : /**
95 : * \returns \p true if the matrix is empty.
96 : */
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;
106 : friend class ConstCouplingRowConstIterator;
107 :
108 : /**
109 : * Coupling matrices are typically either full or very sparse, and
110 : * all values are only zero or one.
111 : *
112 : * We store non-zeros as ranges: the first entry of each range pair
113 : * is the location of the first non-zero, and the second is the
114 : * location of the last subsequent non-zero (*not* the next
115 : * subsequent zero; we drop empty ranges).
116 : *
117 : * We store locations (i,j) as long integers i*_size+j
118 : */
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 :
123 : /**
124 : * The size of the matrix.
125 : */
126 : std::size_t _size;
127 : };
128 :
129 :
130 :
131 : /**
132 : * This accessor class allows simple access to CouplingMatrix values.
133 : */
134 : class ConstCouplingAccessor
135 : {
136 : public:
137 0 : ConstCouplingAccessor(std::size_t loc_in,
138 0 : const CouplingMatrix & mat_in) :
139 0 : _location(loc_in), _mat(mat_in)
140 : {
141 0 : libmesh_assert_less(_location, _mat.size() * _mat.size());
142 0 : }
143 :
144 : operator bool() const
145 : {
146 : const std::size_t max_size = std::numeric_limits<std::size_t>::max();
147 :
148 : // Find the range that might contain i,j
149 : // lower_bound isn't *quite* what we want
150 : CouplingMatrix::rc_type::const_iterator lb = std::upper_bound
151 : (_mat._ranges.begin(), _mat._ranges.end(),
152 : std::make_pair(_location, max_size));
153 : if (lb!=_mat._ranges.begin())
154 : --lb;
155 : else
156 : lb=_mat._ranges.end();
157 :
158 : // If no range might contain i,j then it's 0
159 : if (lb == _mat._ranges.end())
160 : return false;
161 :
162 : const std::size_t lastloc = lb->second;
163 :
164 : #ifdef DEBUG
165 : const std::size_t firstloc = lb->first;
166 : libmesh_assert_less_equal(firstloc, lastloc);
167 : libmesh_assert_less_equal(firstloc, _location);
168 :
169 : CouplingMatrix::rc_type::const_iterator next = lb;
170 : next++;
171 : if (next != _mat._ranges.end())
172 : {
173 : // Ranges should be sorted and should not touch
174 : libmesh_assert_greater(next->first, lastloc+1);
175 : }
176 : #endif
177 :
178 : return (lastloc >= _location);
179 : }
180 :
181 : protected:
182 :
183 : std::size_t _location;
184 : const CouplingMatrix & _mat;
185 : };
186 :
187 :
188 : /**
189 : * This accessor class allows simple setting of CouplingMatrix values.
190 : */
191 : class CouplingAccessor : public ConstCouplingAccessor
192 : {
193 : public:
194 0 : CouplingAccessor(std::size_t loc_in,
195 0 : CouplingMatrix & mat_in) :
196 0 : ConstCouplingAccessor(loc_in, mat_in), _my_mat(mat_in) {}
197 :
198 : template <typename T>
199 0 : CouplingAccessor & operator = (T new_value)
200 : {
201 : // For backwards compatibility we take integer arguments,
202 : // but coupling matrix entries are really all zero or one.
203 0 : const bool as_bool = new_value;
204 0 : libmesh_assert_equal_to(new_value, as_bool);
205 :
206 0 : *this = as_bool;
207 0 : return *this;
208 : }
209 :
210 0 : CouplingAccessor & operator = (bool new_value)
211 : {
212 0 : const std::size_t max_size = std::numeric_limits<std::size_t>::max();
213 :
214 : // Find the range that might contain i,j
215 : // lower_bound isn't *quite* what we want
216 : CouplingMatrix::rc_type::iterator lb =
217 0 : std::upper_bound (_my_mat._ranges.begin(), _my_mat._ranges.end(),
218 0 : std::make_pair(_location, max_size));
219 0 : if (lb!=_my_mat._ranges.begin())
220 0 : --lb;
221 : else
222 0 : lb=_my_mat._ranges.end();
223 :
224 : // If no range might contain i,j then we might need to make a new
225 : // one.
226 0 : if (lb == _my_mat._ranges.end())
227 : {
228 0 : if (new_value == true)
229 0 : _my_mat._ranges.emplace
230 0 : (_my_mat._ranges.begin(), _location, _location);
231 : }
232 : else
233 : {
234 0 : const std::size_t firstloc = lb->first;
235 0 : const std::size_t lastloc = lb->second;
236 0 : libmesh_assert_less_equal(firstloc, lastloc);
237 0 : libmesh_assert_less_equal(firstloc, _location);
238 :
239 : #ifdef DEBUG
240 : {
241 0 : CouplingMatrix::rc_type::const_iterator next = lb;
242 0 : next++;
243 0 : if (next != _my_mat._ranges.end())
244 : {
245 : // Ranges should be sorted and should not touch
246 0 : libmesh_assert_greater(next->first, lastloc+1);
247 : }
248 : }
249 : #endif
250 :
251 : // If we're in this range, we might need to shorten or remove
252 : // or split it
253 0 : if (new_value == false)
254 : {
255 0 : if (_location == firstloc)
256 : {
257 0 : if (_location == lastloc)
258 : {
259 0 : _my_mat._ranges.erase(lb);
260 : }
261 : else
262 : {
263 0 : libmesh_assert_less (lb->first, lastloc);
264 0 : lb->first++;
265 : }
266 : }
267 0 : else if (_location == lastloc)
268 : {
269 0 : libmesh_assert_less (firstloc, lb->second);
270 :
271 0 : lb->second--;
272 : }
273 0 : else if (_location < lastloc)
274 : {
275 0 : libmesh_assert_less_equal(_location+1, lastloc);
276 :
277 0 : lb->first = _location+1;
278 :
279 0 : libmesh_assert_less_equal(firstloc, _location-1);
280 :
281 0 : _my_mat._ranges.emplace(lb, firstloc, _location-1);
282 : }
283 : }
284 :
285 : // If we're not in this range, we might need to extend it or
286 : // join it with its neighbor or add a new one.
287 : else // new_value == true
288 : {
289 0 : CouplingMatrix::rc_type::iterator next = lb;
290 0 : next++;
291 : const std::size_t nextloc =
292 0 : (next == _my_mat._ranges.end()) ?
293 0 : std::numeric_limits<std::size_t>::max() :
294 0 : next->first;
295 :
296 : // Ranges should be sorted and should not touch
297 0 : libmesh_assert_greater(nextloc, lastloc+1);
298 :
299 0 : if (_location > lastloc)
300 : {
301 0 : if (_location == lastloc + 1)
302 : {
303 0 : if (_location == nextloc - 1)
304 : {
305 0 : next->first = firstloc;
306 0 : _my_mat._ranges.erase(lb);
307 : }
308 : else
309 0 : lb->second++;
310 : }
311 : else
312 : {
313 0 : if (_location == nextloc - 1)
314 0 : next->first--;
315 : else
316 0 : _my_mat._ranges.emplace(next, _location, _location);
317 : }
318 : }
319 : }
320 : }
321 :
322 0 : return *this;
323 : }
324 :
325 : private:
326 :
327 : CouplingMatrix & _my_mat;
328 : };
329 :
330 :
331 :
332 : /**
333 : * This proxy class acts like a container of indices from a single
334 : * coupling row
335 : */
336 : class ConstCouplingRow
337 : {
338 : public:
339 360 : ConstCouplingRow(std::size_t row_in,
340 360 : const CouplingMatrix & mat_in) :
341 328 : _row_i(row_in), _mat(mat_in),
342 360 : _end_location(_row_i*_mat.size()+_mat.size()-1)
343 : {
344 0 : libmesh_assert_less(_row_i, _mat.size());
345 :
346 : // Location for i,N, where we'll start looking for our beginning
347 : // range
348 360 : _begin_location = _end_location;
349 :
350 0 : const std::size_t max_size = std::numeric_limits<std::size_t>::max();
351 :
352 : // Find the range that might contain i,N
353 : // lower_bound isn't *quite* what we want
354 : _begin_it = std::upper_bound
355 328 : (_mat._ranges.begin(), _mat._ranges.end(),
356 360 : std::make_pair(_begin_location, max_size));
357 360 : if (_begin_it !=_mat._ranges.begin())
358 0 : --_begin_it;
359 : else
360 0 : _begin_it=_mat._ranges.end();
361 :
362 : // If that range doesn't exist then we're an empty row
363 392 : if (_begin_it == _mat._ranges.end())
364 0 : _begin_location = max_size;
365 : else
366 : {
367 360 : const std::size_t lastloc = _begin_it->second;
368 : #ifdef DEBUG
369 0 : const std::size_t firstloc = _begin_it->first;
370 0 : libmesh_assert_less_equal(firstloc, lastloc);
371 : #endif
372 :
373 : // If that range ends before i,0 then we're an empty row
374 32 : std::size_t zero_location = _row_i*_mat.size();
375 360 : if (zero_location > lastloc)
376 : {
377 90 : _begin_location = max_size;
378 90 : _begin_it = _mat._ranges.end();
379 : }
380 : else
381 : // We have *some* entry(s) in this row, we just need to find
382 : // the earliest
383 : {
384 294 : while (_begin_it != _mat._ranges.begin())
385 : {
386 0 : CouplingMatrix::rc_type::const_iterator prev =
387 : _begin_it;
388 0 : --prev;
389 :
390 180 : if (prev->second < zero_location)
391 0 : break;
392 :
393 0 : _begin_it = prev;
394 : }
395 270 : if (_begin_it->first < zero_location)
396 90 : _begin_location = zero_location;
397 : else
398 180 : _begin_location = _begin_it->first;
399 : }
400 : }
401 360 : }
402 :
403 : /*
404 : * A forward iterator type for looping over indices in this row
405 : */
406 : typedef ConstCouplingRowConstIterator const_iterator;
407 :
408 : /*
409 : * An iterator to the first index in this row, or to end() for an
410 : * empty row
411 : */
412 : const_iterator begin() const;
413 :
414 : /*
415 : * An iterator representing past-the-end of this row
416 : */
417 : const_iterator end() const;
418 :
419 0 : bool operator== (const ConstCouplingRow & other) const
420 : {
421 : // Thinking that rows from different matrix objects are equal is
422 : // not even wrong
423 0 : libmesh_assert(&_mat == &other._mat);
424 :
425 0 : return ((_begin_location == other._begin_location) &&
426 0 : (_begin_it == other._begin_it));
427 : }
428 :
429 : bool operator!= (const ConstCouplingRow & other) const
430 : {
431 : return !(*this == other);
432 : }
433 : protected:
434 :
435 : friend class ConstCouplingRowConstIterator;
436 :
437 : std::size_t _row_i;
438 : const CouplingMatrix & _mat;
439 :
440 : // The location (i*size+j) corresponding to the first entry in this
441 : // row, or numeric_limits<size_t>::max() for an empty row.
442 : std::size_t _begin_location;
443 :
444 : // The location (i*size+j) corresponding to the end of this row
445 : // (regardless of whether that end falls within a range). Cached
446 : // because every ++iterator call needs to use it.
447 : const std::size_t _end_location;
448 :
449 : // Iterator to the range containing the first row element, or
450 : // _row._mat._ranges.end() if no CouplingMatrix values are true for
451 : // this row
452 : CouplingMatrix::rc_type::const_iterator _begin_it;
453 : };
454 :
455 :
456 :
457 : class ConstCouplingRowConstIterator
458 : {
459 : public:
460 0 : ConstCouplingRowConstIterator (const ConstCouplingRow & row_in,
461 : std::size_t loc_in,
462 360 : CouplingMatrix::rc_type::const_iterator it_in) :
463 328 : _location(loc_in),
464 328 : _row(row_in),
465 360 : _it(it_in)
466 : {
467 : #ifndef NDEBUG
468 0 : if (_it != _row._mat._ranges.end())
469 : {
470 0 : libmesh_assert_less_equal(_it->first, _location);
471 0 : libmesh_assert_less_equal(_location, _it->second);
472 : }
473 : else
474 : {
475 0 : libmesh_assert_equal_to
476 : (_location, std::numeric_limits<size_t>::max());
477 : }
478 : #endif
479 0 : }
480 :
481 0 : std::size_t operator* ()
482 : {
483 0 : libmesh_assert_not_equal_to
484 : (_location, std::numeric_limits<std::size_t>::max());
485 450 : return _location % _row._mat.size();
486 : }
487 :
488 450 : ConstCouplingRowConstIterator & operator++ ()
489 : {
490 0 : libmesh_assert_not_equal_to
491 : (_location, std::numeric_limits<std::size_t>::max());
492 0 : libmesh_assert
493 : (_it != _row._mat._ranges.end());
494 :
495 450 : if (_location >= _it->second)
496 : {
497 0 : ++_it;
498 :
499 : // Are we past the end of the matrix?
500 180 : if (_it == _row._mat._ranges.end())
501 90 : _location = std::numeric_limits<std::size_t>::max();
502 : else
503 : {
504 90 : _location = _it->first;
505 : // Is the new range past the end of the row?
506 90 : if (_location > _row._end_location)
507 90 : _location = std::numeric_limits<std::size_t>::max();
508 : }
509 : }
510 : // Would incrementing put us past the end of the row?
511 270 : else if (_location >= _row._end_location)
512 90 : _location = std::numeric_limits<std::size_t>::max();
513 : else
514 180 : ++_location;
515 :
516 450 : return *this;
517 : }
518 :
519 0 : bool operator== (const ConstCouplingRowConstIterator & other) const
520 : {
521 : // Thinking that iterators from different row objects are equal
522 : // is not even wrong
523 0 : libmesh_assert(_row == other._row);
524 :
525 738 : return (_location == other._location);
526 : // Testing for equality of _it is either redundant (for a
527 : // dereferenceable iterator, where _it is defined to be the range
528 : // which includes _location) or wrong (for an end iterator, where
529 : // we don't always bother to set _it)
530 : // && (_it == other._it);
531 : }
532 :
533 0 : bool operator!= (const ConstCouplingRowConstIterator & other) const
534 : {
535 72 : return !(*this == other);
536 : }
537 :
538 : private:
539 : // The location (i*size+j) corresponding to this iterator, or
540 : // numeric_limits<size_t>::max() to signify end()
541 : std::size_t _location;
542 : const ConstCouplingRow & _row;
543 : // The range containing this iterator location, or
544 : // _row._mat._ranges.end() if no range contains this location
545 : CouplingMatrix::rc_type::const_iterator _it;
546 : };
547 :
548 :
549 :
550 : //--------------------------------------------------
551 : // ConstCouplingRow inline methods
552 : inline
553 0 : ConstCouplingRow::const_iterator ConstCouplingRow::begin() const
554 : {
555 360 : return const_iterator (*this, _begin_location, _begin_it);
556 : }
557 :
558 : inline
559 0 : ConstCouplingRow::const_iterator ConstCouplingRow::end() const
560 : {
561 : return const_iterator
562 : (*this, std::numeric_limits<std::size_t>::max(),
563 328 : _mat._ranges.end());
564 : }
565 :
566 :
567 :
568 : //--------------------------------------------------
569 : // CouplingMatrix inline methods
570 : inline
571 : CouplingMatrix::CouplingMatrix (const std::size_t n) :
572 : _ranges(), _size(n)
573 : {
574 : this->resize(n);
575 : }
576 :
577 :
578 :
579 : inline
580 : bool CouplingMatrix::operator() (const std::size_t i,
581 : const std::size_t j) const
582 : {
583 : libmesh_assert_less (i, _size);
584 : libmesh_assert_less (j, _size);
585 :
586 : const std::size_t location = std::size_t(i)*_size + j;
587 :
588 : return bool(ConstCouplingAccessor(location, *this));
589 : }
590 :
591 :
592 :
593 :
594 : inline
595 0 : CouplingAccessor CouplingMatrix::operator() (const std::size_t i,
596 : const std::size_t j)
597 : {
598 0 : const std::size_t location = std::size_t(i)*_size + j;
599 :
600 0 : return CouplingAccessor(location, *this);
601 : }
602 :
603 :
604 :
605 : inline
606 0 : std::size_t CouplingMatrix::size() const
607 : {
608 738 : return _size;
609 : }
610 :
611 :
612 :
613 : inline
614 0 : void CouplingMatrix::resize(const std::size_t n)
615 : {
616 0 : _size = n;
617 :
618 0 : _ranges.clear();
619 0 : }
620 :
621 :
622 :
623 : inline
624 : void CouplingMatrix::clear()
625 : {
626 : this->resize(0);
627 : }
628 :
629 :
630 :
631 : inline
632 0 : bool CouplingMatrix::empty() const
633 : {
634 0 : return (_size == 0);
635 : }
636 :
637 :
638 : } // namespace libMesh
639 :
640 :
641 : #endif // LIBMESH_COUPLING_MATRIX_H
|
