Dune Core Modules (2.8.0)

multitypeblockvector.hh
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
12#include <dune/common/hybridutilities.hh>
14
15#include "istlexception.hh"
16
17// forward declaration
18namespace Dune {
19 template < typename... Args >
20 class MultiTypeBlockVector;
21}
22
23#include "gsetc.hh"
24
25namespace 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
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
A Vector class to support different block types.
Definition: multitypeblockvector.hh:56
Provides the functions dot(a,b) := and dotT(a,b) := .
Traits for type conversions and type information.
Type traits to determine the type of reals (when working with complex numbers)
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:266
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:182
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
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[](const std::integral_constant< size_type, index > indexVariable)
Random-access operator.
Definition: multitypeblockvector.hh:142
FieldTraits< field_type >::real_type two_norm() const
Compute the Euclidean norm.
Definition: multitypeblockvector.hh:249
size_type dim() const
Number of scalar elements.
Definition: multitypeblockvector.hh:113
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
FieldTraits< field_type >::real_type infinity_norm_real() const
Compute the simplified maximum norm (uses 1-norm for complex values)
Definition: multitypeblockvector.hh:284
auto one_norm() const
Compute the 1-norm.
Definition: multitypeblockvector.hh:222
void operator-=(const type &newv)
Definition: multitypeblockvector.hh:181
void operator+=(const type &newv)
Definition: multitypeblockvector.hh:171
FieldTraits< field_type >::real_type infinity_norm() const
Compute the maximum norm.
Definition: multitypeblockvector.hh:253
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
FieldTraits< field_type >::real_type two_norm2() const
Compute the squared Euclidean norm.
Definition: multitypeblockvector.hh:240
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:11
Whether this type has a value of NaN.
Definition: typetraits.hh:212
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)