DUNE-FEM (unstable)

multitypeblockvector.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
6#define DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
7
8#include <cmath>
9#include <iostream>
10#include <tuple>
11
14#include <dune/common/hybridutilities.hh>
16#include <dune/common/std/type_traits.hh>
17
18#include "istlexception.hh"
19
20// forward declaration
21namespace Dune {
22 template < typename... Args >
23 class MultiTypeBlockVector;
24}
25
26#include "gsetc.hh"
27
28namespace Dune {
29
41 template <typename... Args>
42 struct FieldTraits< MultiTypeBlockVector<Args...> >
43 {
44 using field_type = typename MultiTypeBlockVector<Args...>::field_type;
45 using real_type = typename MultiTypeBlockVector<Args...>::real_type;
46 };
56 template < typename... Args >
58 : public std::tuple<Args...>
59 {
61 typedef std::tuple<Args...> TupleType;
62 public:
63
65 using size_type = std::size_t;
66
70 using TupleType::TupleType;
71
75 typedef MultiTypeBlockVector<Args...> type;
76
83
90
91 // make sure that we have an std::common_type: using an assertion produces a more readable error message
92 // than a compiler template instantiation error
93 static_assert ( sizeof...(Args) == 0 or
94 not (std::is_same_v<field_type, Std::nonesuch> or std::is_same_v<real_type, Std::nonesuch>),
95 "No std::common_type implemented for the given field_type/real_type of the Args. Please provide an implementation and check that a FieldTraits class is present for your type.");
96
97
103 static constexpr size_type size()
104 {
105 return sizeof...(Args);
106 }
107
110 static constexpr size_type N()
111 {
112 return sizeof...(Args);
113 }
114
117 {
118 size_type result = 0;
119 Hybrid::forEach(std::make_index_sequence<N()>{},
120 [&](auto i){result += std::get<i>(*this).dim();});
121
122 return result;
123 }
124
143 template< size_type index >
144 typename std::tuple_element<index,TupleType>::type&
145 operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
146 {
147 return std::get<index>(*this);
148 }
149
155 template< size_type index >
156 const typename std::tuple_element<index,TupleType>::type&
157 operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) const
158 {
159 return std::get<index>(*this);
160 }
161
164 template<typename T>
165 void operator= (const T& newval) {
166 Dune::Hybrid::forEach(*this, [&](auto&& entry) {
167 entry = newval;
168 });
169 }
170
174 void operator+= (const type& newv) {
175 using namespace Dune::Hybrid;
176 forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
177 (*this)[i] += newv[i];
178 });
179 }
180
184 void operator-= (const type& newv) {
185 using namespace Dune::Hybrid;
186 forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
187 (*this)[i] -= newv[i];
188 });
189 }
190
192 template<class T,
193 std::enable_if_t< IsNumber<T>::value, int> = 0>
194 void operator*= (const T& w) {
195 Hybrid::forEach(*this, [&](auto&& entry) {
196 entry *= w;
197 });
198 }
199
201 template<class T,
202 std::enable_if_t< IsNumber<T>::value, int> = 0>
203 void operator/= (const T& w) {
204 Hybrid::forEach(*this, [&](auto&& entry) {
205 entry /= w;
206 });
207 }
208
209 field_type operator* (const type& newv) const {
210 using namespace Dune::Hybrid;
211 return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
212 return a + (*this)[i]*newv[i];
213 });
214 }
215
216 field_type dot (const type& newv) const {
217 using namespace Dune::Hybrid;
218 return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
219 return a + (*this)[i].dot(newv[i]);
220 });
221 }
222
225 auto one_norm() const {
226 using namespace Dune::Hybrid;
227 return accumulate(*this, real_type(0), [&](auto&& a, auto&& entry) {
228 return a + entry.one_norm();
229 });
230 }
231
234 auto one_norm_real() const {
235 using namespace Dune::Hybrid;
236 return accumulate(*this, real_type(0), [&](auto&& a, auto&& entry) {
237 return a + entry.one_norm_real();
238 });
239 }
240
244 using namespace Dune::Hybrid;
245 return accumulate(*this, real_type(0), [&](auto&& a, auto&& entry) {
246 return a + entry.two_norm2();
247 });
248 }
249
252 real_type two_norm() const {return sqrt(this->two_norm2());}
253
257 {
258 using namespace Dune::Hybrid;
259 using std::max;
260
261 real_type result = 0.0;
262 // Compute max norm tracking appearing nan values
263 // if the field type supports nan.
264 if constexpr (HasNaN<field_type>()) {
265 // This variable will preserve any nan value
266 real_type nanTracker = 1.0;
267 using namespace Dune::Hybrid; // needed for icc, see issue #31
268 forEach(*this, [&](auto&& entry) {
269 real_type entryNorm = entry.infinity_norm();
270 result = max(entryNorm, result);
271 nanTracker += entryNorm;
272 });
273 // Incorporate possible nan value into result
274 result *= (nanTracker / nanTracker);
275 } else {
276 using namespace Dune::Hybrid; // needed for icc, see issue #31
277 forEach(*this, [&](auto&& entry) {
278 result = max(entry.infinity_norm(), result);
279 });
280 }
281 return result;
282 }
283
287 {
288 using namespace Dune::Hybrid;
289 using std::max;
290
291 real_type result = 0.0;
292 // Compute max norm tracking appearing nan values
293 // if the field type supports nan.
294 if constexpr (HasNaN<field_type>()) {
295 // This variable will preserve any nan value
296 real_type nanTracker = 1.0;
297 using namespace Dune::Hybrid; // needed for icc, see issue #31
298 forEach(*this, [&](auto&& entry) {
299 real_type entryNorm = entry.infinity_norm_real();
300 result = max(entryNorm, result);
301 nanTracker += entryNorm;
302 });
303 // Incorporate possible nan value into result
304 result *= (nanTracker / nanTracker);
305 } else {
306 using namespace Dune::Hybrid; // needed for icc, see issue #31
307 forEach(*this, [&](auto&& entry) {
308 result = max(entry.infinity_norm_real(), result);
309 });
310 }
311 return result;
312 }
313
318 template<typename Ta>
319 void axpy (const Ta& a, const type& y) {
320 using namespace Dune::Hybrid;
321 forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
322 (*this)[i].axpy(a, y[i]);
323 });
324 }
325
326 };
327
328
329
332 template <typename... Args>
333 std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<Args...>& v) {
334 using namespace Dune::Hybrid;
335 forEach(integralRange(Dune::Hybrid::size(v)), [&](auto&& i) {
336 s << "\t(" << i << "):\n" << v[i] << "\n";
337 });
338 return s;
339 }
340
341} // end namespace Dune
342
343namespace std
344{
349 template <size_t i, typename... Args>
350 struct tuple_element<i,Dune::MultiTypeBlockVector<Args...> >
351 {
352 using type = typename std::tuple_element<i, std::tuple<Args...> >::type;
353 };
354
359 template <typename... Args>
360 struct tuple_size<Dune::MultiTypeBlockVector<Args...> >
361 : std::integral_constant<std::size_t, sizeof...(Args)>
362 {};
363}
364
365#endif
A Vector class to support different block types.
Definition: multitypeblockvector.hh:59
Provides the functions dot(a,b) := and dotT(a,b) := .
Type traits to determine the type of reals (when working with complex numbers)
typename detected_or< nonesuch, Op, Args... >::type detected_t
Returns Op<Args...> if that is valid; otherwise returns nonesuch.
Definition: type_traits.hh:174
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:172
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
void operator=(const T &newval)
Assignment operator.
Definition: multitypeblockvector.hh:165
std::size_t size_type
Type used for vector sizes.
Definition: multitypeblockvector.hh:65
static constexpr size_type N()
Number of elements.
Definition: multitypeblockvector.hh:110
Std::detected_t< std::common_type_t, typename FieldTraits< std::decay_t< Args > >::field_type... > field_type
The type used for scalars.
Definition: multitypeblockvector.hh:82
static constexpr size_type size()
Return the number of non-zero vector entries.
Definition: multitypeblockvector.hh:103
std::tuple_element< index, TupleType >::type & operator[](const std::integral_constant< size_type, index > indexVariable)
Random-access operator.
Definition: multitypeblockvector.hh:145
size_type dim() const
Number of scalar elements.
Definition: multitypeblockvector.hh:116
Std::detected_t< std::common_type_t, typename FieldTraits< std::decay_t< Args > >::real_type... > real_type
The type used for real values.
Definition: multitypeblockvector.hh:89
void operator*=(const T &w)
Multiplication with a scalar.
Definition: multitypeblockvector.hh:194
void axpy(const Ta &a, const type &y)
Axpy operation on this vector (*this += a * y)
Definition: multitypeblockvector.hh:319
real_type two_norm() const
Compute the Euclidean norm.
Definition: multitypeblockvector.hh:252
void operator/=(const T &w)
Division by a scalar.
Definition: multitypeblockvector.hh:203
MultiTypeBlockVector< Args... > type
Definition: multitypeblockvector.hh:75
auto one_norm() const
Compute the 1-norm.
Definition: multitypeblockvector.hh:225
real_type two_norm2() const
Compute the squared Euclidean norm.
Definition: multitypeblockvector.hh:243
real_type infinity_norm_real() const
Compute the simplified maximum norm (uses 1-norm for complex values)
Definition: multitypeblockvector.hh:286
void operator-=(const type &newv)
Definition: multitypeblockvector.hh:184
real_type infinity_norm() const
Compute the maximum norm.
Definition: multitypeblockvector.hh:256
void operator+=(const type &newv)
Definition: multitypeblockvector.hh:174
auto one_norm_real() const
Compute the simplified 1-norm (uses 1-norm also for complex values)
Definition: multitypeblockvector.hh:234
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
Whether this type has a value of NaN.
Definition: typetraits.hh:212
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)