Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
tensor4.cc
Go to the documentation of this file.
1
21
22#include "rheolef/tensor4.h"
23namespace rheolef {
24
25// identity tensor:
26// tensor y = ddot(Id,x) ==> x==y
27template<class T>
28tensor4_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
37template<class T>
40{
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}
50template<class T>
51tensor_basic<T>
53{
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}
63template<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}
101template<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}
112template<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}
123template<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}
136template<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}
149template<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}
160template<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}
171template<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}
182template <class T>
183T
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
196template<class T>
197std::ostream&
198tensor4_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) \
219template class tensor4_basic<T>; \
220template tensor_basic<T> ddot (const tensor4_basic<T>&, const tensor_basic<T>&); \
221template tensor_basic<T> ddot (const tensor_basic<T>&, const tensor4_basic<T>&); \
222template T norm2 (const tensor4_basic<T>&);
223
225
226}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
std::ostream & put(std::ostream &out, size_type d=3) const
Definition tensor4.cc:198
tensor4_basic< T > operator+(const tensor4_basic< T > &b) const
Definition tensor4.cc:125
tensor4_basic< T > operator-(const tensor4_basic< T > &b) const
Definition tensor4.cc:138
tensor4_basic< T > & operator-=(const tensor4_basic< T > &)
Definition tensor4.cc:162
tensor4_basic< T > & operator+=(const tensor4_basic< T > &)
Definition tensor4.cc:151
tensor4_basic< T > & operator*=(const T &k)
Definition tensor4.cc:173
tensor4_basic< T > & operator=(const tensor4_basic< T > &a)
Definition tensor4.cc:114
static tensor4_basic< T > eye(size_type d=3)
Definition tensor4.cc:29
T & operator()(size_type i, size_type j, size_type k, size_type l)
Definition tensor4.h:194
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
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
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition vec.h:379
t operator()(const t &a, const t &b)
Definition space.cc:386
STL namespace.