← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/lib/piola-algo.h

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/lib/piola-algo.h
 
1
#ifndef _RHEO_PIOLA_ALGO_H
 
2
#define _RHEO_PIOLA_ALGO_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
#include "rheolef/tensor.h"
 
23
#include "rheolef/geo.h"
 
24
// ---------------------------------------------------------------------------
 
25
// piola transform and its inverse on simplex
 
26
// ---------------------------------------------------------------------------
 
27
// piola transformation
 
28
// F_K : hat_K --> K
 
29
//       hat_x --> x = F_K(hat_x)
 
30
// ------------------------------------------
 
31
// TODO: should be merged with
 
32
//    piola.h  = general piola on all elements (non-simplex)
 
33
// this file is actually used only by ZZ posteriori error estimators
 
34
//
 
35
inline
 
36
point
 
37
piola_e (
 
38
  const point& hat_x, 
 
39
  const point& a, 
 
40
  const point& b,
 
41
  bool is_a_vector = false) 
 
42
{
 
43
  return point((!is_a_vector?a[0]:0)
 
44
               + (b[0]-a[0])*hat_x[0]);
 
45
}
 
46
inline
 
47
point
 
48
piola_t (
 
49
  const point& hat_x, 
 
50
  const point& a, 
 
51
  const point& b,
 
52
  const point& c,
 
53
  bool is_a_vector = false)
 
54
{
 
55
  return ((!is_a_vector)?a:point(0,0))
 
56
          + (b-a)*hat_x[0]
 
57
          + (c-a)*hat_x[1];
 
58
}
 
59
inline
 
60
point
 
61
piola_T (
 
62
  const point& hat_x, 
 
63
  const point& a, 
 
64
  const point& b,
 
65
  const point& c,
 
66
  const point& d,
 
67
  bool is_a_vector = false)
 
68
{
 
69
  return ((!is_a_vector)?a:point(0,0))
 
70
          + (b-a)*hat_x[0] 
 
71
          + (c-a)*hat_x[1]
 
72
          + (d-a)*hat_x[2];
 
73
}
 
74
// ------------------------------------------
 
75
// piola transformation matrix: M_K where
 
76
//  F_K(hat_x) = M_K*hat_x + t_K
 
77
// ------------------------------------------
 
78
inline
 
79
void
 
80
set_piola_matrix_e (
 
81
  tensor& m,
 
82
  const point& a, 
 
83
  const point& b)
 
84
{
 
85
  m(0,0) = b[0]-a[0];
 
86
}
 
87
inline
 
88
void
 
89
set_piola_matrix_t (
 
90
  tensor& m,
 
91
  const point& a, 
 
92
  const point& b,
 
93
  const point& c)
 
94
{
 
95
   m.set_column (b-a, 0, 2);
 
96
   m.set_column (c-a, 1, 2);
 
97
}
 
98
inline
 
99
void
 
100
set_piola_matrix_T (
 
101
  tensor& m,
 
102
  const point& a, 
 
103
  const point& b,
 
104
  const point& c,
 
105
  const point& d)
 
106
{
 
107
   m.set_column (b-a, 0, 3);
 
108
   m.set_column (c-a, 1, 3);
 
109
   m.set_column (d-a, 2, 3);
 
110
}
 
111
inline
 
112
void
 
113
set_piola_matrix (
 
114
  tensor& M_K,
 
115
  const geo& omega,
 
116
  const geo_element& K)
 
117
{
 
118
    switch (K.type()) {
 
119
      case reference_element::e:
 
120
        set_piola_matrix_e (M_K,
 
121
           omega.vertex(K[0]),
 
122
           omega.vertex(K[1]));
 
123
        break;
 
124
      case reference_element::t:
 
125
        set_piola_matrix_t (M_K,
 
126
           omega.vertex(K[0]),
 
127
           omega.vertex(K[1]),
 
128
           omega.vertex(K[2]));
 
129
        break;
 
130
      case reference_element::T:
 
131
        set_piola_matrix_T (M_K,
 
132
           omega.vertex(K[0]),
 
133
           omega.vertex(K[1]),
 
134
           omega.vertex(K[2]),
 
135
           omega.vertex(K[3]));
 
136
        break;
 
137
      default:
 
138
        error_macro ("not yet supported element");
 
139
    }
 
140
}
 
141
// ------------------------------------------
 
142
// inverse piola transformation
 
143
// F_K^{-1} : K --> hat(K)
 
144
//            x --> hat(x)
 
145
// ------------------------------------------
 
146
inline
 
147
point 
 
148
inv_piola_e (
 
149
  const point& x, 
 
150
  const point& a, 
 
151
  const point& b) 
 
152
{
 
153
  return point((x[0]-a[0])/(b[0]-a[0]));
 
154
}
 
155
inline
 
156
point 
 
157
inv_piola_t (
 
158
  const point& x, 
 
159
  const point& a, 
 
160
  const point& b, 
 
161
  const point& c) 
 
162
{
 
163
  Float t9 = 1/(-b[0]*c[1]+b[0]*a[1]+a[0]*c[1]+c[0]*b[1]-c[0]*a[1]-a[0]*b[1]);
 
164
  Float t11 = -a[0]+x[0];
 
165
  Float t15 = -a[1]+x[1];
 
166
  return point((-c[1]+a[1])*t9*t11-(-c[0]+a[0])*t9*t15,
 
167
               (b[1]-a[1])*t9*t11-(b[0]-a[0])*t9*t15);
 
168
}
 
169
inline
 
170
point 
 
171
inv_piola_T (
 
172
  const point& x, 
 
173
  const point& a, 
 
174
  const point& b, 
 
175
  const point& c, 
 
176
  const point& d) 
 
177
{
 
178
  tensor A;
 
179
  point ax;
 
180
  for (size_t i = 0; i < 3; i++) {
 
181
    ax[i]  = x[i]-a[i];
 
182
    A(i,0) = b[i]-a[i];
 
183
    A(i,1) = c[i]-a[i];
 
184
    A(i,2) = d[i]-a[i];
 
185
  }
 
186
  tensor inv_A;
 
187
  bool is_singular = ! invert_3x3 (A, inv_A);
 
188
  check_macro(!is_singular, "inv_piola: singular transformation in tetrahedron");
 
189
  point hat_x = inv_A*ax;
 
190
  return hat_x; 
 
191
}
 
192
#endif // _RHEO_PIOLA_ALGO_H