← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/piola.h

  • 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/plib/piola.h
 
1
#ifndef _RHEOLEF_PIOLA_H
 
2
#define _RHEOLEF_PIOLA_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
// piola on a predefined point set: hat_x[q], q=0..nq
 
25
//
 
26
#include "rheolef/geo.h"
 
27
#include "rheolef/basis_on_lattice.h"
 
28
namespace rheolef { 
 
29
 
 
30
template<class T, class M, class BasisOnPointset>
 
31
point_basic<T>
 
32
piola_transformation (
 
33
  const geo_basic<T,M>&         omega,
 
34
  const BasisOnPointset&        piola_on_pointset,
 
35
  reference_element             hat_K,
 
36
  const std::vector<size_t>&    dis_inod,
 
37
  size_t                        q)
 
38
{
 
39
  typedef typename geo_basic<T,M>::size_type size_type;
 
40
  size_type d = omega.dimension();
 
41
  size_type loc_inod = 0;
 
42
  point_basic<T> xq;
 
43
  for (typename BasisOnPointset::const_iterator
 
44
          iter = piola_on_pointset.begin(hat_K, q),
 
45
          last = piola_on_pointset.end  (hat_K, q);
 
46
          iter != last;
 
47
          iter++, loc_inod++) {
 
48
    const T& cnod = *iter;
 
49
    const point_basic<T>& xnod = omega.dis_node (dis_inod[loc_inod]);
 
50
    for (size_type i = 0; i < d; i++) {
 
51
      xq[i] += cnod*xnod[i];
 
52
    }
 
53
  }
 
54
  return xq;
 
55
}
 
56
template<class T, class M, class BasisOnPointset>
 
57
void
 
58
jacobian_piola_transformation (
 
59
  const geo_basic<T,M>&         omega,
 
60
  const BasisOnPointset&        piola_on_pointset,
 
61
  reference_element             hat_K,
 
62
  const std::vector<size_t>&    dis_inod,
 
63
  size_t                        q,
 
64
  tensor&                       DF)
 
65
{
 
66
  typedef typename geo_basic<T,M>::size_type size_type;
 
67
  size_type d = omega.dimension();
 
68
  size_type map_d = hat_K.dimension();
 
69
  size_type loc_inod = 0;
 
70
  DF.fill (0);
 
71
  for (typename BasisOnPointset::const_iterator_grad
 
72
        first_grad_tr = piola_on_pointset.begin_grad (hat_K, q),
 
73
         last_grad_tr = piola_on_pointset.end_grad   (hat_K, q);
 
74
        first_grad_tr != last_grad_tr;
 
75
        first_grad_tr++, loc_inod++) {
 
76
    const point_basic<T>& xnod = omega.dis_node (dis_inod[loc_inod]);
 
77
    cumul_otimes (DF, xnod, *first_grad_tr, d, map_d);
 
78
  }
 
79
}
 
80
// TODO: template<class T> with tensor and result
 
81
inline
 
82
Float
 
83
det_jacobian_piola_transformation (const tensor& DF, size_t d , size_t map_d)
 
84
{
 
85
    if (d == map_d) {
 
86
      return DF.determinant (map_d);
 
87
    }
 
88
    /* surface jacobian: references:
 
89
     * Spectral/hp element methods for CFD
 
90
     * G. E. M. Karniadakis and S. J. Sherwin
 
91
     * Oxford university press
 
92
     * 1999
 
93
     * page 165
 
94
     */
 
95
    switch (map_d) {
 
96
      case 0: return 1;
 
97
      case 1: return norm(DF.col(0));
 
98
      case 2: return norm(vect(DF.col(0), DF.col(1)));
 
99
      default:
 
100
        error_macro ("det_jacobian_piola_transformation: unsupported element dimension "
 
101
            << map_d << " in " << d << "D mesh.");
 
102
        return 0;
 
103
    }
 
104
}
 
105
// The pseudo inverse extend inv(DF) for face in 3d or edge in 2d
 
106
//  i.e. usefull for Laplacian-Beltrami and others surfacic forms.
 
107
//
 
108
// pinvDF (hat_xq) = inv DF, if tetra in 3d, tri in 2d, etc
 
109
//                  = pseudo-invese, when tri in 3d, edge in 2 or 3d
 
110
// e.g. on 3d face : pinvDF*DF = [1, 0, 0; 0, 1, 0; 0, 0, 0]
 
111
//
 
112
// let DF = [u, v, w], where u, v, w are the column vectors of DF
 
113
// then det(DF) = mixt(u,v,w)
 
114
// and det(DF)*inv(DF)^T = [v^w, w^u, u^v] where u^v = vect(u,v)
 
115
//
 
116
// application:
 
117
// if K=triangle(a,b,c) then u=ab=b-a, v=ac=c-a and w = n = u^v/|u^v|.
 
118
// Thus DF = [ab,ac,n] and det(DF)=|ab^ac|
 
119
// and inv(DF)^T = [ac^n/|ab^ac|, -ab^n/|ab^ac|, n]
 
120
// The pseudo-inverse is obtained by remplacing the last column n by zero.
 
121
//
 
122
// TODO: template<class T> with tensor
 
123
inline
 
124
tensor
 
125
pseudo_inverse_jacobian_piola_transformation (
 
126
    const tensor& DF,
 
127
    size_t d,
 
128
    size_t map_d)
 
129
{
 
130
    if (d == map_d) {
 
131
      return inv(DF, map_d);
 
132
    }
 
133
    tensor invDF;
 
134
    switch (map_d) {
 
135
      case 0: { // point in 1D
 
136
          invDF(0,0) = 1;
 
137
          return invDF;
 
138
      }
 
139
      case 1: { // segment in 2D
 
140
          point t = DF.col(0);
 
141
          invDF.set_row (t/norm2(t), 0, d);
 
142
          return invDF;
 
143
      }
 
144
      case 2: {
 
145
          point t0 = DF.col(0);
 
146
          point t1 = DF.col(1);
 
147
          point n = vect(t0,t1);
 
148
          Float det2 = norm2(n);
 
149
          point v0 =   vect(t1,n)/det2;
 
150
          point v1 = - vect(t0,n)/det2;
 
151
          invDF.set_row (v0, 0, d);
 
152
          invDF.set_row (v1, 1, d);
 
153
          return invDF;
 
154
      }
 
155
      default:
 
156
          error_macro ("pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
 
157
            << map_d << " in " << d << "D mesh.");
 
158
          return invDF;
 
159
    }
 
160
}
 
161
// x --> hat_x  on K
 
162
template<class T, class M>
 
163
point_basic<T>
 
164
inverse_piola_transformation (
 
165
  const geo_basic<T,M>&         omega,
 
166
  reference_element             hat_K,
 
167
  const std::vector<size_t>&    dis_inod,
 
168
  const point_basic<T>&         x);
 
169
 
 
170
// axisymetric weight ?
 
171
// point_basic<T> xq = rheolef::piola_transformation (_omega, _piola_table, K, dis_inod, q);
 
172
template<class T>
 
173
inline
 
174
T
 
175
weight_coordinate_system (space_constant::coordinate_type sys_coord, const point_basic<T>& xq)
 
176
{
 
177
  switch (sys_coord) {
 
178
    case space_constant::axisymmetric_rz: return xq[0];
 
179
    case space_constant::axisymmetric_zr: return xq[1];
 
180
    case space_constant::cartesian:       return 1;
 
181
    default: {
 
182
      fatal_macro ("unsupported coordinate system `"
 
183
        <<  space_constant::coordinate_system_name (sys_coord) << "'");
 
184
      return 0;
 
185
    }
 
186
  }
 
187
}
 
188
 
 
189
}// namespace rheolef
 
190
#endif // _RHEOLEF_PIOLA_H