← Back to branch summary

~ubuntu-branches/ubuntu/quantal/rheolef/quantal

« back to all changes in this revision

Viewing changes to .pc/ld-as-needed.diff/nfem/lib/tensor.h

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2010-12-14 23:20:11 UTC
  • Revision ID: james.westby@ubuntu.com-20101214232011-i2g8dm43k092cxnw
Tags: 5.91-1ubuntu1
Fix build failure with GCC-4.5.

Show diffs side-by-side

added added

removed removed

Lines of Context:
.pc/ld-as-needed.diff/nfem/lib/tensor.h
 
1
# ifndef _RHEO_TENSOR_H
 
2
# define _RHEO_TENSOR_H
 
3
///
 
4
/// This file is part of Rheolef.
 
5
///
 
6
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
7
///
 
8
/// Rheolef is free software; you can redistribute it and/or modify
 
9
/// it under the terms of the GNU General Public License as published by
 
10
/// the Free Software Foundation; either version 2 of the License, or
 
11
/// (at your option) any later version.
 
12
///
 
13
/// Rheolef is distributed in the hope that it will be useful,
 
14
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
15
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
16
/// GNU General Public License for more details.
 
17
///
 
18
/// You should have received a copy of the GNU General Public License
 
19
/// along with Rheolef; if not, write to the Free Software
 
20
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
21
/// 
 
22
/// =========================================================================
 
23
 
 
24
/*Class:tensor
 
25
NAME: @code{tensor} - a N*N tensor, N=1,2,3
 
26
@cindex tensor
 
27
@clindex tensor
 
28
@clindex point
 
29
@clindex field
 
30
SYNOPSYS:
 
31
@noindent
 
32
The @code{tensor} class defines a 3*3 tensor, as the value of
 
33
a tensorial valued field. Basic algebra with scalars, vectors
 
34
of R^3 (i.e. the @code{point} class) and @code{tensor} objects
 
35
are supported.
 
36
AUTHOR: Pierre.Saramito@imag.fr
 
37
DATE:   9 october 2003
 
38
End:
 
39
*/
 
40
 
 
41
#include "rheolef/point.h"
 
42
namespace rheolef { 
 
43
 
 
44
//<tensor:
 
45
class tensor {
 
46
    public:
 
47
 
 
48
        typedef size_t size_type;
 
49
        typedef Float  element_type;
 
50
 
 
51
// allocators:
 
52
 
 
53
        tensor (const Float& init_val = 0);
 
54
        tensor (Float x[3][3]);
 
55
        tensor (const tensor& a);
 
56
 
 
57
// affectation:
 
58
 
 
59
        tensor& operator = (const tensor& a);
 
60
        tensor& operator = (const Float& val);
 
61
 
 
62
// modifiers:
 
63
 
 
64
        void fill (const Float& init_val);
 
65
        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);
 
68
 
 
69
// accessors:
 
70
 
 
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;
 
75
        size_t nrow() const; // = 3, for template matrix compatibility
 
76
        size_t ncol() const;
 
77
 
 
78
// inputs/outputs:
 
79
 
 
80
        friend std::istream& operator >> (std::istream& in, tensor& a);
 
81
        friend std::ostream& operator << (std::ostream& out, const tensor& a);
 
82
        std::ostream& put (std::ostream& s, size_type d = 3) const;
 
83
 
 
84
// algebra:
 
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
 
 
101
// metric and geometric transformations:
 
102
 
 
103
        friend Float dotdot (const tensor&, const tensor&);
 
104
        friend Float norm2 (const tensor& a) { return dotdot(a,a); }
 
105
        friend Float dist2 (const tensor& a, const tensor& b) { return norm2(a-b); }
 
106
        friend Float norm  (const tensor& a) { return ::sqrt(norm2(a)); }
 
107
        friend Float dist  (const tensor& a, const tensor& b) { return norm(a-b); }
 
108
        Float determinant (size_type d = 3) const;
 
109
        friend Float determinant (const tensor& A, size_t d = 3);
 
110
        friend bool invert_3x3 (const tensor& A, tensor& result);
 
111
 
 
112
// spectral:
 
113
        // eigenvalues & eigenvectors:
 
114
        // a = q*d*q^T
 
115
        // a may be symmetric
 
116
        // where q=(q1,q2,q3) are eigenvectors in rows (othonormal matrix)
 
117
        // and   d=(d1,d2,d3) are eigenvalues, sorted in decreasing order d1 >= d2 >= d3
 
118
        // return d
 
119
        point eig (tensor& q, size_t dim = 3) const;
 
120
        point eig (size_t dim = 3) const;
 
121
 
 
122
        // singular value decomposition:
 
123
        // a = u*s*v^T
 
124
        // a can be unsymmetric
 
125
        // where u=(u1,u2,u3) are left pseudo-eigenvectors in rows (othonormal matrix)
 
126
        //       v=(v1,v2,v3) are right pseudo-eigenvectors in rows (othonormal matrix)
 
127
        // and   s=(s1,s2,s3) are eigenvalues, sorted in decreasing order s1 >= s2 >= s3
 
128
        // return s
 
129
        point svd (tensor& u, tensor& v, size_t dim = 3) const;
 
130
 
 
131
// data:
 
132
        Float _x[3][3];
 
133
// internal:
 
134
        std::istream& get (std::istream&);
 
135
};
 
