Dune Core Modules (2.4.1)

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#if HAVE_DUNE_BOOST
7#ifdef HAVE_BOOST_FUSION
8
9#include <cmath>
10#include <iostream>
11
14
15#include "istlexception.hh"
16
17#include <boost/fusion/sequence.hpp>
18#include <boost/fusion/container.hpp>
19#include <boost/fusion/iterator.hpp>
20#include <boost/typeof/typeof.hpp>
21#include <boost/fusion/algorithm.hpp>
22
23namespace mpl=boost::mpl;
24namespace fusion=boost::fusion;
25
26// forward decl
27namespace Dune {
28 template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
29 typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
30 typename T8=fusion::void_, typename T9=fusion::void_>
31 class MultiTypeBlockVector;
32}
33
34#include "gsetc.hh"
35
36namespace Dune {
37
61 template<int current_element, int remaining_elements, typename TVec>
62 class MultiTypeBlockVector_Print {
63 public:
67 static void print(const TVec& v) {
68 std::cout << "\t(" << current_element << "):\n" << fusion::at_c<current_element>(v) << "\n";
69 MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v); //next element
70 }
71 };
72 template<int current_element, typename TVec> //recursion end (remaining_elements=0)
73 class MultiTypeBlockVector_Print<current_element,0,TVec> {
74 public:
75 static void print(const TVec& v) {std::cout << "\n";}
76 };
77
78
79
91 template<int count, typename T1, typename T2>
92 class MultiTypeBlockVector_Ident {
93 public:
94
99 static void equalize(T1& a, const T2& b) {
100 fusion::at_c<count-1>(a) = b; //equalize current elements
101 MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b); //next elements
102 }
103 };
104 template<typename T1, typename T2> //recursion end (count=0)
105 class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1& a, const T2& b) {} };
106
107
108
109
110
111
117 template<int count, typename T>
118 class MultiTypeBlockVector_Add {
119 public:
120
124 static void add (T& a, const T& b) { //add vector elements
125 fusion::at_c<(count-1)>(a) += fusion::at_c<(count-1)>(b);
126 MultiTypeBlockVector_Add<count-1,T>::add(a,b);
127 }
128
132 static void sub (T& a, const T& b) { //sub vector elements
133 fusion::at_c<(count-1)>(a) -= fusion::at_c<(count-1)>(b);
134 MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
135 }
136 };
137 template<typename T> //recursion end; specialization for count=0
138 class MultiTypeBlockVector_Add<0,T> {public: static void add (T& a, const T& b) {} static void sub (T& a, const T& b) {} };
139
140
141
147 template<int count, typename TVec, typename Ta>
148 class MultiTypeBlockVector_AXPY {
149 public:
150
154 static void axpy(TVec& x, const Ta& a, const TVec& y) {
155 fusion::at_c<(count-1)>(x).axpy(a,fusion::at_c<(count-1)>(y));
156 MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
157 }
158 };
159 template<typename TVec, typename Ta> //specialization for count=0
160 class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec& x, const Ta& a, const TVec& y) {} };
161
162
168 template<int count, typename TVec, typename Ta>
169 class MultiTypeBlockVector_Mulscal {
170 public:
171
175 static void mul(TVec& x, const Ta& a) {
176 fusion::at_c<(count-1)>(x) *= a;
177 MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
178 }
179 };
180 template<typename TVec, typename Ta> //specialization for count=0
181 class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec& x, const Ta& a) {} };
182
183
184
193 template<int count, typename TVec>
194 class MultiTypeBlockVector_Mul {
195 public:
196 static typename TVec::field_type mul(const TVec& x, const TVec& y) { return (fusion::at_c<count-1>(x) * fusion::at_c<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y); }
197 static typename TVec::field_type dot(const TVec& x, const TVec& y) { return (Dune::dot(fusion::at_c<count-1>(x),fusion::at_c<count-1>(y))) + MultiTypeBlockVector_Mul<count-1,TVec>::dot(x,y); }
198 };
199 template<typename TVec>
200 class MultiTypeBlockVector_Mul<0,TVec> {
201 public:
202 static typename TVec::field_type mul(const TVec& x, const TVec& y) {return 0;}
203 static typename TVec::field_type dot(const TVec& x, const TVec& y) {return 0;}
204 };
205
206
207
208
209
216 template<int count, typename T>
217 class MultiTypeBlockVector_Norm {
218 public:
219 typedef typename T::field_type field_type;
220 typedef typename FieldTraits<field_type>::real_type real_type;
221
225 static real_type result (const T& a) { //result = sum of all elements' 2-norms
226 return fusion::at_c<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
227 }
228 };
229
230 template<typename T> //recursion end: no more vector elements to add...
231 class MultiTypeBlockVector_Norm<0,T> {
232 public:
233 typedef typename T::field_type field_type;
234 typedef typename FieldTraits<field_type>::real_type real_type;
235 static real_type result (const T& a) {return 0.0;}
236 };
237
246 template<typename T1, typename T2, typename T3, typename T4,
247 typename T5, typename T6, typename T7, typename T8, typename T9>
248 class MultiTypeBlockVector : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
249
250 public:
251
255 typedef MultiTypeBlockVector<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
256
257 typedef typename T1::field_type field_type;
258
262 int count() {return mpl::size<type>::value;}
263
282 template< std::size_t index >
283 auto
284 operator[] ( const std::integral_constant< std::size_t, index > indexVariable )
285 -> decltype(fusion::at_c<index>(*this))
286 {
287 DUNE_UNUSED_PARAMETER(indexVariable);
288 return fusion::at_c<index>(*this);
289 }
290
296 template< std::size_t index >
297 auto
298 operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const
299 -> decltype(fusion::at_c<index>(*this))
300 {
301 DUNE_UNUSED_PARAMETER(indexVariable);
302 return fusion::at_c<index>(*this);
303 }
304
308 template<typename T>
309 void operator= (const T& newval) {MultiTypeBlockVector_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
310
314 void operator+= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::add(*this,newv);}
315
319 void operator-= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::sub(*this,newv);}
320
321 void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const int>::mul(*this,w);}
322 void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const float>::mul(*this,w);}
323 void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const double>::mul(*this,w);}
324
325 field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::mul(*this,newv);}
326 field_type dot (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::dot(*this,newv);}
327
331 typename FieldTraits<field_type>::real_type two_norm2() const {return MultiTypeBlockVector_Norm<mpl::size<type>::value,type>::result(*this);}
332
336 typename FieldTraits<field_type>::real_type two_norm() const {return sqrt(this->two_norm2());}
337
341 template<typename Ta>
342 void axpy (const Ta& a, const type& y) {
343 MultiTypeBlockVector_AXPY<mpl::size<type>::value,type,Ta>::axpy(*this,a,y);
344 }
345
346 };
347
348
349
355 template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
356 std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9>& v) {
357 MultiTypeBlockVector_Print<0,mpl::size<MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value,MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(v);
358 return s;
359 }
360
361
362
363} // end namespace
364
365#endif // end HAVE_BOOST_FUSION
366#endif // end HAVE_DUNE_BOOST
367
368#endif
Provides the functions dot(a,b) := and dotT(a,b) := .
Type traits to determine the type of reals (when working with complex numbers)
std::ostream & operator<<(std::ostream &s, const array< T, N > &e)
Output operator for array.
Definition: array.hh:26
enable_if<!IsVector< A >::value &&!is_same< typenameFieldTraits< A >::field_type, typenameFieldTraits< A >::real_type >::value, typenamePromotionTraits< A, B >::PromotedType >::type dot(const A &a, const B &b)
computes the dot product for fundamental data types according to Petsc's VectDot function: dot(a,...
Definition: dotproduct.hh:44
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignment.hh:10
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)