libMesh
coupling_matrix.C
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 #include "libmesh/coupling_matrix.h"
20 
21 
22 namespace libMesh {
23 
25 {
26  const std::size_t max_size = std::numeric_limits<std::size_t>::max();
27 
28  rc_type::iterator start_range = this->_ranges.begin();
29 
30  rc_type::const_iterator other_range = other._ranges.begin();
31  const rc_type::const_iterator other_end = other._ranges.end();
32 
33  for (; other_range != other_end; ++other_range)
34  {
35  std::size_t other_range_start = other_range->first;
36  std::size_t other_range_end = other_range->second;
37 
38  // Find our range that might contain the start of the other
39  // range.
40  // lower_bound isn't *quite* what we want.
41  // Because the other._ranges is sorted, we can contract this
42  // search as we proceed, beginning with lb rather than at
43  // begin() every time.
44  rc_type::iterator lb =
45  std::upper_bound (start_range, this->_ranges.end(),
46  std::make_pair(other_range_start, max_size));
47  if (lb!=start_range)
48  --lb;
49  else
50  lb=this->_ranges.end();
51 
52  start_range = lb;
53 
54  // If no range might contain the start of the new range then
55  // we can just break out of here and start appending any
56  // remaining ranges.
57  if (lb == this->_ranges.end())
58  break;
59 
60  // We did find a range which might contain the start of the new
61  // range.
62  const std::size_t lastloc = lb->second;
63  libmesh_assert_less_equal(lb->first, lastloc);
64  libmesh_assert_less_equal(lb->first, other_range_start);
65 
66 #ifdef DEBUG
67  {
68  CouplingMatrix::rc_type::const_iterator next = lb;
69  next++;
70  if (next != this->_ranges.end())
71  {
72  // Ranges should be sorted and should not touch
73  libmesh_assert_greater(next->first, lastloc+1);
74  }
75  }
76 #endif
77 
78  CouplingMatrix::rc_type::iterator next = lb;
79  next++;
80 
81  // We might need to extend this range or add a new range.
82  // Merge contiguous ranges
83  if (other_range_start <= lastloc)
84  lb->second = other_range_end;
85  // Or insert a new range. This invalidates existing iterators,
86  // but that's okay; the only iterator we need is for the newly
87  // inserted range.
88  else
89  start_range = lb = this->_ranges.insert
90  (next, std::make_pair(other_range_start, other_range_end));
91 
92  // At this point we have a range lb that may potentially overlap
93  // subsequent existing ranges, in which case we need to merge
94  // some.
95 
96  // First expand our range as necessary while finding ranges
97  // which will be redundant later
98  for (const std::size_t nextloc =
99  (next == this->_ranges.end()) ?
101  nextloc <= lb->second; ++next)
102  {
103  // Ranges should be sorted and should not have been touching
104  // initially
105  libmesh_assert_greater(nextloc, lastloc+1);
106 
107  lb->second = std::max(lb->second, next->second);
108  }
109 
110  CouplingMatrix::rc_type::iterator oldnext = lb;
111  oldnext++;
112 
113  // Finally remove the redundant ranges
114  this->_ranges.erase(oldnext, next);
115  }
116 
117  // If we broke out early but we still have more other_ranges then
118  // we can safely just append them to our ranges.
119  for (; other_range != other_end; ++other_range)
120  this->_ranges.push_back(*other_range);
121 
122  // Behave like a standard modification operator
123  return *this;
124 }
125 
126 }
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
CouplingMatrix & operator&=(const CouplingMatrix &other)
static const unsigned int next[3]
A lookup table for the increment modulo 3 operation, for iterating through the three nodes per elemen...
This class defines a coupling matrix.