136
// t = a otimes b
 
137
tensor otimes (const point& a, const point& b, size_t na = 3);
 
138
// t += a otimes b
 
139
void cumul_otimes (tensor& t, const point& a, const point& b, size_t na = 3);
 
140
void cumul_otimes (tensor& t, const point& a, const point& b, size_t na, size_t nb);
 
141
//>tensor:
 
142
 
 
143
// -----------------------------------------------------------------------
 
144
// inlined
 
145
// -----------------------------------------------------------------------
 
146
 
 
147
inline
 
148
void
 
149
tensor::fill (const Float& init_val)
 
150
{ 
 
151
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
 
152
      _x[i][j] = init_val;
 
153
}
 
154
inline
 
155
void
 
156
tensor::reset ()
 
157
{
 
158
    fill (0);
 
159
}
 
160
inline
 
161
tensor::tensor (const Float& init_val)
 
162
{ 
 
163
    fill (init_val);
 
164
}
 
165
inline
 
166
tensor::tensor (Float x[3][3])
 
167
{
 
168
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
 
169
      _x[i][j] = x[i][j];
 
170
}
 
171
inline
 
172
tensor::tensor (const tensor& a)
 
173
{
 
174
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
 
175
      _x[i][j] = a._x[i][j];
 
176
}
 
177
inline
 
178
tensor::tensor&
 
179
tensor::operator = (const tensor& a)
 
180
{
 
181
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
 
182
      _x[i][j] = a._x[i][j];
 
183
    return *this;
 
184
}
 
185
inline
 
186
tensor::tensor&
 
187
tensor::operator = (const Float& val)
 
188
{
 
189
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
 
190
      _x[i][j] = val;
 
191
    return *this;
 
192
}
 
193
inline
 
194
size_t
 
195
tensor::nrow() const
 
196
{
 
197
  return 3;
 
198
}
 
199
inline
 
200
size_t
 
201
tensor::ncol() const
 
202
{
 
203
  return 3;
 
204
}
 
205
inline
 
206
Float&
 
207
tensor::operator()(size_type i, size_type j)  
 
208
{
 
209
    return _x[i%3][j%3];
 
210
}
 
211
inline
 
212
Float
 
213
tensor::operator()(size_type i, size_type j) const
 
214
{
 
215
    return _x[i%3][j%3];
 
216
}
 
217
inline
 
218
std::istream& operator >> (std::istream& in, tensor& a)
 
219
{
 
220
    return a.get (in);
 
221
}
 
222
inline
 
223
std::ostream& operator << (std::ostream& out, const tensor& a)
 
224
{
 
225
    return a.put (out); 
 
226
}
 
227
inline
 
228
tensor
 
229
operator * (const tensor& a, Float k)
 
230
{
 
231
    return k*a;
 
232
}
 
233
inline
 
234
tensor
 
235
operator / (const tensor& a, Float k)
 
236
{
 
237
    return (1/k)*a;
 
238
}
 
239
inline
 
240
point
 
241
tensor::trans_mult (const point& x) const
 
242
{
 
243
    return x*(*this);
 
244
}
 
245
inline
 
246
void
 
247
cumul_otimes (tensor& t, const point& a, const point& b, size_t n)
 
248
{
 
249
    cumul_otimes (t, a, b, n, n);
 
250
}
 
251
inline
 
252
tensor
 
253
otimes (const point& a, const point& b, size_t n)
 
254
{
 
255
    tensor t;
 
256
    cumul_otimes (t, a, b, n, n);
 
257
    return t;
 
258
}
 
259
inline
 
260
Float
 
261
determinant (const tensor& A, size_t d)
 
262
{ 
 
263
    return A.determinant (d);
 
264
}
 
265
inline
 
266
tensor
 
267
diag (const point& d)
 
268
{
 
269
  tensor a;
 
270
  a(0,0) = d[0];
 
271
  a(1,1) = d[1];
 
272
  a(2,2) = d[2];
 
273
  return a;
 
274
}
 
275
inline
 
276
void
 
277
tensor::set_column (const point& c, size_t j, size_t d)
 
278
{
 
279
  for (size_t i = 0; i < d; i++)
 
280
    operator()(i,j) = c[i];
 
281
}
 
282
inline
 
283
void
 
284
tensor::set_row (const point& r, size_t i, size_t d)
 
285
{
 
286
  for (size_t j = 0; j < d; j++)
 
287
    operator()(i,j) = r[j];
 
288
}
 
289
}// namespace rheolef
 
290
# endif /* _RHEO_TENSOR_H */