← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/form_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/form_element.cc
20
20
/// =========================================================================
21
21
#include "rheolef/form_element.h"
22
22
#include "rheolef/tensor.h"
 
23
#include "rheolef/geo_domain.h"
 
24
#include "rheolef/piola.h"
 
25
#include "rheolef/ublas_matrix_range.h"
23
26
 
24
27
namespace rheolef {
25
28
using namespace std;
29
32
// ==========================================================================
30
33
template<class T, class M>
31
34
void
32
 
form_element_rep<T,M>::initialize_all ()
 
35
form_element_rep<T,M>::initialize_all () const
33
36
{
34
37
  initialize();
35
38
#ifdef TODO
36
39
  check_macro (_X.get_global_geo() == _Y.get_global_geo(),
37
40
        "form assembly: global spaces based on different meshes: unsupported");
38
41
#endif // TODO
39
 
 
40
 
  size_type quad_order = get_first_basis().degree() + get_second_basis().degree();
41
 
#ifdef TODO
42
 
  if (coordinate_system_type() != fem_helper::cartesian) quad_order++;
43
 
  if (_is_weighted) quad_order += get_weight_basis().degree();
44
 
#endif // TODO
45
 
  quad_order -= n_derivative();
46
 
  _quad.set_order (quad_order);
 
42
  if (_qopt.get_order() == 0 ||
 
43
      _qopt.get_order() == std::numeric_limits<quadrature_option_type::size_type>::max()) {
 
44
    size_t k1 = get_first_basis().degree();
 
45
    size_t k2 = get_second_basis().degree();
 
46
    size_type quad_order = k1 + k2 - n_derivative() + 1;
 
47
    if (is_weighted()) {
 
48
      size_type k3 = get_weight().get_space().get_numbering().get_basis().degree();
 
49
      quad_order += k3;
 
50
    } else {
 
51
    }
 
52
    if (coordinate_system() != space_constant::cartesian) quad_order++;
 
53
    if (_omega.get_background_geo().sizes().ownership_by_variant[reference_element::q].dis_size() != 0 ||
 
54
        _omega.get_background_geo().sizes().ownership_by_variant[reference_element::P].dis_size() != 0 ||
 
55
        _omega.get_background_geo().sizes().ownership_by_variant[reference_element::H].dis_size() != 0) {
 
56
          // integrate exactly ??
 
57
          quad_order += 2;
 
58
    }
 
59
    _qopt.set_order (quad_order);
 
60
  } else {
 
61
  }
 
62
  if (_qopt.get_family() == quadrature_option_type::max_family) {
 
63
    _qopt.set_family (quadrature_option_type::gauss);
 
64
  }
 
65
  _quad.set_order  (_qopt.get_order());
 
66
  _quad.set_family (_qopt.get_family());
 
67
trace_macro ("quadrature : " << _quad.get_family_name() << " order " << _quad.get_order());
47
68
  _bx_table.set (_quad, get_first_basis());
48
69
  _by_table.set (_quad, get_second_basis());
49
 
  _tr_p1_table.set (_quad, get_p1_transformation());
 
70
  _piola_table.set (_quad, get_piola_basis());
50
71
#ifdef TODO
51
72
  if (_is_weighted) {
52
73
    check_macro (get_weight().get_space().get_global_geo() == get_global_geo(),
57
78
  _initialized = true;
58
79
}
59
80
// ==========================================================================
60
 
// part 2) piola transform
 
81
// part 2) axisymmetric coord system handler
61
82
// ==========================================================================
62
83
template<class T, class M>
63
84
T
64
 
form_element_rep<T,M>::det_jacobian_piola_transformation (const tensor& DF, size_type map_dim) const
 
85
form_element_rep<T,M>::weight_coordinate_system (
 
86
    const geo_element&            K,
 
87
    const std::vector<size_type>& dis_inod,
 
88
    size_type                     q) const
65
89
{
66
 
    if (coordinate_dimension() == map_dim) {
67
 
      return DF.determinant (map_dim);
68
 
    }
69
 
    /* surface jacobian: references:
70
 
     * Spectral/hp element methods for CFD
71
 
     * G. E. M. Karniadakis and S. J. Sherwin
72
 
     * Oxford university press
73
 
     * 1999
74
 
     * page 165
75
 
     */
76
 
    switch (map_dim) {
77
 
      case 0: return 1;
78
 
      case 1: return norm(DF.col(0));
79
 
      case 2: return norm(vect(DF.col(0), DF.col(1)));
80
 
      default:
81
 
        error_macro ("det_jacobian_piola_transformation: unsupported element dimension "
82
 
            << map_dim << " in " << coordinate_dimension() << "D mesh.");
 
90
  if (coordinate_system() == space_constant::cartesian || ! use_coordinate_system_weight()) {
 
91
    return 1;
 
92
  }
 
93
  point_basic<T> xq = rheolef::piola_transformation (_omega, _piola_table, K, dis_inod, q);
 
94
  switch (coordinate_system()) {
 
95
    case space_constant::axisymmetric_rz: 
 
96
          return use_coordinate_system_dual_weight() ? 1/xq[0] : xq[0];
 
97
    case space_constant::axisymmetric_zr:
 
98
          return use_coordinate_system_dual_weight() ? 1/xq[1] : xq[1];
 
99
    default: {
 
100
        fatal_macro ("unsupported coordinate system `" 
 
101
        <<  space_constant::coordinate_system_name (coordinate_system()) << "'");
83
102
        return 0;
84
103
    }
85
 
}
86
 
// ----------------------------------------------------------
87
 
// DF (hat_xq)
88
 
// ----------------------------------------------------------
89
 
template<class T, class M>
90
 
void
91
 
form_element_rep<T,M>::jacobian_piola_transformation (const geo_element& K, size_type q, tensor& DF) const
92
 
{
93
 
    size_type nt = get_p1_transformation().size(K);
94
 
    check_macro (nt == K.size(), "unsupported transformation " << get_p1_transformation().name());
95
 
    DF.fill (0);
96
 
    basis_on_quadrature::const_iterator_grad grad_tr
97
 
     = _tr_p1_table.begin_grad (K, q);
98
 
    size_type coord_d = coordinate_dimension();
99
 
    for (size_type k = 0; k < nt; k++, grad_tr++)
100
 
      my_cumul_otimes (DF, dis_vertex(K[k]), *grad_tr, coord_d, K.dimension());
101
 
}
102
 
// The pseudo inverse extend inv(DF) for face in 3d or edge in 2d
103
 
//  i.e. usefull for Laplacian-Beltrami and others surfacic forms.
104
 
//
105
 
// pinvDF (hat_xq) = inv DF, if tetra in 3d, tri in 2d, etc
106
 
//                  = pseudo-invese, when tri in 3d, edge in 2 or 3d
107
 
// e.g. on 3d face : pinvDF*DF = [1, 0, 0; 0, 1, 0; 0, 0, 0]
108
 
//
109
 
// let DF = [u, v, w], where u, v, w are the column vectors of DF
110
 
// then det(DF) = mixt(u,v,w)
111
 
// and det(DF)*inv(DF)^T = [v^w, w^u, u^v] where u^v = vect(u,v)
112
 
//
113
 
// application:
114
 
// if K=triangle(a,b,c) then u=ab=b-a, v=ac=c-a and w = n = u^v/|u^v|.
115
 
// Thus DF = [ab,ac,n] and det(DF)=|ab^ac|
116
 
// and inv(DF)^T = [ac^n/|ab^ac|, -ab^n/|ab^ac|, n]
117
 
// The pseudo-inverse is obtained by remplacing the last column n by zero.
118
 
//
119
 
template<class T, class M>
120
 
tensor
121
 
form_element_rep<T,M>::pseudo_inverse_jacobian_piola_transformation (
122
 
    const tensor& DF,
123
 
    size_type map_dim) const
124
 
{
125
 
    size_t d = coordinate_dimension();
126
 
    if (d == map_dim) {
127
 
      return inv(DF, map_dim);
128
 
    }
129
 
    tensor invDF;
130
 
    switch (map_dim) {
131
 
      case 0: { // point in 1D
132
 
          invDF(0,0) = 1;
133
 
          return invDF;
134
 
      }
135
 
      case 1: { // segment in 2D
136
 
          point t = DF.col(0);
137
 
          invDF.set_row (t/norm2(t), 0, d);
138
 
          return invDF;
139
 
      }
140
 
      case 2: {
141
 
          point t0 = DF.col(0);
142
 
          point t1 = DF.col(1);
143
 
          point n = vect(t0,t1);
144
 
          T det2 = norm2(n);
145
 
          point v0 =   vect(t1,n)/det2;
146
 
          point v1 = - vect(t0,n)/det2;
147
 
          invDF.set_row (v0, 0, d);
148
 
          invDF.set_row (v1, 1, d);
149
 
          return invDF;
150
 
      }
151
 
      default:
152
 
          error_macro ("pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
153
 
            << map_dim << " in " << coordinate_dimension() << "D mesh.");
154
 
          return invDF;
155
 
    }
 
104
  }
156
105
}
157
106
// ==========================================================================
158
107
// part 3) some few algebra on elementary matrices
170
119
  tensor::size_type map_d)
171
120
{
172
121
  typedef tensor::size_type size_type;
173
 
  Float sum = 0;
 
122
  typedef Float T; // TODO tensor_basic<T>
 
123
  T sum = 0;
174
124
  for (size_type k = 0; k < map_d; k++) {
175
125
    for (size_type l = 0; l < map_d; l++) {
176
126
      sum += v[k]*t(k,l)*u[l];
233
183
void
234
184
form_element_rep<T,M>::build_scalar_mass (const geo_element& K, ublas::matrix<T>& m) const
235
185
{
236
 
  reference_element hat_K (K.type());
 
186
  std::vector<size_type> dis_inod;
 
187
  _omega.dis_inod (K, dis_inod);
 
188
  reference_element hat_K (K.variant());
237
189
  m.resize (get_second_basis().size(hat_K), get_first_basis().size(hat_K));
238
190
  m.clear();
239
191
  quadrature::const_iterator first = _quad.begin(hat_K);
241
193
  tensor DF;
242
194
  for (size_type q = 0; first != last; first++, q++) {
243
195
    T wq = 1;
244
 
#ifdef TODO
245
 
    wq *= is_weighted() ? weight(K,q,weight_i_comp) : T(1);
246
 
    wq *= weight_coordinate_system (K, q);
247
 
#endif // TODO
248
 
    jacobian_piola_transformation (K, q, DF);
249
 
    wq *= det_jacobian_piola_transformation (DF, K.dimension());
 
196
    wq *= is_weighted() ? weight(K,q) : T(1);
 
197
    wq *= weight_coordinate_system (K, dis_inod, q);
 
198
    rheolef::jacobian_piola_transformation (_omega, _piola_table, K, dis_inod, q, DF);
 
199
    wq *= rheolef::det_jacobian_piola_transformation (DF, coordinate_dimension(), K.dimension());
250
200
    wq *= (*first).w; 
 
201
#ifdef TO_CLEAN
 
202
    const point& hat_xq = (*first).x; 
 
203
    cerr << setprecision(15) 
 
204
         << "DF(hat_x"<<q<<"="; hat_xq.put (cerr,2); cerr << ")="; DF.put(cerr,2); cerr << endl;
 
205
#endif // TO_CLEAN
251
206
    cumul_otimes (m, wq, 
252
207
        _by_table.begin (hat_K, q), _by_table.end (hat_K, q),
253
208
        _bx_table.begin (hat_K, q), _bx_table.end (hat_K, q));
254
209
  }
255
210
}
 
211
// --------------------------------------------------------------------
 
212
// general local mass matrix : multi-component spaces and weight
 
213
// --------------------------------------------------------------------
 
214
template<class T, class M>
 
215
void
 
216
form_element_rep<T,M>::build_general_mass (const geo_element& K, ublas::matrix<T>& m) const
 
217
{
 
218
  space_constant::valued_type space_valued  = get_first_space().valued_tag();
 
219
  if (space_valued == space_constant::scalar) {
 
220
    build_scalar_mass (K, m);
 
221
    return;
 
222
  }
 
223
  // space is vectorial or tensorial: weight could be scalar, tensor or tensor4
 
224
  size_type n_comp = get_first_space().size();
 
225
  size_type loc_n1dof = get_second_basis().size(K);
 
226
  size_type loc_n2dof = get_first_basis().size(K);
 
227
  m.resize (n_comp*loc_n1dof, n_comp*loc_n2dof);
 
228
  m.clear();
 
229
  space_constant::valued_type weight_valued = space_constant::scalar;
 
230
  if (is_weighted()) {
 
231
    weight_valued = _wh.get_space().get_constitution().valued_tag();
 
232
  }
 
233
  if (weight_valued == space_constant::scalar) {
 
234
    // --------------------------------------------------------------------------
 
235
    // all components have the same weight and the local matrix is block-diagonal
 
236
    // --------------------------------------------------------------------------
 
237
    ublas::matrix<T> m_ij;
 
238
    build_scalar_mass (K, m_ij);
 
239
    switch (space_valued) {
 
240
      case space_constant::scalar:
 
241
      case space_constant::vector: {
 
242
        for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
 
243
          mr_set (m, ublas::range(i_comp*loc_n1dof,(i_comp+1)*loc_n1dof), ublas::range(i_comp*loc_n2dof,(i_comp+1)*loc_n2dof), m_ij);
 
244
        }
 
245
        break;
 
246
      }
 
247
      case space_constant::tensor: {
 
248
        // symmetric tensor => 2 factor for non-diagonal subscripts
 
249
        for (size_type ij_comp = 0; ij_comp < n_comp; ij_comp++) {
 
250
          space_constant::coordinate_type sys_coord = coordinate_system();
 
251
          std::pair<size_type,size_type> ij = space_constant::tensor_subscript (space_valued, sys_coord, ij_comp);
 
252
          size_t i_comp = ij.first;
 
253
          size_t j_comp = ij.second;
 
254
          T factor_ij = ((i_comp == j_comp) ? 1 : 2); // symmetry => multiplicity factor: when i_comp != j_comp : 2*w_ij
 
255
          mr_set (m, ublas::range(ij_comp*loc_n1dof,(ij_comp+1)*loc_n1dof), ublas::range(ij_comp*loc_n2dof,(ij_comp+1)*loc_n2dof), factor_ij*m_ij);
 
256
        }
 
257
        break;
 
258
      }
 
259
      default: {
 
260
        error_macro ("unexpected `" << space_constant::valued_name(space_valued) << "'-valued space for the `mass' form.");
 
261
      }
 
262
    }
 
263
  } else { // weight_valued != fem_helper::scalar
 
264
    // ------------------------------------------------------------------------------
 
265
    // components have different weight : the local matrix is no more block-diagonal
 
266
    // ------------------------------------------------------------------------------
 
267
#ifdef TODO
 
268
    switch (space_valued) {
 
269
      case fem_helper::vectorial: {
 
270
        // first case : vector space and tensor weight
 
271
        //     m(u,v) = int_Omega w_ij u_j v_i dx
 
272
        //     and w_ij = w_ji
 
273
        ublas::matrix<T> m_ij;
 
274
        for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
 
275
          for (size_type j_comp = 0; j_comp < n_comp; j_comp++) {
 
276
            size_type ij_comp = fem_helper::tensor_index (weight_valued, sys_coord, i_comp, j_comp);
 
277
            build_scalar_mass (K, m_ij, ij_comp);
 
278
            mr_set (m, ublas::range(i_comp*loc_n1dof,(i_comp+1)*loc_n1dof), ublas::range(j_comp*loc_n2dof,(j_comp+1)*loc_n2dof), m_ij);
 
279
          }
 
280
        }
 
281
        break;
 
282
      }
 
283
      case fem_helper::tensorial: {
 
284
        // second case : tensor space and tensor4 weight
 
285
        //     m(tau,gamma) = int_Omega w_ijkl tau_kj gamma_ij dx
 
286
        //     and w_ijkl = w_jikl, w_ijkl = w_ijlk and w_ijkl = w_klij
 
287
        // symmetry => multiplicity factor: 4 when i!=j and k!=l or 2 when i!=j xor k!=l
 
288
        ublas::matrix<T> m_ijkl;
 
289
        for (size_type ij_comp = 0; ij_comp < n_comp; ij_comp++) {
 
290
          fem_helper::pair_size_type ij = fem_helper::tensor_subscript (space_valued, sys_coord, ij_comp);
 
291
          size_t i = ij.first;
 
292
          size_t j = ij.second;
 
293
          T factor_ij = ((i == j) ? 1 : 2);
 
294
          for (size_type kl_comp = 0; kl_comp < n_comp; kl_comp++) {
 
295
            fem_helper::pair_size_type kl = fem_helper::tensor_subscript (space_valued, sys_coord, kl_comp);
 
296
            size_t k = kl.first;
 
297
            size_t l = kl.second;
 
298
            T factor_kl = ((k == l) ? 1 : 2);
 
299
            T factor_ijkl = factor_ij*factor_kl;
 
300
            size_type ijkl_comp = fem_helper::tensor4_index (weight_valued, sys_coord, i, j, k, l);
 
301
            build_scalar_mass (K, m_ijkl, ijkl_comp);
 
302
            mr_set (m, ublas::range(ij_comp*loc_n1dof,(ij_comp+1)*loc_n1dof), ublas::range(kl_comp*loc_n2dof,(kl_comp+1)*loc_n2dof), factor_ijkl*m_ijkl);
 
303
          }
 
304
        }
 
305
        break;
 
306
      }
 
307
      default: {
 
308
        error_macro ("unexpected " << fem_helper::valued_field_name(space_valued) << "-valued space and `"
 
309
                << fem_helper::valued_field_name(weight_valued)
 
310
                << "'-valued weight for the `mass' form.");
 
311
      }
 
312
    }
 
313
#endif// TODO
 
314
  error_macro ("mass: multi-component with non-scalar weight: not yet");
 
315
  }
 
316
}
256
317
// ----------------------------------------------------------
257
318
// elementary grad_grad form:
258
319
// --------------------------------------------------------------------
260
321
void 
261
322
form_element_rep<T,M>::build_scalar_grad_grad (const geo_element& K, ublas::matrix<T>& a) const
262
323
{
263
 
  reference_element hat_K (K.type());
 
324
  bool scalar_diffusion = true;
 
325
  std::vector<size_type> dis_inod;
 
326
  _omega.dis_inod (K, dis_inod);
 
327
  reference_element hat_K (K.variant());
264
328
  size_type map_d = K.dimension();
265
329
  a.resize (get_second_basis().size(hat_K), get_first_basis().size(hat_K));
266
330
  a.clear();
271
335
  tensor Dw;
272
336
  tensor W;
273
337
  for (size_type q = 0; first != last; first++, q++) {
274
 
    jacobian_piola_transformation (K, q, DF);
275
 
    invDF = pseudo_inverse_jacobian_piola_transformation (DF, map_d);
 
338
    rheolef::jacobian_piola_transformation (_omega, _piola_table, K, dis_inod, q, DF);
 
339
    invDF = pseudo_inverse_jacobian_piola_transformation (DF, coordinate_dimension(), map_d);
276
340
    T wq = 1;
277
 
#ifdef TODO
278
 
    wq *= weight_coordinate_system (K, q);
279
 
#endif // TODO
280
 
    wq *= det_jacobian_piola_transformation (DF, K.dimension());
 
341
    wq *= weight_coordinate_system (K, dis_inod, q);
 
342
    wq *= rheolef::det_jacobian_piola_transformation (DF, coordinate_dimension(), K.dimension());
281
343
    wq *= (*first).w;
282
 
#ifdef TODO
283
344
    if (scalar_diffusion) {
284
345
       wq *= (is_weighted() ? weight(K,q) : T(1));
285
 
#endif // TODO
286
346
       Dw = wq*invDF*trans(invDF);
287
 
#ifdef TODO
288
347
    } else { // tensorial diffusion
 
348
       error_macro ("tensorial diffusion: not yet");
 
349
#ifdef TODO
289
350
       weight (K, q, W);
290
351
       Dw = wq*invDF*W*trans(invDF);
 
352
#endif // TODO
291
353
    }
292
 
#endif // TODO
293
354
    cumul_otimes (a, Dw,
294
355
      _by_table.begin_grad (hat_K, q), _by_table.end_grad (hat_K, q),
295
356
      _bx_table.begin_grad (hat_K, q), _bx_table.end_grad (hat_K, q),
296
357
      map_d);
297
358
  }
298
359
}
 
360
// --------------------------------------------------------------------
 
361
// (d_dxi phi, psi)
 
362
// Mij += trans(grad_phi_i) * Dw * grad_psi_j
 
363
// --------------------------------------------------------------------
 
364
template<class T, class M>
 
365
void 
 
366
form_element_rep<T,M>::build_d_dx (const geo_element& K, ublas::matrix<T>& m, size_type i_comp) const
 
367
{
 
368
  check_macro (get_first_basis().family() == element_constant::Lagrange &&
 
369
              get_second_basis().family() == element_constant::Lagrange, 
 
370
        "Only Lagrange elements supported");
 
371
  size_type map_d = K.dimension();
 
372
  check_macro (i_comp < map_d, "unexpected `d_dx" << i_comp << "' on element type "
 
373
        << "`" << K.name() << "'");
 
374
  check_macro (coordinate_dimension() == K.dimension(),
 
375
    "unsupported `d_dx" << i_comp << "' form on element type `" << K.name()
 
376
    << " in " << coordinate_dimension() << "D geometry");
 
377
  reference_element hat_K = K;
 
378
  m.resize (get_second_basis().size(hat_K), get_first_basis().size(hat_K));
 
379
  m.clear();
 
380
  std::vector<size_type> dis_inod;
 
381
  _omega.dis_inod (K, dis_inod);
 
382
  quadrature::const_iterator first = _quad.begin(hat_K);
 
383
  quadrature::const_iterator last  = _quad.end  (hat_K);
 
384
  tensor DF;
 
385
  tensor invDF;
 
386
  for (size_type q = 0; first != last; first++, q++) {
 
387
    rheolef::jacobian_piola_transformation (_omega, _piola_table, K, dis_inod, q, DF);
 
388
    invDF = pseudo_inverse_jacobian_piola_transformation (DF, coordinate_dimension(), map_d);
 
389
    invDF = inv(DF,map_d);
 
390
    T wq = is_weighted() ? weight(K,q) : T(1);
 
391
    wq *= weight_coordinate_system (K, dis_inod, q);
 
392
    wq *= rheolef::det_jacobian_piola_transformation (DF, coordinate_dimension(), K.dimension());
 
393
    wq *= (*first).w;
 
394
    point invDFt_icomp_wq = invDF.col(i_comp) * wq;
 
395
  
 
396
    // put in a subroutine:
 
397
    basis_on_quadrature::const_iterator phi      = _by_table.begin(hat_K, q);
 
398
    basis_on_quadrature::const_iterator last_phi = _by_table.end(hat_K, q);
 
399
    basis_on_quadrature::const_iterator_grad first_grad_psi = _bx_table.begin_grad (hat_K, q);
 
400
    basis_on_quadrature::const_iterator_grad last_grad_psi  = _bx_table.end_grad (hat_K, q);
 
401
    for (size_type i = 0; phi != last_phi; phi++, i++) {
 
402
      basis_on_quadrature::const_iterator_grad grad_psi = first_grad_psi;
 
403
      T phi_i = *phi;
 
404
      for (size_type j = 0; grad_psi != last_grad_psi; grad_psi++, j++) {
 
405
        T sum = 0;
 
406
        const point& grad_psi_j = *grad_psi;
 
407
        for (size_type k = 0; k < map_d; k++) {
 
408
          sum += phi_i * invDFt_icomp_wq[k] * grad_psi_j[k];
 
409
        }
 
410
        m(i,j) += sum;
 
411
      }
 
412
    }
 
413
  }
 
414
}
 
415
// --------------------------------------------------------------------
 
416
// int_K (gradt u : grad v) w dK = int_K [(M U).V] w dK,
 
417
//              with u = Sum_j,l U_{nx*l+j} psi_j
 
418
//                       v = Sum_i,k V_{ny*k+i} phi_i
 
419
// so M_{ny*k+i,nx*l+j} = int_K d psi_j/dk * d phi_i/dl w dK 
 
420
// --------------------------------------------------------------------
 
421
template<class T, class M>
 
422
void 
 
423
form_element_rep<T,M>::build_gradt_grad (const geo_element& K, ublas::matrix<T>& m) const
 
424
{
 
425
  check_macro (get_first_basis().family() == element_constant::Lagrange &&
 
426
              get_second_basis().family() == element_constant::Lagrange, 
 
427
        "Only Lagrange elements supported");
 
428
  size_type map_d = K.dimension();
 
429
  check_macro (coordinate_dimension() == K.dimension(),
 
430
    "unsupported `gradt_grad' form on element type `" << K.name()
 
431
    << " in " << coordinate_dimension() << "D geometry");
 
432
  size_type nx = get_first_basis().size(K);
 
433
  size_type ny = get_second_basis().size(K);
 
434
  m.resize (map_d*ny, map_d*nx);
 
435
  m.clear();
 
436
  std::vector<size_type> dis_inod;
 
437
  _omega.dis_inod (K, dis_inod);
 
438
  quadrature::const_iterator first = _quad.begin(K);
 
439
  quadrature::const_iterator last  = _quad.end  (K);
 
440
  tensor DF;
 
441
  tensor invDF;
 
442
  tensor invDFt;
 
443
  tensor Dw;
 
444
  for (size_type q = 0; first != last; first++, q++) {
 
445
    rheolef::jacobian_piola_transformation (_omega, _piola_table, K, dis_inod, q, DF);
 
446
    invDF = inv(DF,map_d);
 
447
    invDFt = trans(invDF);
 
448
    T wq = is_weighted() ? weight(K,q) : T(1);
 
449
    wq *= weight_coordinate_system (K, dis_inod, q);
 
450
    wq *= rheolef::det_jacobian_piola_transformation (DF, coordinate_dimension(), K.dimension());
 
451
    wq *= (*first).w;
 
452
    
 
453
    basis_on_quadrature::const_iterator_grad grad_phi       = _by_table.begin_grad(K, q);
 
454
    basis_on_quadrature::const_iterator_grad last_grad_phi  = _by_table.end_grad(K, q);
 
455
    basis_on_quadrature::const_iterator_grad first_grad_psi = _bx_table.begin_grad (K, q);
 
456
    basis_on_quadrature::const_iterator_grad last_grad_psi  = _bx_table.end_grad (K, q);
 
457
    for (size_type i = 0; grad_phi != last_grad_phi; grad_phi++, i++) {
 
458
      const point& hat_grad_phi_i = *grad_phi;
 
459
      point grad_phi_i = invDFt*hat_grad_phi_i;
 
460
      basis_on_quadrature::const_iterator_grad grad_psi = first_grad_psi;
 
461
      for (size_type j = 0; grad_psi != last_grad_psi; grad_psi++, j++) {
 
462
        const point& hat_grad_psi_j = *grad_psi;
 
463
        point grad_psi_j = invDFt*hat_grad_psi_j;
 
464
        for (size_type k = 0; k < map_d; k++) {
 
465
          for (size_type l = 0; l < map_d; l++) {
 
466
            m(i+ny*k,j+nx*l) += grad_psi_j[k]*grad_phi_i[l]*wq;
 
467
          }
 
468
        }
 
469
      }
 
470
    }
 
471
  }
 
472
}
 
473
// --------------------------------------------------------------------
 
474
// vector-valued : <div_div phi, psi> = (div phi, div psi)
 
475
// Mij += trans(grad_phi_i) * Dw * grad_psi_j
 
476
// --------------------------------------------------------------------
 
477
template<class T, class M>
 
478
void 
 
479
form_element_rep<T,M>::build_div_div (const geo_element& K, ublas::matrix<T>& m) const
 
480
{
 
481
  check_macro (get_first_basis().family() == element_constant::Lagrange &&
 
482
              get_second_basis().family() == element_constant::Lagrange, 
 
483
        "Only Lagrange elements supported");
 
484
  size_type map_d = K.dimension();
 
485
  reference_element hat_K = K;
 
486
  check_macro (coordinate_dimension() == K.dimension(),
 
487
    "unsupported `gradt_grad' form on element type `" << K.name()
 
488
    << "' in " << coordinate_dimension() << "D geometry");
 
489
  size_type nx = get_first_basis().size(hat_K);
 
490
  size_type ny = get_second_basis().size(hat_K);
 
491
  m.resize (map_d*ny, map_d*nx);
 
492
  m.clear();
 
493
  std::vector<size_type> dis_inod;
 
494
  _omega.dis_inod (K, dis_inod);
 
495
  quadrature::const_iterator first = _quad.begin(hat_K);
 
496
  quadrature::const_iterator last  = _quad.end  (hat_K);
 
497
  tensor DF;
 
498
  tensor invDF;
 
499
  tensor invDFt;
 
500
  tensor Dw;
 
501
  for (size_type q = 0; first != last; first++, q++) {
 
502
    rheolef::jacobian_piola_transformation (_omega, _piola_table, K, dis_inod, q, DF);
 
503
    invDF = inv(DF,map_d);
 
504
    invDFt = trans(invDF);
 
505
    T wq = is_weighted() ? weight(K,q) : T(1);
 
506
    wq *= weight_coordinate_system (K, dis_inod, q);
 
507
    wq *= rheolef::det_jacobian_piola_transformation (DF, coordinate_dimension(), K.dimension());
 
508
    wq *= (*first).w;
 
509
    
 
510
    basis_on_quadrature::const_iterator_grad grad_phi       = _by_table.begin_grad(hat_K, q);
 
511
    basis_on_quadrature::const_iterator_grad last_grad_phi  = _by_table.end_grad(hat_K, q);
 
512
    basis_on_quadrature::const_iterator_grad first_grad_psi = _bx_table.begin_grad (hat_K, q);
 
513
    basis_on_quadrature::const_iterator_grad last_grad_psi  = _bx_table.end_grad (hat_K, q);
 
514
    for (size_type i = 0; grad_phi != last_grad_phi; grad_phi++, i++) {
 
515
      const point& hat_grad_phi_i = *grad_phi;
 
516
      point grad_phi_i = invDFt*hat_grad_phi_i;
 
517
      basis_on_quadrature::const_iterator_grad grad_psi = first_grad_psi;
 
518
      for (size_type j = 0; grad_psi != last_grad_psi; grad_psi++, j++) {
 
519
        const point& hat_grad_psi_j = *grad_psi;
 
520
        point grad_psi_j = invDFt*hat_grad_psi_j;
 
521
        for (size_type k = 0; k < map_d; k++) {
 
522
          for (size_type l = 0; l < map_d; l++) {
 
523
            m(i+ny*k,j+nx*l) += grad_psi_j[l]*grad_phi_i[k]*wq;
 
524
          }
 
525
        }
 
526
      }
 
527
    }
 
528
  }
 
529
}
 
530
// ----------------------------------------------------------
 
531
// elementary grad_grad form:
 
532
// --------------------------------------------------------------------
 
533
// weighted forms: i-th component on node x_q in K
 
534
template<class T, class M>
 
535
void
 
536
form_element_rep<T,M>::set_weight (const field_basic<T,M>& wh) const
 
537
{
 
538
  _wh = wh;
 
539
  _wh.dis_dof_update();
 
540
  _bw_table.set (_quad, _wh.get_space().get_numbering().get_basis());
 
541
  _is_weighted = true;
 
542
  _qopt.set_order (0); // will be re-computed at re-initailization
 
543
  initialize_all ();
 
544
}
 
545
template<class T, class M>
 
546
T
 
547
form_element_rep<T,M>::weight (const geo_element& K, size_type q) const
 
548
{
 
549
    // get wh value at idof elements
 
550
    const basis& bw = _wh.get_space().get_numbering().get_basis();
 
551
    size_t loc_ndof = bw.size (K.variant()) ;
 
552
    std::vector<size_type> dis_idof1 (loc_ndof);
 
553
    _wh.get_space().dis_idof (K, dis_idof1);
 
554
 
 
555
    std::vector<T> dof (loc_ndof);
 
556
    for (size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
 
557
      dof [loc_idof] = _wh.dis_dof (dis_idof1 [loc_idof]);
 
558
    }
 
559
    basis_on_quadrature::const_iterator iter_b = _bw_table.begin (K, q);
 
560
    T wq = 0;
 
561
    for (size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++, iter_b++) {
 
562
      wq += dof [loc_idof] * (*iter_b); // sum_i w_coef(i)*hat_phi(x_i)
 
563
    }
 
564
    return wq;
 
565
}
299
566
// ==========================================================================
300
567
// instanciation in library
301
568
// ==========================================================================