← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/form_assembly.icc

  • Committer: Bazaar Package Importer
  • Author(s): Pierre Saramito
  • Date: 2011-03-23 11:14:26 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20110323111426-cjvhey7lxt6077ty
Tags: 5.93-1
* New upstream release (minor changes):
  - some extra warning message deleted in heap_allocator
  - graphic output with mayavi2 fixed
  - add doc refman in .info and .pdf format

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/plib/form_assembly.icc
27
27
//
28
28
# include "rheolef/form.h"
29
29
namespace rheolef { 
30
 
using namespace std;
31
30
 
32
 
template <class T>
 
31
template <class T, class M>
33
32
void
34
 
form_basic<T>::assembly (
35
 
    const geo_type&          omega,
36
 
    const form_element<T>&   form_e) 
 
33
form_basic<T,M>::assembly (
 
34
    const form_element<T,M>&   form_e,
 
35
    const geo_basic<T,M>&      X_geo,
 
36
    const geo_basic<T,M>&      Y_geo,
 
37
    bool                       X_geo_is_background)
37
38
{
 
39
    using namespace std;
38
40
    //
39
41
    // use temporary asr = dynamic matrix structures
40
42
    // 
41
 
    asr<T> auu (_Y.iu_ownership(), _X.iu_ownership());
42
 
    asr<T> aub (_Y.iu_ownership(), _X.ib_ownership());
43
 
    asr<T> abu (_Y.ib_ownership(), _X.iu_ownership());
44
 
    asr<T> abb (_Y.ib_ownership(), _X.ib_ownership());
45
 
    ublas::matrix<Float> ak;
46
 
    std::vector<size_type> par_idy, par_jdx;
 
43
    asr<T,M> auu (_Y.iu_ownership(), _X.iu_ownership());
 
44
    asr<T,M> aub (_Y.iu_ownership(), _X.ib_ownership());
 
45
    asr<T,M> abu (_Y.ib_ownership(), _X.iu_ownership());
 
46
    asr<T,M> abb (_Y.ib_ownership(), _X.ib_ownership());
 
47
    ublas::matrix<T> ak;
 
48
    std::vector<size_type> dis_idy, dis_jdx;
47
49
 
48
 
    for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
49
 
      const geo_element& K = omega [ie];
 
50
    for (size_type ie = 0, ne = X_geo.size(); ie < ne; ie++) {
 
51
      const geo_element& Kx = X_geo [ie];
 
52
      const geo_element& Ky = Y_geo [ie];
50
53
      // --------------------------------
51
54
      // compute elementary matrix ak
52
55
      // --------------------------------
53
 
      form_e (K, ak);
54
 
      _X.set_par_dof (K, par_jdx);
55
 
      _Y.set_par_dof (K, par_idy);
 
56
      if (X_geo_is_background) {
 
57
          form_e (Kx, ak);
 
58
      } else {
 
59
          form_e (Ky, ak); 
 
60
      }
 
61
      _X.get_dis_idof (Kx, dis_jdx);
 
62
      _Y.get_dis_idof (Ky, dis_idy);
56
63
      // ------------------------------
57
64
      // assembly ak in sparse matrix a
58
65
      // ------------------------------
60
67
        for (size_type loc_jdof = 0, nx = ak.size2(); loc_jdof < nx; loc_jdof++) {
61
68
 
62
69
            const T& value = ak (loc_idof, loc_jdof);
63
 
            if (value == Float(0)) continue;
64
 
 
65
 
            size_type par_idof = par_idy [loc_idof];
66
 
            size_type par_jdof = par_jdx [loc_jdof];
67
 
 
68
 
            size_type par_iub = _Y.par_iub (par_idof);
69
 
            size_type par_jub = _X.par_iub (par_jdof);
70
 
 
71
 
            if   (_Y.par_is_blocked(par_idof))
72
 
              if (_X.par_is_blocked(par_jdof)) abb.entry (par_iub, par_jub) += value;
73
 
              else                             abu.entry (par_iub, par_jub) += value;
 
70
#ifdef TO_CLEAN
 
71
            // not performed: conserve the sparsity pattern, even with some zero coefs
 
72
            if (value == T(0)) continue;
 
73
#endif // TO_CLEAN
 
74
 
 
75
            size_type dis_idof = dis_idy [loc_idof];
 
76
            size_type dis_jdof = dis_jdx [loc_jdof];
 
77
 
 
78
            size_type dis_iub = _Y.dis_iub (dis_idof);
 
79
            size_type dis_jub = _X.dis_iub (dis_jdof);
 
80
 
 
81
            if   (_Y.dis_is_blocked(dis_idof))
 
82
              if (_X.dis_is_blocked(dis_jdof)) abb.dis_entry (dis_iub, dis_jub) += value;
 
83
              else                             abu.dis_entry (dis_iub, dis_jub) += value;
74
84
            else 
75
 
              if (_X.par_is_blocked(par_jdof)) aub.entry (par_iub, par_jub) += value;
76
 
              else                             auu.entry (par_iub, par_jub) += value;
 
85
              if (_X.dis_is_blocked(dis_jdof)) aub.dis_entry (dis_iub, dis_jub) += value;
 
86
              else                             auu.dis_entry (dis_iub, dis_jub) += value;
77
87
        }
78
88
      }
79
89
    }
80
90
    //
81
 
    // since all is local, axx.assembly() compute only axx.par_nnz
 
91
    // since all is local, axx.dis_entry_assembly() compute only axx.dis_nnz
82
92
    //
83
 
    auu.assembly();
84
 
    aub.assembly();
85
 
    abu.assembly();
86
 
    abb.assembly();
 
93
    auu.dis_entry_assembly();
 
94
    aub.dis_entry_assembly();
 
95
    abu.dis_entry_assembly();
 
96
    abb.dis_entry_assembly();
87
97
    //
88
98
    // convert dynamic matrix asr to fixed-size one csr
89
99
    //
90
 
    _uu = csr<Float>(auu);
91
 
    _ub = csr<Float>(aub);
92
 
    _bu = csr<Float>(abu);
93
 
    _bb = csr<Float>(abb);
 
100
    uu = csr<T,M>(auu);
 
101
    ub = csr<T,M>(aub);
 
102
    bu = csr<T,M>(abu);
 
103
    bb = csr<T,M>(abb);
 
104
    //
 
105
    // set pattern dimension to uu:
 
106
    // => used by solvers, for efficiency: direct(d<3) or iterative(d=3)
 
107
    //
 
108
    uu.set_pattern_dimension (X_geo.map_dimension());
 
109
    //
 
110
    // when form_e is symmetric, uu & bb are also
 
111
    // => used by solvers, for efficiency: LDL^t or LU, CG or GMRES
 
112
    //
 
113
    bool is_symm = form_e.is_symmetric();
 
114
    uu.set_symmetry (is_symm);
 
115
    bb.set_symmetry (is_symm);
94
116
}
95
117
 
96
118
}// namespace rheolef