← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rheolef/saucy-proposed

« back to all changes in this revision

Viewing changes to nfem/plib/field_element.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/plib/field_element.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
#include "rheolef/field_element.h"
 
22
#include "rheolef/field_evaluate.h"
 
23
namespace rheolef {
 
24
 
 
25
// cstor by specifying the quadrature formulae
 
26
template<class T, class M>
 
27
void
 
28
field_element<T,M>::initialize (const space_basic<T,M>& X, const quadrature_option_type& qopt) const {
 
29
  _X = X;
 
30
  if (qopt.get_order() != std::numeric_limits<quadrature_option_type::size_type>::max()) {
 
31
    _quad.set_order (qopt.get_order());
 
32
  } else {
 
33
    size_type k = X.get_numbering().degree();
 
34
    _quad.set_order (2*k-1); // see Raviart & Thomas, Masson, 1988, page 130, theorem 5.3.1
 
35
  }
 
36
  _quad.set_family (qopt.get_family());
 
37
 
 
38
trace_macro ("quadrature : " << _quad.get_family_name() << " order " << _quad.get_order());
 
39
 
 
40
  _basis_on_quad.set (_quad, get_basis());
 
41
  _piola_on_quad.set (_quad, get_piola_basis());
 
42
  _initialized = true;
 
43
}
 
44
// compute Mk(i,q) such that mk(phi_i,f) = sum_q Mk(i,q)*f(q)*wq
 
45
template<class T, class M>
 
46
void
 
47
field_element<T,M>::eval (
 
48
    const geo_element&    K,
 
49
    const std::vector<T>& fq,
 
50
    std::vector<T>&       uk) const
 
51
{
 
52
  std::vector<size_type> dis_inod;
 
53
  get_space().get_geo().dis_inod (K, dis_inod);
 
54
  size_type d = get_space().get_geo().dimension();
 
55
  size_type map_d = K.dimension();
 
56
  size_t b_size = get_basis().size(K);
 
57
  uk.resize (b_size);
 
58
  fill (uk.begin(), uk.end(), T(0.0));
 
59
  quadrature::const_iterator first_quad = _quad.begin(K);
 
60
  quadrature::const_iterator last_quad  = _quad.end  (K);
 
61
  tensor DF;
 
62
  for (size_type q = 0; first_quad != last_quad; first_quad++, q++) {
 
63
    // TODO: passer (hat_x, K) a` field_o_characteristic: ca evite de localiser
 
64
    // si on a vh sur le meme mailage
 
65
    jacobian_piola_transformation (get_space().get_geo(), _piola_on_quad, K, dis_inod, q, DF);
 
66
    T det_DF = det_jacobian_piola_transformation (DF, d, map_d);
 
67
    T wq = 1;
 
68
#ifdef TODO
 
69
    wq *= weight_coordinate_system (K, q); // TODO axi
 
70
#endif // TODO
 
71
    wq *= det_DF;
 
72
    wq *= (*first_quad).w;
 
73
    Float fwq = fq[q]*wq;
 
74
    basis_on_quadrature::const_iterator first_phi = _basis_on_quad.begin (K, q);
 
75
    basis_on_quadrature::const_iterator last_phi  = _basis_on_quad.end   (K, q);
 
76
    for (size_t i = 0; first_phi != last_phi; first_phi++, i++) {
 
77
      Float phi_wq = (*first_phi)*wq;
 
78
      uk[i] += (*first_phi)*fwq;
 
79
    }
 
80
  }
 
81
}
 
82
// idem for vector-valued functions
 
83
template<class T, class M>
 
84
void
 
85
field_element<T,M>::eval (
 
86
    const geo_element&                  K,
 
87
    const std::vector<point_basic<T> >& fq,
 
88
    std::vector<T>&                     uk) const
 
89
{
 
90
  std::vector<size_type> dis_inod;
 
91
  get_space().get_geo().dis_inod (K, dis_inod);
 
92
  size_type map_d = K.dimension();
 
93
  size_type d = _X.get_geo().dimension();
 
94
  size_t b_size = get_basis().size(K);
 
95
  size_type n_comp = d;
 
96
  uk.resize (n_comp*b_size);
 
97
  fill (uk.begin(), uk.end(), Float(0.0));
 
98
  quadrature::const_iterator first_quad = _quad.begin(K);
 
99
  quadrature::const_iterator last_quad  = _quad.end  (K);
 
100
  tensor DF;
 
101
  for (size_type q = 0; first_quad != last_quad; first_quad++, q++) {
 
102
    // TODO: passer (hat_x, K) a` field_o_characteristic: ca evite de localiser
 
103
    // si on a vh sur le meme mailage
 
104
    jacobian_piola_transformation (get_space().get_geo(), _piola_on_quad, K, dis_inod, q, DF);
 
105
    T det_DF = det_jacobian_piola_transformation (DF, d, map_d);
 
106
    Float wq = 1;
 
107
#ifdef TODO
 
108
    wq *= weight_coordinate_system (K, q); // TODO axi
 
109
#endif // TODO
 
110
    wq *= det_DF;
 
111
    wq *= (*first_quad).w; 
 
112
    point_basic<T> fwq = fq[q]*wq;
 
113
    basis_on_quadrature::const_iterator first_phi = _basis_on_quad.begin (K, q);
 
114
    basis_on_quadrature::const_iterator last_phi  = _basis_on_quad.end (K, q);
 
115
    for (size_t i = 0; first_phi != last_phi; first_phi++, i++) {
 
116
      Float phi_wq = (*first_phi)*wq;
 
117
      for (size_t i_comp = 0; i_comp < n_comp; i_comp++) {
 
118
        uk[i_comp*b_size + i] += (*first_phi)*fwq[i_comp];
 
119
      }
 
120
    }
 
121
  }
 
122
}
 
123
#ifdef TODO
 
124
template<>
 
125
void 
 
126
field_element::operator() (const geo_element& K,
 
127
        const field_o_characteristic& f,
 
128
        std::vector<Float>& uk) const
 
129
{
 
130
  switch (f.get_field().get_valued_type()) {
 
131
    case fem_helper::scalar: {
 
132
      std::vector<Float> fq(_quad.size(K));
 
133
      quadrature::const_iterator first = _quad.begin(K);
 
134
      quadrature::const_iterator last  = _quad.end  (K);
 
135
      for (size_type q = 0; first != last; first++, q++) {
 
136
        point xq = get_space().get_geo().dehatter ((*first).x, K.index());
 
137
        fq[q] = f(xq);
 
138
      }
 
139
      eval (K, fq, uk);
 
140
      break;
 
141
    }
 
142
    case fem_helper::vectorial: {
 
143
      std::vector<point> fq(_quad.size(K));
 
144
      quadrature::const_iterator first = _quad.begin(K);
 
145
      quadrature::const_iterator last  = _quad.end  (K);
 
146
      for (size_type q = 0; first != last; first++, q++) {
 
147
        point xq = get_space().get_geo().dehatter ((*first).x, K.index());
 
148
        fq[q] = f.vector_evaluate(xq);
 
149
      }
 
150
      eval (K, fq, uk);
 
151
      break;
 
152
    }
 
153
    default:
 
154
      error_macro ("unsupported");
 
155
  }
 
156
}
 
157
#endif // TODO
 
158
template<class T, class M>
 
159
void 
 
160
field_element<T,M>::operator() (
 
161
        const geo_element&      K,
 
162
        const field_basic<T,M>& fh,
 
163
        std::vector<T>&         mfk) const
 
164
{
 
165
#ifdef TODO
 
166
  switch (fh.get_valued_type()) {
 
167
    case fem_helper::scalar: {
 
168
#endif // TODO
 
169
      // 1) evaluate f at all integration points xq ; fem basis is evaluated at xq one time for all
 
170
      std::vector<T> fq(_quad.size(K));
 
171
      std::vector<size_type> dis_idof;
 
172
      const numbering<T,M>& fem = get_space().get_numbering();
 
173
      fem.dis_idof (get_space().get_geo().sizes(), K, dis_idof);
 
174
      for (size_type q = 0, nq = _quad.size(K); q < nq; q++) {
 
175
        fq[q] = field_evaluate (fh, _basis_on_quad, K, dis_idof, q);
 
176
      }
 
177
      // 2) integrate f using quadrature formulae (xq,wq)
 
178
      eval (K, fq, mfk);
 
179
#ifdef TODO
 
180
      break;
 
181
    }
 
182
    case fem_helper::vectorial: {
 
183
      size_t d = fh.dimension();
 
184
      std::vector<point> fq(_quad.size(K));
 
185
      quadrature::const_iterator first = _quad.begin(K);
 
186
      quadrature::const_iterator last  = _quad.end  (K);
 
187
      for (size_type q = 0; first != last; first++, q++) {
 
188
        meshpoint xq (K.index(), (*first).x);
 
189
        for (size_t i = 0; i < d; i++) {
 
190
          fq[q][i] = fh.evaluate (xq, i);
 
191
        }
 
192
      }
 
193
      eval (K, fq, mfk);
 
194
      break;
 
195
    }
 
196
    default:
 
197
      error_macro ("unsupported yet");
 
198
  }
 
199
#endif // TODO
 
200
}
 
201
// ==========================================================================
 
202
// instanciation in library
 
203
// ==========================================================================
 
204
template class field_element<Float,sequential>;
 
205
 
 
206
#ifdef _RHEOLEF_HAVE_MPI
 
207
template class field_element<Float,distributed>;
 
208
#endif // _RHEOLEF_HAVE_MPI
 
209
 
 
210
}// namespace rheolef