Rheolef  7.1
an efficient C++ finite element environment
tensor4.cc
Go to the documentation of this file.
1 
22 #include "rheolef/tensor4.h"
23 namespace rheolef {
24 
25 // identity tensor:
26 // tensor y = ddot(Id,x) ==> x==y
27 template<class T>
28 tensor4_basic<T>
30 {
32  for (size_t i = 0; i < d; i++)
33  for (size_t j = 0; j < d; j++) a(i,j,i,j) = 1;
34  return a;
35 }
36 // algebra
37 template<class T>
40 {
41  typedef typename tensor4_basic<T>::size_type size_type;
43  for (size_type i = 0; i < 3; i++)
44  for (size_type j = 0; j < 3; j++)
45  for (size_type k = 0; k < 3; k++)
46  for (size_type l = 0; l < 3; l++)
47  y(i,j) += a(i,j,k,l)*x(k,l);
48  return y;
49 }
50 template<class T>
51 tensor_basic<T>
53 {
54  typedef typename tensor4_basic<T>::size_type size_type;
56  for (size_type i = 0; i < 3; i++)
57  for (size_type j = 0; j < 3; j++)
58  for (size_type k = 0; k < 3; k++)
59  for (size_type l = 0; l < 3; l++)
60  x(k,l) += y(i,j)*a(i,j,k,l);
61  return x;
62 }
63 template<class T>
65  const std::initializer_list<std::initializer_list<
66  std::initializer_list<std::initializer_list<T> > > >& il)
67  : _x (tensor_basic<T>(T()))
68 {
69  typedef std::initializer_list<T> L1;
70  typedef std::initializer_list<L1> L2;
71  typedef std::initializer_list<L2> L3;
72  typedef std::initializer_list<L3> L4;
73  typedef typename std::initializer_list<L4>::const_iterator const_iterator;
74  typedef typename std::initializer_list<L3>::const_iterator const_iterator_row3;
75  typedef typename std::initializer_list<L2>::const_iterator const_iterator_row2;
76  typedef typename std::initializer_list<L1>::const_iterator const_iterator_row1;
77  typedef typename std::initializer_list<T>::const_iterator const_iterator_elt;
78  this->operator= (T());
79  check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
80  size_type i = 0;
81  for (const_iterator_row3 iter = il.begin(); iter != il.end(); ++iter, ++i) {
82  const L3& row3 = *iter;
83  check_macro (row3.size() <= 3, "unexpected initializer list size=" << row3.size() << " > 3");
84  size_type j = 0;
85  for (const_iterator_row2 jter = row3.begin(); jter != row3.end(); ++jter, ++j) {
86  const L2& row2 = *jter;
87  check_macro (row2.size() <= 3, "unexpected initializer list size=" << row2.size() << " > 3");
88  size_type k = 0;
89  for (const_iterator_row1 kter = row2.begin(); kter != row2.end(); ++kter, ++k) {
90  const L1& row1 = *kter;
91  check_macro (row1.size() <= 3, "unexpected initializer list size=" << row1.size() << " > 3");
92  size_type l = 0;
93  for (const_iterator_elt lter = row1.begin(); lter != row1.end(); ++lter, ++l) {
94  const T& elt = *lter;
95  operator() (i,j,k,l) = elt;
96  }
97  }
98  }
99  }
100 }
101 template<class T>
104 {
105  for (size_type i = 0; i < 3; i++)
106  for (size_type j = 0; j < 3; j++)
107  for (size_type k = 0; k < 3; k++)
108  for (size_type l = 0; l < 3; l++)
109  operator() (i,j,k,l) = val;
110  return *this;
111 }
112 template<class T>
115 {
116  for (size_type i = 0; i < 3; i++)
117  for (size_type j = 0; j < 3; j++)
118  for (size_type k = 0; k < 3; k++)
119  for (size_type l = 0; l < 3; l++)
120  operator() (i,j,k,l) = a(i,j,k,l);
121  return *this;
122 }
123 template<class T>
126 {
127  typedef typename tensor4_basic<T>::size_type size_type;
129  for (size_type i = 0; i < 3; i++)
130  for (size_type j = 0; j < 3; j++)
131  for (size_type k = 0; k < 3; k++)
132  for (size_type l = 0; l < 3; l++)
133  c(i,j,k,l) = operator() (i,j,k,l) + b(i,j,k,l);
134  return c;
135 }
136 template<class T>
139 {
140  typedef typename tensor4_basic<T>::size_type size_type;
142  for (size_type i = 0; i < 3; i++)
143  for (size_type j = 0; j < 3; j++)
144  for (size_type k = 0; k < 3; k++)
145  for (size_type l = 0; l < 3; l++)
146  c(i,j,k,l) = operator() (i,j,k,l) - b(i,j,k,l);
147  return c;
148 }
149 template<class T>
152 {
153  for (size_type i = 0; i < 3; i++)
154  for (size_type j = 0; j < 3; j++)
155  for (size_type k = 0; k < 3; k++)
156  for (size_type l = 0; l < 3; l++)
157  operator() (i,j,k,l) += b(i,j,k,l);
158  return *this;
159 }
160 template<class T>
163 {
164  for (size_type i = 0; i < 3; i++)
165  for (size_type j = 0; j < 3; j++)
166  for (size_type k = 0; k < 3; k++)
167  for (size_type l = 0; l < 3; l++)
168  operator() (i,j,k,l) -= b(i,j,k,l);
169  return *this;
170 }
171 template<class T>
174 {
175  for (size_type i = 0; i < 3; i++)
176  for (size_type j = 0; j < 3; j++)
177  for (size_type k = 0; k < 3; k++)
178  for (size_type l = 0; l < 3; l++)
179  operator() (i,j,k,l) *= fact;
180  return *this;
181 }
182 template <class T>
183 T
185 {
186  typedef typename tensor4_basic<T>::size_type size_type;
187  T sum = 0;
188  for (size_type i = 0; i < 3; i++)
189  for (size_type j = 0; j < 3; j++)
190  for (size_type k = 0; k < 3; k++)
191  for (size_type l = 0; l < 3; l++)
192  sum += sqr(a(i,j,k,l));
193  return sum;
194 }
195 // output
196 template<class T>
197 std::ostream&
198 tensor4_basic<T>::put (std::ostream& out, size_type d) const
199 {
200  using namespace std;
201  switch (d) {
202  case 0 : return out;
203  case 1 : return out << _x(0,0)(0,0);
204  case 2 : return out << "[[" << _x(0,0)(0,0) << ", " << _x(0,0)(0,1) << ";" << endl
205  << " " << _x(0,0)(1,0) << ", " << _x(0,0)(1,1) << "]," << endl
206  << " [" << _x(0,1)(0,0) << ", " << _x(0,1)(0,1) << ";" << endl
207  << " " << _x(0,1)(1,0) << ", " << _x(0,1)(1,1) << "];" << endl
208  << " [" << _x(1,0)(0,0) << ", " << _x(1,0)(0,1) << ";" << endl
209  << " " << _x(1,0)(1,0) << ", " << _x(1,0)(1,1) << "]," << endl
210  << " [" << _x(1,1)(0,0) << ", " << _x(1,1)(0,1) << ";" << endl
211  << " " << _x(1,1)(1,0) << ", " << _x(1,1)(1,1) << "]]" << endl;
212  default: fatal_macro("put(d="<<d<<"): not yet, sorry"); return out;
213  }
214 }
215 // ----------------------------------------------------------------------------
216 // instanciation in library
217 // ----------------------------------------------------------------------------
218 #define _RHEOLEF_instanciation(T) \
219 template class tensor4_basic<T>; \
220 template tensor_basic<T> ddot (const tensor4_basic<T>&, const tensor_basic<T>&); \
221 template tensor_basic<T> ddot (const tensor_basic<T>&, const tensor4_basic<T>&); \
222 template T norm2 (const tensor4_basic<T>&);
223 
225 
226 }// namespace rheolef
rheolef::tensor4_basic::tensor4_basic
tensor4_basic()
Definition: tensor4.h:160
rheolef::io::out
Definition: rheostream.h:167
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::tensor4_basic::operator()
T & operator()(size_type i, size_type j, size_type k, size_type l)
Definition: tensor4.h:194
rheolef::tensor4_basic
Definition: tensor4.h:80
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::tensor4_basic::operator*=
tensor4_basic< T > & operator*=(const T &k)
Definition: tensor4.cc:173
rheolef::tensor_basic
Definition: tensor.h:90
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
rheolef::tensor4_basic::put
std::ostream & put(std::ostream &out, size_type d=3) const
Definition: tensor4.cc:198
rheolef::tensor4_basic::operator-
tensor4_basic< T > operator-(const tensor4_basic< T > &b) const
Definition: tensor4.cc:138
rheolef::tensor4_basic::operator-=
tensor4_basic< T > & operator-=(const tensor4_basic< T > &)
Definition: tensor4.cc:162
rheolef::ddot
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
a
Definition: diffusion_isotropic.h:25
rheolef::tensor4_basic::operator=
tensor4_basic< T > & operator=(const tensor4_basic< T > &a)
Definition: tensor4.cc:114
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::tensor4_basic::operator+
tensor4_basic< T > operator+(const tensor4_basic< T > &b) const
Definition: tensor4.cc:125
Float
see the Float page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::tensor4_basic::eye
static tensor4_basic< T > eye(size_type d=3)
Definition: tensor4.cc:29
rheolef::tensor4_basic::operator+=
tensor4_basic< T > & operator+=(const tensor4_basic< T > &)
Definition: tensor4.cc:151
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
rheolef::const_iterator
Definition: field_expr_recursive.h:552
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::tensor4_basic::size_type
size_t size_type
Definition: tensor4.h:83
rheolef::std
Definition: vec_expr_v2.h:391
T
Expr1::float_type T
Definition: field_expr.h:218