← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/form_element/2D_D.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/form_element/2D_D.cc
19
19
///
20
20
/// =========================================================================
21
21
#include "2D_D.h"
 
22
#include "rheolef/tensor4.h"
22
23
#include "rheolef/ublas_matrix_range.h"
23
 
using namespace rheolef;
 
24
#include "map_projector.h"
 
25
 
 
26
namespace rheolef {
24
27
using namespace std;
25
28
using namespace ublas;
26
29
 
27
 
//          [    u0      u1      u2   ]
28
 
//
29
 
// [ v0 ]   [ d0^T*d0 d0^T*d1 d0^T*d2 ]
30
 
// [ v1 ] = [ d1^T*d0 d1^T*d1 d1^T*d2 ]
31
 
// [ v2 ]   [ d2^T*d0 d2^T*d1 d2^T*d2 ]
32
 
//
33
 
//          + (sum_{k=0}^2 dk^T*dk)  Id
34
 
//
35
 
//        = grad^T*grad + (sum_{k=0}^2 dk^T*dk)  Id
36
 
//
37
 
//  In the axisymmetric case (r,z):
38
 
//      2D_D(0,0) += 2 ur/r vr/r r dr dz 
39
 
//                 = 2 ur vr (1/r) dr dz
40
 
//
41
 
void
42
 
_2D_D::operator() (const geo_element& K, matrix<Float>& dd) const
43
 
{
44
 
  build_gradt_grad (K, dd);
 
30
void
 
31
compute_tensor_c (
 
32
    const geo_element& K, 
 
33
    const tensor&      G, 
 
34
    const tensor&      P, 
 
35
    size_t             d,
 
36
    size_t             map_d,
 
37
    tensor4&           C)
 
38
{
 
39
 typedef _2D_D::size_type size_type;
 
40
 C = 0;
 
41
 if (d > map_d) {
 
42
   // ------------------------------------------
 
43
   // cas surfacique
 
44
   // ------------------------------------------
 
45
   // calcul de B 
 
46
   tensor4 B;
 
47
   for (size_type i = 0; i < d; i++)
 
48
   for (size_type j = 0; j < d; j++) {
 
49
       for (size_type k = 0; k < d;     k++)
 
50
       for (size_type r = 0; r < map_d; r++) {
 
51
           Float sum = 0;
 
52
           for (size_type l = 0; l < d; l++) {
 
53
              sum += (P(i,k)*P(l,j) + P(i,l)*P(k,j))*G(l,r);
 
54
           }
 
55
           B(i,j,k,r) = sum;
 
56
        }
 
57
   }
 
58
   // calcul de C cas surfacique  
 
59
   for (size_type k = 0; k < d;     k++)
 
60
   for (size_type r = 0; r < map_d; r++) {
 
61
       for (size_type p = 0; p < d;     p++)
 
62
       for (size_type t = 0; t < map_d; t++) {
 
63
           for (size_type i = 0; i < d; i++)
 
64
           for (size_type j = 0; j < d; j++) {
 
65
               C(i,j,k,r) += B(i,j,k,r)*B(i,j,p,t);
 
66
           }
 
67
       }
 
68
   }
 
69
 } else {
 
70
   // ------------------------------------------
 
71
   //  calcul de C cas classique;
 
72
   // ------------------------------------------
 
73
   for (size_type r = 0; r < map_d; r++) 
 
74
   for (size_type t = 0; t < map_d; t++) {
 
75
     for (size_type k = 0; k < d; k++) {
 
76
       Float sum = 0;
 
77
       for (size_type i = 0; i < d; i++) {
 
78
         sum += 2*G(i,r)*G(i,t);
 
79
       }
 
80
       C(k,r,k,t) += sum;
 
81
     }
 
82
     for (size_type k = 0; k < d; k++)
 
83
     for (size_type p = 0; p < d; p++) {
 
84
       for (size_type i = 0; i < d; i++) {
 
85
         C(k,r,p,t) += 2*G(p,r)*G(k,t);
 
86
       }
 
87
     }
 
88
   }
 
89
 }
 
90
}
 
91
void
 
92
_2D_D::operator() (const geo_element& K, matrix<Float>& ak) const
 
93
{
45
94
  size_type nj = get_first_basis().size(K);
46
95
  size_type ni = get_second_basis().size(K);
 
96
#undef USE_NEW_VERSION
 
97
#ifdef USE_NEW_VERSION
 
98
  size_type d     = coordinate_dimension();
 
99
  size_type map_d = K.dimension();
 
100
  reference_element hat_K = K;
 
101
  size_type nx = get_first_basis().size(hat_K);
 
102
  size_type ny = get_second_basis().size(hat_K);
 
103
  ak.resize (d*ny, d*nx);
 
104
  ak.clear();
 
105
  quadrature::const_iterator first = _quad.begin(hat_K);
 
106
  quadrature::const_iterator last  = _quad.end  (hat_K);
 
107
  tensor DF;
 
108
  tensor invDF;
 
109
  tensor G;
 
110
  tensor P; // P = I - nxn
 
111
  tensor4 C;
 
112
  for (size_type q = 0; first != last; first++, q++) {
 
113
    jacobian_piola_transformation (K,q,DF);
 
114
    invDF = pseudo_inverse_jacobian_piola_transformation (DF, map_d);
 
115
    Float wq = is_weighted() ? weight(K,q) : Float(1);
 
116
    wq *= weight_coordinate_system (K, q);
 
117
    wq *= det_jacobian_piola_transformation (DF, K.dimension());
 
118
    wq *= (*first).w;
 
119
    jacobian_piola_transformation (K,q,DF);
 
120
    G = trans(invDF);
 
121
    map_projector (DF, d, map_d, P);
 
122
    compute_tensor_c (K, G, P, d, map_d, C);
 
123
    basis_on_quadrature::const_iterator_grad grad_phi       = _by_table.begin_grad(hat_K, q);
 
124
    basis_on_quadrature::const_iterator_grad last_grad_phi  = _by_table.end_grad(hat_K, q);
 
125
    basis_on_quadrature::const_iterator_grad first_grad_psi = _bx_table.begin_grad (hat_K, q);
 
126
    basis_on_quadrature::const_iterator_grad last_grad_psi  = _bx_table.end_grad (hat_K, q);
 
127
    for (size_type i = 0; grad_phi != last_grad_phi; grad_phi++, i++) {
 
128
      const point& grad_phi_i = *grad_phi;
 
129
      basis_on_quadrature::const_iterator_grad grad_psi = first_grad_psi;
 
130
      for (size_type j = 0; grad_psi != last_grad_psi; grad_psi++, j++) {
 
131
        const point& grad_psi_j = *grad_psi;
 
132
        for (size_type k = 0; k < d; k++)
 
133
        for (size_type p = 0; p < d; p++)
 
134
        for (size_type r = 0; r < map_d; r++)
 
135
        for (size_type t = 0; t < map_d; t++) {
 
136
          ak(i+k*ny,j+p*nx) += C(k,r,p,t)*wq*grad_phi_i[r]*grad_psi_j[t];
 
137
        }
 
138
      }
 
139
    }
 
140
  }
 
141
#else // USE_OLD_VERSION
 
142
  build_gradt_grad (K, ak);
47
143
  switch (K.dimension()) {
48
144
    case 1:
49
 
      dd *= 2;
 
145
      ak *= 2;
50
146
      break;
51
147
    case 2:
52
148
    case 3: {
53
149
      matrix<Float> trace (ni, nj);
54
150
      trace.clear();
55
151
      for (size_type k = 0; k < K.dimension(); k++) {
56
 
        matrix_range<matrix<Float> > dk (dd, range(k*ni,(k+1)*ni),range(k*nj,(k+1)*nj));
 
152
        matrix_range<matrix<Float> > dk (ak, range(k*ni,(k+1)*ni),range(k*nj,(k+1)*nj));
57
153
        trace += dk;
58
154
      }
59
155
      for (size_type k = 0; k < K.dimension(); k++) {
60
 
        matrix_range<matrix<Float> > dk (dd, range(k*ni,(k+1)*ni), range(k*nj,(k+1)*nj));
 
156
        matrix_range<matrix<Float> > dk (ak, range(k*ni,(k+1)*ni), range(k*nj,(k+1)*nj));
61
157
        dk += trace;
62
158
     }
63
159
     break;
67
163
      break;
68
164
    }
69
165
  }
 
166
#endif // USE_OLD_VERSION
70
167
  if (coordinate_system_type() == fem_helper::cartesian) return;
71
168
  // --------------------------------------
72
169
  // otherwise, add the (ur/r) (vr/r) r dr dz term
79
176
  set_use_coordinate_system_dual_weight (false);
80
177
  size_type k = 0;
81
178
  if (coordinate_system_type() == fem_helper::axisymmetric_zr) k++;
82
 
  matrix_range<matrix<Float> > dk (dd, range(k*ni,(k+1)*ni), range(k*nj,(k+1)*nj));
 
179
  matrix_range<matrix<Float> > dk (ak, range(k*ni,(k+1)*ni), range(k*nj,(k+1)*nj));
83
180
  dk += 2.*m_invr;
84
181
}
 
182
 
85
183
_2D_D::size_type
86
184
_2D_D::n_derivative() const
87
185
{
88
186
  if (coordinate_system_type() == fem_helper::cartesian) return 2;
89
187
  else return 0;
90
188
}
 
189
 
91
190
void
92
191
_2D_D::check_after_initialize () const
93
192
{
97
196
        d == get_second_space().n_component(),
98
197
        "unsupported non-vectorial space for `2D_D' form");
99
198
}
 
199
 
 
200
} // namespace rheolef