← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to nfem/geo_element/tensor.h

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito, Pierre Saramito, Sylvestre Ledru
  • Date: 2012-05-14 14:02:09 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20120514140209-dzbdlidkotyflf9e
Tags: 6.1-1
http://bugs.debian.org/671996
[ Pierre Saramito ]
* New upstream release 6.1 (minor changes):
  - support arbitrarily polynomial order Pk
  - source code supports g++-4.7 (closes: #671996)

[ Sylvestre Ledru ]
* update of the watch file

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/tensor.h
42
42
namespace rheolef {
43
43
 
44
44
//<tensor:
45
 
class tensor {
 
45
template<class T>
 
46
class tensor_basic {
46
47
    public:
47
48
 
48
49
        typedef size_t size_type;
49
 
        typedef Float  element_type;
 
50
        typedef T  element_type;
50
51
 
51
52
// allocators:
52
53
 
53
 
        tensor (const Float& init_val = 0);
54
 
        tensor (Float x[3][3]);
55
 
        tensor (const tensor& a);
 
54
        tensor_basic (const T& init_val = 0);
 
55
        tensor_basic (T x[3][3]);
 
56
        tensor_basic (const tensor_basic<T>& a);
56
57
 
57
58
// affectation:
58
59
 
59
 
        tensor& operator = (const tensor& a);
60
 
        tensor& operator = (const Float& val);
 
60
        tensor_basic<T>& operator = (const tensor_basic<T>& a);
 
61
        tensor_basic<T>& operator = (const T& val);
61
62
 
62
63
// modifiers:
63
64
 
64
 
        void fill (const Float& init_val);
 
65
        void fill (const T& init_val);
65
66
        void reset ();
66
 
        void set_row    (const point& r, size_t i, size_t d = 3);
67
 
        void set_column (const point& c, size_t j, size_t d = 3);
 
67
        void set_row    (const point_basic<T>& r, size_t i, size_t d = 3);
 
68
        void set_column (const point_basic<T>& c, size_t j, size_t d = 3);
68
69
 
69
70
// accessors:
70
71
 
71
 
        Float& operator()(size_type i, size_type j);
72
 
        Float  operator()(size_type i, size_type j) const;
73
 
        point  row(size_type i) const;
74
 
        point  col(size_type i) const;
 
72
        T& operator()(size_type i, size_type j);
 
73
        T  operator()(size_type i, size_type j) const;
 
74
        point_basic<T>  row(size_type i) const;
 
75
        point_basic<T>  col(size_type i) const;
75
76
        size_t nrow() const; // = 3, for template matrix compatibility
76
77
        size_t ncol() const;
77
78
 
78
79
// inputs/outputs:
79
80
 
80
 
        friend std::istream& operator >> (std::istream& in, tensor& a);
81
 
        friend std::ostream& operator << (std::ostream& out, const tensor& a);
82
81
        std::ostream& put (std::ostream& s, size_type d = 3) const;
 
82
        std::istream& get (std::istream&);
83
83
 
84
84
// algebra:
85
85
 
86
 
        friend bool  operator == (const tensor&, const tensor&);
87
 
        friend tensor operator - (const tensor&);
88
 
        friend tensor operator + (const tensor&, const tensor&);
89
 
        friend tensor operator - (const tensor&, const tensor&);
90
 
        friend tensor operator * (Float, const tensor&);
91
 
        friend tensor operator * (const tensor& a, Float k);
92
 
        friend tensor operator / (const tensor& a, Float k);
93
 
        friend point  operator * (const tensor& a, const point& x);
94
 
        friend point  operator * (const point& yt, const tensor& a);
95
 
               point  trans_mult (const point& x) const;
96
 
        friend tensor trans      (const tensor& a, size_type d = 3);
97
 
        friend tensor operator * (const tensor& a, const tensor& b);
98
 
        friend tensor inv        (const tensor& a, size_type d = 3);
99
 
        friend tensor diag (const point& d);
100
 
        friend tensor identity (size_t d=3);
101
 
        friend tensor dyadic (const point& u, const point& v, size_t d=3);
 
86
        bool operator== (const tensor_basic<T>&) const;
 
87
        bool operator!= (const tensor_basic<T>& b) const { return ! operator== (b); }
 
88
        template <class U>
 
89
        friend tensor_basic<U> operator- (const tensor_basic<U>&);
 
90
        template <class U>
 
91
        friend tensor_basic<U> operator+ (const tensor_basic<U>&, const tensor_basic<U>&);
 
92
        template <class U>
 
93
        friend tensor_basic<U> operator- (const tensor_basic<U>&, const tensor_basic<U>&);
 
94
        template <class U>
 
95
        friend tensor_basic<U> operator* (int k, const tensor_basic<U>& a);
 
96
        template <class U>
 
97
        friend tensor_basic<U> operator* (const U& k, const tensor_basic<U>& a);
 
98
        template <class U>
 
99
        friend tensor_basic<U> operator* (const tensor_basic<U>& a, int k);
 
100
        template <class U>
 
101
        friend tensor_basic<U> operator* (const tensor_basic<U>& a, const U& k);
 
102
        template <class U>
 
103
        friend tensor_basic<U> operator/ (const tensor_basic<U>& a, int k);
 
104
        template <class U>
 
105
        friend tensor_basic<U> operator/ (const tensor_basic<U>& a, const U& k);
 
106
        template <class U>
 
107
        friend point_basic<U>  operator* (const tensor_basic<U>&, const point_basic<U>&);
 
108
        template <class U>
 
109
        friend point_basic<U>  operator* (const point_basic<U>& yt, const tensor_basic<U>& a);
 
110
               point_basic<T>  trans_mult (const point_basic<T>& x) const;
 
111
        template <class U>
 
112
        friend tensor_basic<U> trans      (const tensor_basic<U>& a, size_t d = 3);
 
113
        template <class U>
 
114
        friend tensor_basic<U> operator* (const tensor_basic<U>& a, const tensor_basic<U>& b);
 
115
        template <class U>
 
116
        friend void prod (const tensor_basic<U>& a, const tensor_basic<U>& b, tensor_basic<U>& result,
 
117
                size_t di=3, size_t dj=3, size_t dk=3);
 
118
        template <class U>
 
119
        friend tensor_basic<U> inv        (const tensor_basic<U>& a, size_t d = 3);
 
120
        template <class U>
 
121
        friend tensor_basic<U> diag (const point_basic<U>& d);
 
122
        template <class U>
 
123
        friend tensor_basic<U> identity (size_t d=3);
 
124
        template <class U>
 
125
        friend tensor_basic<U> dyadic (const point_basic<U>& u, const point_basic<U>& v, size_t d=3);
102
126
 
103
127
// metric and geometric transformations:
104
128
 
105
 
        friend Float dotdot (const tensor&, const tensor&);
106
 
        friend Float norm2 (const tensor& a) { return dotdot(a,a); }
107
 
        friend Float dist2 (const tensor& a, const tensor& b) { return norm2(a-b); }
108
 
        friend Float norm  (const tensor& a) { return ::sqrt(norm2(a)); }
109
 
        friend Float dist  (const tensor& a, const tensor& b) { return norm(a-b); }
110
 
        Float determinant (size_type d = 3) const;
111
 
        friend Float determinant (const tensor& A, size_t d = 3);
112
 
        friend bool invert_3x3 (const tensor& A, tensor& result);
 
129
        template <class U>
 
130
        friend U dotdot (const tensor_basic<U>&, const tensor_basic<U>&);
 
131
        template <class U>
 
132
        friend U norm2 (const tensor_basic<U>& a) { return dotdot(a,a); }
 
133
        template <class U>
 
134
        friend U dist2 (const tensor_basic<U>& a, const tensor_basic<U>& b) { return norm2(a-b); }
 
135
        template <class U>
 
136
        friend U norm  (const tensor_basic<U>& a) { return ::sqrt(norm2(a)); }
 
137
        template <class U>
 
138
        friend U dist  (const tensor_basic<U>& a, const tensor_basic<U>& b) { return norm(a-b); }
 
139
        T determinant (size_type d = 3) const;
 
140
        template <class U>
 
141
        friend U determinant (const tensor_basic<U>& A, size_t d = 3);
 
142
        template <class U>
 
143
        friend bool invert_3x3 (const tensor_basic<U>& A, tensor_basic<U>& result);
113
144
 
114
145
// spectral:
115
146
        // eigenvalues & eigenvectors:
118
149
        // where q=(q1,q2,q3) are eigenvectors in rows (othonormal matrix)
119
150
        // and   d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3
120
151
        // return d
121
 
        point eig (tensor& q, size_t dim = 3) const;
122
 
        point eig (size_t dim = 3) const;
 
152
        point_basic<T> eig (tensor_basic<T>& q, size_t dim = 3) const;
 
153
        point_basic<T> eig (size_t dim = 3) const;
123
154
 
124
155
        // singular value decomposition:
125
156
        // a = u*s*v^T
128
159
        //       v=(v1,v2,v3) are right pseudo-eigenvectors in rows (othonormal matrix)
129
160
        // and   s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3
130
161
        // return s
131
 
        point svd (tensor& u, tensor& v, size_t dim = 3) const;
 
162
        point_basic<T> svd (tensor_basic<T>& u, tensor_basic<T>& v, size_t dim = 3) const;
132
163
 
133
164
// data:
134
 
        Float _x[3][3];
135
 
// internal:
136
 
        std::istream& get (std::istream&);
 
165
        T _x[3][3];
137
166
};
 
167
typedef tensor_basic<Float> tensor;
 
168
 
 
169
// inputs/outputs:
 
170
template<class T>
 
171
inline
 
172
std::istream& operator>> (std::istream& in, tensor_basic<T>& a)
 
173
{
 
174
    return a.get (in);
 
175
}
 
176
template<class T>
 
177
inline
 
178
std::ostream& operator<< (std::ostream& out, const tensor_basic<T>& a)
 
179
{
 
180
    return a.put (out); 
 
181
}
138
182
// t = a otimes b
139
 
tensor otimes (const point& a, const point& b, size_t na = 3);
 
183
template<class T>
 
184
tensor_basic<T> otimes (const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
140
185
// t += a otimes b
141
 
void cumul_otimes (tensor& t, const point& a, const point& b, size_t na = 3);
142
 
void cumul_otimes (tensor& t, const point& a, const point& b, size_t na, size_t nb);
 
186
template<class T>
 
187
void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na = 3);
 
188
template<class T>
 
189
void cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t na, size_t nb);
143
190
//>tensor:
144
191
 
145
192
// -----------------------------------------------------------------------
146
193
// inlined
147
194
// -----------------------------------------------------------------------
148
195
 
 
196
template<class T>
149
197
inline
150
198
void
151
 
tensor::fill (const Float& init_val)
 
199
tensor_basic<T>::fill (const T& init_val)
152
200
{ 
153
201
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
154
202
      _x[i][j] = init_val;
155
203
}
 
204
template<class T>
156
205
inline
157
206
void
158
 
tensor::reset ()
 
207
tensor_basic<T>::reset ()
159
208
{
160
209
    fill (0);
161
210
}
 
211
template<class T>
162
212
inline
163
 
tensor::tensor (const Float& init_val)
 
213
tensor_basic<T>::tensor_basic (const T& init_val)
164
214
{ 
165
215
    fill (init_val);
166
216
}
 
217
template<class T>
167
218
inline
168
 
tensor::tensor (Float x[3][3])
 
219
tensor_basic<T>::tensor_basic (T x[3][3])
169
220
{
170
221
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
171
222
      _x[i][j] = x[i][j];
172
223
}
 
224
template<class T>
173
225
inline
174
 
tensor::tensor (const tensor& a)
 
226
tensor_basic<T>::tensor_basic (const tensor_basic<T>& a)
175
227
{
176
228
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
177
229
      _x[i][j] = a._x[i][j];
178
230
}
 
231
template<class T>
179
232
inline
180
 
tensor&
181
 
tensor::operator = (const tensor& a)
 
233
tensor_basic<T>&
 
234
tensor_basic<T>::operator= (const tensor_basic<T>& a)
182
235
{
183
236
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
184
237
      _x[i][j] = a._x[i][j];
185
238
    return *this;
186
239
}
 
240
template<class T>
187
241
inline
188
 
tensor&
189
 
tensor::operator = (const Float& val)
 
242
tensor_basic<T>&
 
243
tensor_basic<T>::operator= (const T& val)
190
244
{
191
245
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
192
246
      _x[i][j] = val;
193
247
    return *this;
194
248
}
195
 
inline
196
 
size_t
197
 
tensor::nrow() const
198
 
{
199
 
  return 3;
200
 
}
201
 
inline
202
 
size_t
203
 
tensor::ncol() const
204
 
{
205
 
  return 3;
206
 
}
207
 
inline
208
 
Float&
209
 
tensor::operator()(size_type i, size_type j)  
210
 
{
211
 
    return _x[i%3][j%3];
212
 
}
213
 
inline
214
 
Float
215
 
tensor::operator()(size_type i, size_type j) const
216
 
{
217
 
    return _x[i%3][j%3];
218
 
}
219
 
inline
220
 
std::istream& operator >> (std::istream& in, tensor& a)
221
 
{
222
 
    return a.get (in);
223
 
}
224
 
inline
225
 
std::ostream& operator << (std::ostream& out, const tensor& a)
226
 
{
227
 
    return a.put (out); 
228
 
}
229
 
inline
230
 
tensor
231
 
operator * (const tensor& a, Float k)
 
249
template<class T>
 
250
inline
 
251
size_t
 
252
tensor_basic<T>::nrow() const
 
253
{
 
254
  return 3;
 
255
}
 
256
template<class T>
 
257
inline
 
258
size_t
 
259
tensor_basic<T>::ncol() const
 
260
{
 
261
  return 3;
 
262
}
 
263
template<class T>
 
264
inline
 
265
T&
 
266
tensor_basic<T>::operator()(size_type i, size_type j)  
 
267
{
 
268
    return _x[i%3][j%3];
 
269
}
 
270
template<class T>
 
271
inline
 
272
T
 
273
tensor_basic<T>::operator()(size_type i, size_type j) const
 
274
{
 
275
    return _x[i%3][j%3];
 
276
}
 
277
template<class T>
 
278
inline
 
279
tensor_basic<T>
 
280
operator* (int k, const tensor_basic<T>& a)
 
281
{
 
282
    return T(k)*a;
 
283
}
 
284
template<class T>
 
285
inline
 
286
tensor_basic<T>
 
287
operator* (const tensor_basic<T>& a, int k)
 
288
{
 
289
    return T(k)*a;
 
290
}
 
291
template<class T>
 
292
inline
 
293
tensor_basic<T>
 
294
operator* (const tensor_basic<T>& a, const T& k)
232
295
{
233
296
    return k*a;
234
297
}
235
 
inline
236
 
tensor
237
 
operator / (const tensor& a, Float k)
 
298
template<class T>
 
299
inline
 
300
tensor_basic<T>
 
301
operator/ (const tensor_basic<T>& a, int k)
 
302
{
 
303
    return (1.0/T(k))*a;
 
304
}
 
305
template<class T>
 
306
inline
 
307
tensor_basic<T>
 
308
operator/ (const tensor_basic<T>& a, const T& k)
238
309
{
239
310
    return (1/k)*a;
240
311
}
 
312
template<class T>
241
313
inline
242
 
point
243
 
tensor::trans_mult (const point& x) const
 
314
point_basic<T>
 
315
tensor_basic<T>::trans_mult (const point_basic<T>& x) const
244
316
{
245
317
    return x*(*this);
246
318
}
 
319
template<class T>
247
320
inline
248
321
void
249
 
cumul_otimes (tensor& t, const point& a, const point& b, size_t n)
 
322
cumul_otimes (tensor_basic<T>& t, const point_basic<T>& a, const point_basic<T>& b, size_t n)
250
323
{
251
324
    cumul_otimes (t, a, b, n, n);
252
325
}
 
326
template<class T>
253
327
inline
254
 
tensor
255
 
otimes (const point& a, const point& b, size_t n)
 
328
tensor_basic<T>
 
329
otimes (const point_basic<T>& a, const point_basic<T>& b, size_t n)
256
330
{
257
 
    tensor t;
 
331
    tensor_basic<T> t;
258
332
    cumul_otimes (t, a, b, n, n);
259
333
    return t;
260
334
}
 
335
template<class T>
261
336
inline
262
 
Float
263
 
determinant (const tensor& A, size_t d)
 
337
T
 
338
determinant (const tensor_basic<T>& A, size_t d)
264
339
{ 
265
340
    return A.determinant (d);
266
341
}
 
342
template<class T>
267
343
inline
268
 
tensor
269
 
diag (const point& d)
 
344
tensor_basic<T>
 
345
diag (const point_basic<T>& d)
270
346
{
271
 
  tensor a;
 
347
  tensor_basic<T> a;
272
348
  a(0,0) = d[0];
273
349
  a(1,1) = d[1];
274
350
  a(2,2) = d[2];
275
351
  return a;
276
352
}
 
353
template<class T>
277
354
inline
278
 
tensor
 
355
tensor_basic<T>
279
356
identity (size_t d)
280
357
{
281
 
  tensor a;
 
358
  tensor_basic<T> a;
282
359
  for (size_t i = 0; i < d; i++) a(i,i) = 1;
283
360
  return a;
284
361
}
 
362
template<class T>
285
363
inline
286
 
tensor 
287
 
dyadic (const point& u, const point& v, size_t d)
 
364
tensor_basic<T> 
 
365
dyadic (const point_basic<T>& u, const point_basic<T>& v, size_t d)
288
366
{
289
 
  tensor a;
 
367
  tensor_basic<T> a;
290
368
  for (size_t i = 0; i < d; i++) 
291
369
    for (size_t j = 0; j < d; j++) 
292
370
      a(i,j) = u[i]*v[j];
293
371
  return a;
294
372
}
 
373
template<class T>
295
374
inline
296
375
void
297
 
tensor::set_column (const point& c, size_t j, size_t d)
 
376
tensor_basic<T>::set_column (const point_basic<T>& c, size_t j, size_t d)
298
377
{
299
378
  for (size_t i = 0; i < d; i++)
300
379
    operator()(i,j) = c[i];
301
380
}
 
381
template<class T>
302
382
inline
303
383
void
304
 
tensor::set_row (const point& r, size_t i, size_t d)
 
384
tensor_basic<T>::set_row (const point_basic<T>& r, size_t i, size_t d)
305
385
{
306
386
  for (size_t j = 0; j < d; j++)
307
387
    operator()(i,j) = r[j];