dune-istl  2.8.0
multitypeblockvector.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
4 #define DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
5 
6 #include <cmath>
7 #include <iostream>
8 #include <tuple>
9 
10 #include <dune/common/dotproduct.hh>
11 #include <dune/common/ftraits.hh>
12 #include <dune/common/hybridutilities.hh>
13 #include <dune/common/typetraits.hh>
14 
15 #include "istlexception.hh"
16 
17 // forward declaration
18 namespace Dune {
19  template < typename... Args >
20  class MultiTypeBlockVector;
21 }
22 
23 #include "gsetc.hh"
24 
25 namespace Dune {
26 
38  template <typename Arg0, typename... Args>
39  struct FieldTraits< MultiTypeBlockVector<Arg0, Args...> >
40  {
41  typedef typename FieldTraits<Arg0>::field_type field_type;
42  typedef typename FieldTraits<Arg0>::real_type real_type;
43  };
53  template < typename... Args >
55  : public std::tuple<Args...>
56  {
58  typedef std::tuple<Args...> TupleType;
59  public:
60 
62  using size_type = std::size_t;
63 
67  using TupleType::TupleType;
68 
72  typedef MultiTypeBlockVector<Args...> type;
73 
81  typedef double field_type;
82 
88  static constexpr size_type size()
89  {
90  return sizeof...(Args);
91  }
92 
95  static constexpr size_type N()
96  {
97  return sizeof...(Args);
98  }
99 
106  [[deprecated("Use method 'N' instead")]]
107  int count() const
108  {
109  return sizeof...(Args);
110  }
111 
113  size_type dim() const
114  {
115  size_type result = 0;
116  Hybrid::forEach(std::make_index_sequence<N()>{},
117  [&](auto i){result += std::get<i>(*this).dim();});
118 
119  return result;
120  }
121 
140  template< size_type index >
141  typename std::tuple_element<index,TupleType>::type&
142  operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
143  {
144  return std::get<index>(*this);
145  }
146 
152  template< size_type index >
153  const typename std::tuple_element<index,TupleType>::type&
154  operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) const
155  {
156  return std::get<index>(*this);
157  }
158 
161  template<typename T>
162  void operator= (const T& newval) {
163  Dune::Hybrid::forEach(*this, [&](auto&& entry) {
164  entry = newval;
165  });
166  }
167 
171  void operator+= (const type& newv) {
172  using namespace Dune::Hybrid;
173  forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
174  (*this)[i] += newv[i];
175  });
176  }
177 
181  void operator-= (const type& newv) {
182  using namespace Dune::Hybrid;
183  forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
184  (*this)[i] -= newv[i];
185  });
186  }
187 
189  template<class T,
190  std::enable_if_t< IsNumber<T>::value, int> = 0>
191  void operator*= (const T& w) {
192  Hybrid::forEach(*this, [&](auto&& entry) {
193  entry *= w;
194  });
195  }
196 
198  template<class T,
199  std::enable_if_t< IsNumber<T>::value, int> = 0>
200  void operator/= (const T& w) {
201  Hybrid::forEach(*this, [&](auto&& entry) {
202  entry /= w;
203  });
204  }
205 
206  field_type operator* (const type& newv) const {
207  using namespace Dune::Hybrid;
208  return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
209  return a + (*this)[i]*newv[i];
210  });
211  }
212 
213  field_type dot (const type& newv) const {
214  using namespace Dune::Hybrid;
215  return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
216  return a + (*this)[i].dot(newv[i]);
217  });
218  }
219 
222  auto one_norm() const {
223  using namespace Dune::Hybrid;
224  return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
225  return a + entry.one_norm();
226  });
227  }
228 
231  auto one_norm_real() const {
232  using namespace Dune::Hybrid;
233  return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
234  return a + entry.one_norm_real();
235  });
236  }
237 
240  typename FieldTraits<field_type>::real_type two_norm2() const {
241  using namespace Dune::Hybrid;
242  return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
243  return a + entry.two_norm2();
244  });
245  }
246 
249  typename FieldTraits<field_type>::real_type two_norm() const {return sqrt(this->two_norm2());}
250 
253  typename FieldTraits<field_type>::real_type infinity_norm() const
254  {
255  using namespace Dune::Hybrid;
256  using std::max;
257  using real_type = typename FieldTraits<field_type>::real_type;
258 
259  real_type result = 0.0;
260  // Compute max norm tracking appearing nan values
261  // if the field type supports nan.
262  if constexpr (HasNaN<field_type>()) {
263  // This variable will preserve any nan value
264  real_type nanTracker = 1.0;
265  using namespace Dune::Hybrid; // needed for icc, see issue #31
266  forEach(*this, [&](auto&& entry) {
267  real_type entryNorm = entry.infinity_norm();
268  result = max(entryNorm, result);
269  nanTracker += entryNorm;
270  });
271  // Incorporate possible nan value into result
272  result *= (nanTracker / nanTracker);
273  } else {
274  using namespace Dune::Hybrid; // needed for icc, see issue #31
275  forEach(*this, [&](auto&& entry) {
276  result = max(entry.infinity_norm(), result);
277  });
278  }
279  return result;
280  }
281 
284  typename FieldTraits<field_type>::real_type infinity_norm_real() const
285  {
286  using namespace Dune::Hybrid;
287  using std::max;
288  using real_type = typename FieldTraits<field_type>::real_type;
289 
290  real_type result = 0.0;
291  // Compute max norm tracking appearing nan values
292  // if the field type supports nan.
293  if constexpr (HasNaN<field_type>()) {
294  // This variable will preserve any nan value
295  real_type nanTracker = 1.0;
296  using namespace Dune::Hybrid; // needed for icc, see issue #31
297  forEach(*this, [&](auto&& entry) {
298  real_type entryNorm = entry.infinity_norm_real();
299  result = max(entryNorm, result);
300  nanTracker += entryNorm;
301  });
302  // Incorporate possible nan value into result
303  result *= (nanTracker / nanTracker);
304  } else {
305  using namespace Dune::Hybrid; // needed for icc, see issue #31
306  forEach(*this, [&](auto&& entry) {
307  result = max(entry.infinity_norm_real(), result);
308  });
309  }
310  return result;
311  }
312 
317  template<typename Ta>
318  void axpy (const Ta& a, const type& y) {
319  using namespace Dune::Hybrid;
320  forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
321  (*this)[i].axpy(a, y[i]);
322  });
323  }
324 
325  };
326 
327 
328 
331  template <typename... Args>
332  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<Args...>& v) {
333  using namespace Dune::Hybrid;
334  forEach(integralRange(Dune::Hybrid::size(v)), [&](auto&& i) {
335  s << "\t(" << i << "):\n" << v[i] << "\n";
336  });
337  return s;
338  }
339 
340 
341 
342 } // end namespace
343 
344 #endif
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
FieldTraits< Arg0 >::field_type field_type
Definition: multitypeblockvector.hh:41
FieldTraits< Arg0 >::real_type real_type
Definition: multitypeblockvector.hh:42
void operator=(const T &newval)
Assignment operator.
Definition: multitypeblockvector.hh:162
std::size_t size_type
Type used for vector sizes.
Definition: multitypeblockvector.hh:62
int count() const
Definition: multitypeblockvector.hh:107
static constexpr size_type N()
Number of elements.
Definition: multitypeblockvector.hh:95
static constexpr size_type size()
Return the number of non-zero vector entries.
Definition: multitypeblockvector.hh:88
std::tuple_element< index, TupleType >::type & operator[]([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
Random-access operator.
Definition: multitypeblockvector.hh:142
size_type dim() const
Number of scalar elements.
Definition: multitypeblockvector.hh:113
field_type dot(const type &newv) const
Definition: multitypeblockvector.hh:213
void operator*=(const T &w)
Multiplication with a scalar.
Definition: multitypeblockvector.hh:191
void axpy(const Ta &a, const type &y)
Axpy operation on this vector (*this += a * y)
Definition: multitypeblockvector.hh:318
void operator/=(const T &w)
Division by a scalar.
Definition: multitypeblockvector.hh:200
MultiTypeBlockVector< Args... > type
Definition: multitypeblockvector.hh:72
auto one_norm() const
Compute the 1-norm.
Definition: multitypeblockvector.hh:222
void operator-=(const type &newv)
Definition: multitypeblockvector.hh:181
FieldTraits< field_type >::real_type infinity_norm_real() const
Compute the simplified maximum norm (uses 1-norm for complex values)
Definition: multitypeblockvector.hh:284
field_type operator*(const type &newv) const
Definition: multitypeblockvector.hh:206
FieldTraits< field_type >::real_type infinity_norm() const
Compute the maximum norm.
Definition: multitypeblockvector.hh:253
void operator+=(const type &newv)
Definition: multitypeblockvector.hh:171
FieldTraits< field_type >::real_type two_norm2() const
Compute the squared Euclidean norm.
Definition: multitypeblockvector.hh:240
FieldTraits< field_type >::real_type two_norm() const
Compute the Euclidean norm.
Definition: multitypeblockvector.hh:249
auto one_norm_real() const
Compute the simplified 1-norm (uses 1-norm also for complex values)
Definition: multitypeblockvector.hh:231
double field_type
The type used for scalars.
Definition: multitypeblockvector.hh:81
Definition: allocator.hh:9
std::ostream & operator<<(std::ostream &s, const BlockVector< K, A > &v)
Send BlockVector to an output stream.
Definition: bvector.hh:588
A Vector class to support different block types.
Definition: multitypeblockvector.hh:56