2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
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.
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.
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
20
/// =========================================================================
21
#include "rheolef/field_element.h"
22
#include "rheolef/field_evaluate.h"
25
// cstor by specifying the quadrature formulae
26
template<class T, class M>
28
field_element<T,M>::initialize (const space_basic<T,M>& X, const quadrature_option_type& qopt) const {
30
if (qopt.get_order() != std::numeric_limits<quadrature_option_type::size_type>::max()) {
31
_quad.set_order (qopt.get_order());
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
36
_quad.set_family (qopt.get_family());
38
trace_macro ("quadrature : " << _quad.get_family_name() << " order " << _quad.get_order());
40
_basis_on_quad.set (_quad, get_basis());
41
_piola_on_quad.set (_quad, get_piola_basis());
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>
47
field_element<T,M>::eval (
49
const std::vector<T>& fq,
50
std::vector<T>& uk) const
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);
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);
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);
69
wq *= weight_coordinate_system (K, q); // TODO axi
72
wq *= (*first_quad).w;
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;
82
// idem for vector-valued functions
83
template<class T, class M>
85
field_element<T,M>::eval (
87
const std::vector<point_basic<T> >& fq,
88
std::vector<T>& uk) const
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);
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);
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);
108
wq *= weight_coordinate_system (K, q); // TODO axi
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];
126
field_element::operator() (const geo_element& K,
127
const field_o_characteristic& f,
128
std::vector<Float>& uk) const
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());
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);
154
error_macro ("unsupported");
158
template<class T, class M>
160
field_element<T,M>::operator() (
161
const geo_element& K,
162
const field_basic<T,M>& fh,
163
std::vector<T>& mfk) const
166
switch (fh.get_valued_type()) {
167
case fem_helper::scalar: {
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);
177
// 2) integrate f using quadrature formulae (xq,wq)
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);
197
error_macro ("unsupported yet");
201
// ==========================================================================
202
// instanciation in library
203
// ==========================================================================
204
template class field_element<Float,sequential>;
206
#ifdef _RHEOLEF_HAVE_MPI
207
template class field_element<Float,distributed>;
208
#endif // _RHEOLEF_HAVE_MPI
210
}// namespace rheolef