← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/basis/tensor.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/basis/tensor.cc
1
 
///
2
 
/// This file is part of Rheolef.
3
 
///
4
 
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
5
 
///
6
 
/// Rheolef is free software; you can redistribute it and/or modify
7
 
/// it under the terms of the GNU General Public License as published by
8
 
/// the Free Software Foundation; either version 2 of the License, or
9
 
/// (at your option) any later version.
10
 
///
11
 
/// Rheolef is distributed in the hope that it will be useful,
12
 
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
 
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14
 
/// GNU General Public License for more details.
15
 
///
16
 
/// You should have received a copy of the GNU General Public License
17
 
/// along with Rheolef; if not, write to the Free Software
18
 
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
19
 
///
20
 
/// =========================================================================
21
 
//
22
 
//    point : set of coordinates
23
 
//
24
 
//    author: Pierre.Saramito@imag.fr
25
 
//
26
 
#include "rheolef/tensor.h"
27
 
#include "rheolef/tensor-eig3.h"
28
 
#include "rheolef/tensor-svd-algo.h"
29
 
using namespace rheolef;
30
 
using namespace std;
31
 
namespace rheolef { 
32
 
 
33
 
 
34
 
// output
35
 
 
36
 
ostream& 
37
 
tensor::put (ostream& out, size_type d) const
38
 
{
39
 
    switch (d) {
40
 
    case 0 : return out;
41
 
    case 1 : return out << _x[0][0];
42
 
    case 2 : return out << "[" << _x[0][0] << ", " << _x[0][1] << ";\n" 
43
 
                        << " " << _x[1][0] << ", " << _x[1][1] << "]";
44
 
    default: return out << "[" << _x[0][0] << ", " << _x[0][1] << ", " << _x[0][2] << ";\n" 
45
 
                        << " " << _x[1][0] << ", " << _x[1][1] << ", " << _x[1][2] << ";\n"
46
 
                        << " " << _x[2][0] << ", " << _x[2][1] << ", " << _x[2][2] << "]";
47
 
    }
48
 
}
49
 
istream&
50
 
tensor::get (istream& in)
51
 
{
52
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
53
 
        _x[i][j] = 0;
54
 
    char c;
55
 
    if (!in.good()) return in;
56
 
    in >> ws >> c;
57
 
    // no tensor specifier: 
58
 
    // 1d case: "3.14"
59
 
    if (c != '[') {
60
 
        in.unget();
61
 
        return in >> _x[0][0];
62
 
    }
63
 
    // has a tensor specifier, as:
64
 
    //          [  23, 17; 25, 2]
65
 
    in >> _x[0][0] >> ws >> c;
66
 
    if (c == ']') return in; // 1d case
67
 
    if (c != ',') error_macro ("invalid tensor input: expect `,'");
68
 
    in >> _x[0][1];
69
 
    in >> ws >> c;
70
 
    size_type d = 3;
71
 
    if (c == ',') d = 3;
72
 
    else if (c == ';') d = 2;
73
 
    else error_macro ("invalid tensor input: expect `,' or `;'");
74
 
    if (d == 3) {
75
 
        in >> _x[0][2] >> ws >> c;
76
 
        if (c != ';') error_macro ("invalid tensor input: expect `;'");
77
 
    }
78
 
    in >> _x[1][0] >> ws >> c;
79
 
    if (c != ',') error_macro ("invalid tensor input: expect `,'");
80
 
    in >> _x[1][1] >> ws >> c;
81
 
    if (d == 2) {
82
 
        if (c != ']') error_macro ("invalid tensor input: expect `]'");
83
 
        return in;
84
 
    }
85
 
    // d == 3
86
 
    if (c != ',') error_macro ("invalid tensor input: expect `,'");
87
 
    in >> _x[1][2] >> ws >> c;
88
 
    if (c != ';') error_macro ("invalid tensor input: expect `;'");
89
 
    in >> _x[2][0] >> ws >> c;
90
 
    if (c != ',') error_macro ("invalid tensor input: expect `,'");
91
 
    in >> _x[2][1] >> ws >> c;
92
 
    if (c != ',') error_macro ("invalid tensor input: expect `,'");
93
 
    in >> _x[2][2] >> ws >> c;
94
 
    if (c != ']') error_macro ("invalid tensor input: expect `]'");
95
 
    return in;
96
 
}
97
 
bool operator == (const tensor& a, const tensor& b)
98
 
{
99
 
    typedef tensor::size_type size_type;
100
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
101
 
        if (a._x[i][j] != b._x[i][j]) return false;
102
 
    return true;
103
 
}
104
 
tensor operator - (const tensor& a)
105
 
{
106
 
    tensor b;
107
 
    typedef tensor::size_type size_type;
108
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
109
 
        b._x[i][j] = - a._x[i][j];
110
 
    return b;
111
 
}
112
 
tensor operator + (const tensor& a, const tensor& b)
113
 
{
114
 
    tensor c;
115
 
    typedef tensor::size_type size_type;
116
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
117
 
        c._x[i][j] = a._x[i][j] + b._x[i][j];
118
 
    return c;
119
 
}
120
 
tensor operator - (const tensor& a, const tensor& b)
121
 
{
122
 
    tensor c;
123
 
    typedef tensor::size_type size_type;
124
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
125
 
        c._x[i][j] = a._x[i][j] - b._x[i][j];
126
 
    return c;
127
 
}
128
 
tensor operator * (Float k, const tensor& a)
129
 
{
130
 
    tensor b;
131
 
    typedef tensor::size_type size_type;
132
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
133
 
        b._x[i][j] = k * a._x[i][j];
134
 
    return b;
135
 
}
136
 
point operator * (const tensor& a, const point& x)
137
 
{
138
 
    point y;
139
 
    typedef tensor::size_type size_type;
140
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
141
 
        y [i] += a._x[i][j] * x[j];
142
 
    return y;
143
 
}
144
 
point operator * (const point& x, const tensor& a)
145
 
{
146
 
    point y;
147
 
    typedef tensor::size_type size_type;
148
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
149
 
        y [i] += a._x[j][i] * x[j];
150
 
    return y;
151
 
}
152
 
tensor
153
 
trans (const tensor& a, tensor::size_type d) {
154
 
    tensor b;
155
 
    typedef tensor::size_type size_type;
156
 
    for (size_type i = 0; i < d; i++) for (size_type j = 0; j < d; j++)
157
 
        b._x[i][j] = a._x[j][i];
158
 
    return b;
159
 
}
160
 
tensor
161
 
inv (const tensor& a, tensor::size_type d) {
162
 
  tensor b;
163
 
  typedef tensor::size_type size_type;
164
 
  switch (d) {
165
 
    case 0: 
166
 
      break;
167
 
    case 1: 
168
 
      b(0,0) = 1/a(0,0);
169
 
      break;
170
 
    case 2: {
171
 
      Float det = a.determinant(2);
172
 
      b(0,0) =   a(1,1)/det;
173
 
      b(1,1) =   a(0,0)/det;
174
 
      b(0,1) = - a(0,1)/det;
175
 
      b(1,0) = - a(1,0)/det;
176
 
      break;
177
 
    }
178
 
    case 3:
179
 
    default: {
180
 
      Float det = a.determinant(3);
181
 
      b(0,0) = (a(1,1)*a(2,2) - a(1,2)*a(2,1))/det;
182
 
      b(0,1) = (a(0,2)*a(2,1) - a(0,1)*a(2,2))/det;
183
 
      b(0,2) = (a(0,1)*a(1,2) - a(0,2)*a(1,1))/det;
184
 
      b(1,0) = (a(1,2)*a(2,0) - a(1,0)*a(2,2))/det;
185
 
      b(1,1) = (a(0,0)*a(2,2) - a(0,2)*a(2,0))/det;
186
 
      b(1,2) = (a(0,2)*a(1,0) - a(0,0)*a(1,2))/det;
187
 
      b(2,0) = (a(1,0)*a(2,1) - a(1,1)*a(2,0))/det;
188
 
      b(2,1) = (a(0,1)*a(2,0) - a(0,0)*a(2,1))/det;
189
 
      b(2,2) = (a(0,0)*a(1,1) - a(0,1)*a(1,0))/det;
190
 
      break;
191
 
    }
192
 
  }
193
 
  return b;
194
 
}
195
 
tensor operator * (const tensor& a, const tensor& b)
196
 
{
197
 
    tensor c;
198
 
    typedef tensor::size_type size_type;
199
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
200
 
      for (size_type k = 0; k < 3; k++)
201
 
        c._x[i][j] += a._x[i][k] * b._x[k][j];
202
 
    return c;
203
 
}
204
 
Float dotdot (const tensor& a, const tensor & b) {
205
 
    Float r = 0;
206
 
    typedef tensor::size_type size_type;
207
 
    for (size_type i = 0; i < 3; i++) for (size_type j = 0; j < 3; j++)
208
 
        r += a._x[i][j] * b._x[i][j];
209
 
    return r;
210
 
}
211
 
Float
212
 
tensor::determinant (size_type d) const
213
 
{
214
 
    switch (d) {
215
 
      case 1 : return _x[0][0];
216
 
      case 2 : return _x[0][0]*_x[1][1] - _x[0][1]*_x[1][0];
217
 
      case 3 : return _x[0][0]*(_x[1][1]*_x[2][2] - _x[1][2]*_x[2][1])
218
 
                    - _x[0][1]*(_x[1][0]*_x[2][2] - _x[1][2]*_x[2][0])
219
 
                    + _x[0][2]*(_x[1][0]*_x[2][1] - _x[1][1]*_x[2][0]);
220
 
      default: {
221
 
        error_macro ("determinant: unexpected dimension " << d);
222
 
      }
223
 
    }
224
 
}
225
 
// t += a otimes b
226
 
void
227
 
cumul_otimes (tensor& t, const point& a, const point& b, size_t na, size_t nb)
228
 
{
229
 
  for (size_t i = 0; i < na; i++)
230
 
    for (size_t j = 0; j < nb; j++)
231
 
      t(i,j) += a[i] * b[j];
232
 
}
233
 
point
234
 
tensor::row(size_type i) const
235
 
{
236
 
    point r;
237
 
    r[0] = _x[i][0];
238
 
    r[1] = _x[i][1];
239
 
    r[2] = _x[i][2];
240
 
    return r;
241
 
}
242
 
point
243
 
tensor::col(size_type j) const
244
 
{
245
 
    point c;
246
 
    c[0] = _x[0][j];
247
 
    c[1] = _x[1][j];
248
 
    c[2] = _x[2][j];
249
 
    return c;
250
 
}
251
 
bool
252
 
invert_3x3 (const tensor& A, tensor& result) {
253
 
  Float det = determinant(A);
254
 
  if (1+det == 1) return false;
255
 
  Float invdet = 1.0/det;
256
 
  result(0,0) =  (A(1,1)*A(2,2)-A(2,1)*A(1,2))*invdet;
257
 
  result(1,0) = -(A(0,1)*A(2,2)-A(0,2)*A(2,1))*invdet;
258
 
  result(2,0) =  (A(0,1)*A(1,2)-A(0,2)*A(1,1))*invdet;
259
 
  result(0,1) = -(A(1,0)*A(2,2)-A(1,2)*A(2,0))*invdet;
260
 
  result(1,1) =  (A(0,0)*A(2,2)-A(0,2)*A(2,0))*invdet;
261
 
  result(2,1) = -(A(0,0)*A(1,2)-A(1,0)*A(0,2))*invdet;
262
 
  result(0,2) =  (A(1,0)*A(2,1)-A(2,0)*A(1,1))*invdet;
263
 
  result(1,2) = -(A(0,0)*A(2,1)-A(2,0)*A(0,1))*invdet;
264
 
  result(2,2) =  (A(0,0)*A(1,1)-A(1,0)*A(0,1))*invdet;
265
 
  return true;
266
 
}
267
 
// ------------------------------------------
268
 
// eigenvalues: the 3D case
269
 
// ------------------------------------------
270
 
static point eig3x3 (const tensor& a, tensor& q)
271
 
{
272
 
    Float A[3][3];
273
 
    for (size_t i = 0; i < 3; i++)
274
 
      for (size_t j = 0; j < 3; j++)
275
 
        A[i][j] = a(i,j);
276
 
    Float Q[3][3];
277
 
    Float D[3];
278
 
    eigen_decomposition(A, Q, D);
279
 
    point d;
280
 
    for (size_t i = 0; i < 3; i++) {
281
 
      d[i] = D[i];
282
 
      for (size_t j = 0; j < 3; j++)
283
 
        q(i,j) = Q[i][j];
284
 
    }
285
 
    return d;
286
 
}
287
 
// ------------------------------------------
288
 
// eigenvalues: the 2D case
289
 
// ------------------------------------------
290
 
static point eig2x2 (const tensor& a, tensor& q)
291
 
{
292
 
    // TODO: check if a is symetric !
293
 
    point d (0,0,0);
294
 
    if (a(0,1) == 0) {
295
 
      // a is already diagonal
296
 
      if (a(0,0) >= a(1,1)) {
297
 
        d[0] = a(0,0);
298
 
        d[1] = a(1,1);
299
 
        q(0,0) = q(1,1) = 1;
300
 
        q(0,1) = q(1,0) = 0;
301
 
      } else { // swap column 0 & 1:
302
 
        d[1] = a(0,0);
303
 
        d[0] = a(1,1);
304
 
        q(0,0) = q(1,1) = 0;
305
 
        q(0,1) = q(1,0) = 1;
306
 
      }
307
 
      return d;
308
 
    }
309
 
    // here a(0,1) != 0
310
 
    Float discr = sqr(a(1,1) - a(0,0)) + 4*sqr(a(0,1));
311
 
    Float trace = a(0,0) + a(1,1);
312
 
    d[0] = (trace + ::sqrt(discr))/2; 
313
 
    d[1] = (trace - ::sqrt(discr))/2;
314
 
    Float lost_precision = 1e-6;
315
 
    if (fabs(d[0]) + lost_precision*fabs(d[1]) == fabs(d[0])) { // d[1] == 0 up to machine precision
316
 
        d[1] = 0;
317
 
    }
318
 
    Float c0 = (d[0]-a(0,0))/a(0,1);
319
 
    Float n0 = ::sqrt(1+sqr(c0));
320
 
    q(0,0) =  1/n0;
321
 
    q(1,0) = c0/n0;
322
 
    Float c1 = (d[1]-a(0,0))/a(0,1);
323
 
    Float n1 = ::sqrt(1+sqr(c1));
324
 
    q(0,1) =  1/n1;
325
 
    q(1,1) = c1/n1;
326
 
    return d;
327
 
}
328
 
// ------------------------------------------
329
 
// eigenvalues: general case
330
 
// ------------------------------------------
331
 
point
332
 
tensor::eig (tensor& q, size_t d) const
333
 
{
334
 
  switch (d) {
335
 
    case 1 : {
336
 
      point d;
337
 
      q(0,0) = 1; d[0] = operator()(0,0); return d;
338
 
    }
339
 
    case 2 : return eig2x2 (*this, q);
340
 
    default: return eig3x3 (*this, q);
341
 
  }
342
 
}
343
 
point 
344
 
tensor::eig (size_t d) const
345
 
{
346
 
  tensor q;
347
 
  return eig (q, d);
348
 
}
349
 
// ------------------------------------------
350
 
// svd: the 3D case
351
 
// ------------------------------------------
352
 
point
353
 
tensor::svd (tensor& u, tensor& v, size_t dim) const
354
 
{
355
 
  Float U[3][3];
356
 
  for (size_t i = 0; i < dim; i++)
357
 
    for (size_t j = 0; j < dim; j++)
358
 
      U[i][j] = operator()(i,j);
359
 
  Float V[3][3];
360
 
  Float S[3];
361
 
  svdcmp (U, S, V, dim, dim, Float(0));
362
 
  // Sort eigenvalues and corresponding vectors.
363
 
  for (size_t i = 0; i < dim-1; i++) {
364
 
    size_t k = i;
365
 
    Float p = S[i];
366
 
    for (size_t j = i+1; j < dim; j++) {
367
 
      if (S[j] > p) {
368
 
        k = j;
369
 
        p = S[j];
370
 
      }
371
 
    }
372
 
    if (k != i) {
373
 
      S[k] = S[i];
374
 
      S[i] = p;
375
 
      for (size_t j = 0; j < dim; j++) {
376
 
        p = U[j][i];
377
 
        U[j][i] = U[j][k];
378
 
        U[j][k] = p;
379
 
        p = V[j][i];
380
 
        V[j][i] = V[j][k];
381
 
        V[j][k] = p;
382
 
      }
383
 
    }
384
 
  }
385
 
  // copy to tensor and point
386
 
  point s;
387
 
  for (size_t i = 0; i < dim; i++) {
388
 
    s[i] = S[i];
389
 
    for (size_t j = 0; j < dim; j++) {
390
 
      u(i,j) = U[i][j];
391
 
      v(i,j) = V[i][j];
392
 
    }
393
 
  }
394
 
  return s;
395
 
}
396
 
}// namespace rheolef