← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/space_seq.cc

  • 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/space_seq.cc
24
24
namespace rheolef {
25
25
 
26
26
// ---------------------------------------------------------------
27
 
// space_rep
 
27
// space_base_rep
28
28
// ---------------------------------------------------------------
29
29
template <class T, class M>
30
 
inline
31
 
space_rep<T,M>::space_rep (
32
 
    const geo_basic<T,M>& omega,
 
30
space_base_rep<T,M>::space_base_rep (
 
31
    const geo_basic<T,M>& omega_in,
33
32
    std::string approx)
34
 
  : _constit(omega,approx), 
 
33
  : _constit(omega_in,approx), 
35
34
    _is_freezed(false),
36
35
    _idof2blk_iub(),
37
36
    _iu_ownership(),
38
37
    _ib_ownership()
39
38
{
 
39
    // omega_in is compressed by space_constitution() allocator
 
40
    // when it is a domain: then it becomes a geo_domain, with compressed numbering
 
41
    // so, omega_in can be different from omega:
 
42
    const geo_basic<T,M>& omega = get_geo();
 
43
 
40
44
    if (approx == "") return; // empty element => default cstor
41
 
    // implements only P1 without blocked domains:
42
 
    // P1 ownership follows vertex ownership
43
 
    warning_macro ("cstor: par_n_vertex  = " << get_geo().par_n_vertex() << ", n_vertex  = " << get_geo().n_vertex());
44
 
    size_type par_ndof = get_geo().par_n_vertex();
45
 
    size_type ndof     = get_geo().n_vertex();
46
 
    distributor idof_ownership (par_ndof, comm(), ndof);
47
 
    warning_macro ("cstor: 0) par_ndof  = " << par_ndof << ", ndof  = " << ndof);
48
 
    warning_macro ("cstor: 1) par_size  = " << idof_ownership.par_size() << ", size  = " << idof_ownership.size());
 
45
    size_type ndof     = _constit.get_element().get_numbering().ndof (
 
46
                                get_geo().map_dimension(),
 
47
                                get_geo().size_by_dimension(),
 
48
                                get_geo().size_by_variant());
 
49
    size_type dis_ndof = _constit.get_element().get_numbering().ndof (
 
50
                                get_geo().map_dimension(),
 
51
                                get_geo().dis_size_by_dimension(),
 
52
                                get_geo().dis_size_by_variant());
 
53
 
 
54
    distributor idof_ownership (dis_ndof, comm(), ndof);
49
55
    _idof2blk_iub.resize (idof_ownership);
50
 
    warning_macro ("cstor: 2) par_size  = " << ownership().par_size() << ", size  = " << ownership().size());
51
 
    check_macro (idof_ownership.par_size() == par_ndof, "invalid par_ndof count");
 
56
    check_macro (idof_ownership.dis_size() == dis_ndof, "invalid dis_ndof count");
52
57
    check_macro (idof_ownership.size()     == ndof,     "invalid ndof count");
53
58
}
54
59
template <class T, class M>
55
60
void
56
 
space_rep<T,M>::base_freeze_body () const
 
61
space_base_rep<T,M>::base_freeze_body () const
57
62
{
58
 
    warning_macro ("base_freeze_body...");
59
63
    check_macro (get_element().name() != "", "space with undefined finite element method cannot be used");
60
64
    // -----------------------------------------------------------------------
61
65
    // loop on domains: mark blocked dofs
62
66
    // -----------------------------------------------------------------------
63
 
    warning_macro ("base_freeze_body[1]: n_dom="<<_constit.size());
64
67
    const geo_basic<T,M>& omega = get_geo();
65
 
    size_type first_par_idof = ownership().first_index();
66
 
    size_type last_par_idof  = ownership().last_index();
67
 
    warning_macro ("first_par_idof = " << first_par_idof);
68
 
    warning_macro ("last_par_idof  = " << last_par_idof);
 
68
    size_type first_dis_idof = ownership().first_index();
 
69
    size_type last_dis_idof  = ownership().last_index();
69
70
    array<size_type,M> blocked_flag (ownership(), 0); // array<bool,M> not supported
70
71
    for (typename space_constitution<T,M>::const_iterator
71
72
        iter = _constit.begin(),
72
73
        last = _constit.end(); iter != last; iter++) {
73
 
      const domain_basic<T,M>&          dom = (*iter).get_domain();
74
 
      typename space_act<T,M>::act_type act = (*iter).get_act();
75
 
      size_type dom_dim = dom.dimension();
76
 
      distributor isg_ownership = omega.subgeo_ownership (dom_dim);
77
 
      size_type   first_par_isg = isg_ownership.first_index();
78
 
      size_type    last_par_isg = isg_ownership.last_index();
79
 
      bool blk = (act == space_act<T,M>::block) ? true : false;
80
 
      warning_macro ("block domain["<<dom.name()<<"]...");
81
 
      std::vector<size_type> par_idof;
 
74
      const domain_indirect_basic<M>&     dom = (*iter).get_domain_indirect();
 
75
      typename space_act<M>::act_type act = (*iter).get_act();
 
76
      size_type dom_dim = dom.map_dimension();
 
77
      distributor ige_ownership = omega.geo_element_ownership (dom_dim);
 
78
      size_type   first_dis_ige = ige_ownership.first_index();
 
79
      size_type    last_dis_ige = ige_ownership.last_index();
 
80
      bool blk = (act == space_act<M>::block) ? true : false;
 
81
      std::vector<size_type> dis_idof;
82
82
      geo_element_p P;
83
 
      for (size_type ioisg = 0, noisg = dom.size(); ioisg < noisg; ioisg++) {
84
 
        size_type par_isg = dom.oisg (ioisg).index();
85
 
        check_macro (par_isg >= first_par_isg && par_isg < last_par_isg, "par_isg="<<par_isg<<" out of range");
86
 
        // TODO: shift par_isg to isg into domain
87
 
        size_type isg = par_isg - first_par_isg;
 
83
      for (size_type ioige = 0, noige = dom.size(); ioige < noige; ioige++) {
 
84
        size_type ige     = dom.oige (ioige).index();
 
85
        size_type dis_ige = ige + first_dis_ige;
 
86
#ifdef TO_CLEAN
88
87
        if (dom_dim > 0) {
89
 
          const geo_element& S = omega.subgeo (dom_dim, isg);
90
 
          warning_macro ("block domain["<<dom.name()<<"].side["<<ioisg<<"] = element["<<par_isg<<"] = " << S);
91
 
          set_par_dof (S, par_idof);
 
88
#endif // TO_CLEAN
 
89
          const geo_element& S = omega.get_geo_element (dom_dim, ige);
 
90
          get_dis_idof (S, dis_idof);
 
91
#ifdef TO_CLEAN
92
92
        } else {
93
93
          // fix for the d=0 domains of vertices:
94
 
          P[0] = par_isg;
95
 
          set_par_dof (P, par_idof);
 
94
          P[0] = dis_ige;
 
95
          get_dis_idof (P, dis_idof);
96
96
        }
97
 
        warning_macro ("block domain["<<dom.name()<<"].side["<<ioisg<<"].nloc = " << par_idof.size());
98
 
        for (size_type iloc = 0, nloc = par_idof.size(); iloc < nloc; iloc++) {
99
 
          if (par_idof[iloc] >= first_par_idof && par_idof[iloc] < last_par_idof) {
100
 
            size_type idof = par_idof [iloc] - first_par_idof;
101
 
            warning_macro ("block domain["<<dom.name()<<"].side["<<ioisg<<"]["<<iloc<<"]: block local par_idof="<<par_idof[iloc]);
 
97
#endif // TO_CLEAN
 
98
        for (size_type iloc = 0, nloc = dis_idof.size(); iloc < nloc; iloc++) {
 
99
          if (dis_idof[iloc] >= first_dis_idof && dis_idof[iloc] < last_dis_idof) {
 
100
            size_type idof = dis_idof [iloc] - first_dis_idof;
102
101
            blocked_flag[idof] = blk;
103
102
          } else {
104
 
            warning_macro ("block domain["<<dom.name()<<"].side["<<ioisg<<"]["<<iloc<<"]: block NONLOCAL par_idof="<<par_idof[iloc]);
105
 
            blocked_flag.set (par_idof[iloc], blk);
 
103
            blocked_flag.dis_entry (dis_idof[iloc]) = blk;
106
104
          }
107
105
        }
108
106
      }
109
107
    }
110
 
    blocked_flag.assembly();
 
108
    blocked_flag.dis_entry_assembly();
111
109
    // copy the blocked_flag into the iub array, as the "blocked" bit:
112
110
    for (size_type idof = 0, ndof = blocked_flag.size(); idof < ndof; idof++) {
113
111
      _idof2blk_iub [idof].set_blocked (blocked_flag[idof]);
115
113
    // -----------------------------------------------------------------------
116
114
    // init numbering
117
115
    // -----------------------------------------------------------------------
118
 
    warning_macro ("base_freeze_body[2]: init numbering...");
119
116
    size_type n_unknown = 0;
120
117
    size_type n_blocked = 0;
121
118
    for (size_type idof = 0, ndof = _idof2blk_iub.size(); idof < ndof; idof++) {
128
125
          _idof2blk_iub[idof].set_iub (n_blocked);
129
126
          n_blocked++;
130
127
        }
131
 
        warning_macro ("par_idof2blk_iub["<<first_par_idof + idof<<"] = " << _idof2blk_iub[idof]);
132
 
    }
133
 
    size_type par_n_unknown = mpi::all_reduce (comm(), n_unknown, std::plus<T>());
134
 
    size_type par_n_blocked = mpi::all_reduce (comm(), n_blocked, std::plus<T>());
135
 
    _iu_ownership = distributor (par_n_unknown, comm(), n_unknown);
136
 
    _ib_ownership = distributor (par_n_blocked, comm(), n_blocked);
137
 
 
138
 
#ifndef TO_CLEAN
139
 
    warning_macro ("par_n  = " << ownership().par_size() << ", n  = " << ownership().size());
140
 
    warning_macro ("par_nu = " << par_n_unknown << ", nu = " << n_unknown);
141
 
    warning_macro ("par_nb = " << par_n_blocked << ", nb = " << n_blocked);
142
 
    
143
 
    check_macro (ownership().par_size() == _iu_ownership.par_size() + _ib_ownership.par_size(),
144
 
                "invalid par count: "
145
 
                <<   ownership().par_size() << " = "
146
 
                << _iu_ownership.par_size() << " + "
147
 
                << _ib_ownership.par_size());
148
 
 
149
 
    check_macro (ownership().size() == _iu_ownership.size() + _ib_ownership.size(),
150
 
                "invalid count: "
151
 
                <<   ownership().size() << " = "
152
 
                << _iu_ownership.size() << " + "
153
 
                << _ib_ownership.size());
154
 
 
155
 
#endif // TO_CLEAN
156
 
 
157
 
    warning_macro ("base_freeze_body done");
 
128
    }
 
129
    size_type dis_n_unknown = n_unknown;
 
130
    size_type dis_n_blocked = n_blocked;
 
131
#ifdef _RHEOLEF_HAVE_MPI
 
132
    if (is_distributed<M>::value) {
 
133
      dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
 
134
      dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
 
135
    }
 
136
#endif // // _RHEOLEF_HAVE_MPI
 
137
    _iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
 
138
    _ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
158
139
}
159
140
// ----------------------------------------------------------------------------
160
141
// used by field & form_assembly
161
142
// ----------------------------------------------------------------------------
162
143
template <class T, class M>
163
144
void
164
 
space_rep<T,M>::set_par_dof (const geo_element& K, std::vector<size_type>& par_idof) const
 
145
space_base_rep<T,M>::get_dis_idof (const geo_element& K, std::vector<size_type>& dis_idof) const
165
146
{
166
 
    // TODO: only for P1 element here
167
147
    freeze_guard();
168
 
    par_idof.resize (K.size());
169
 
    for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
170
 
        par_idof [iloc] = K[iloc];
 
148
    size_type nloc = _constit.get_element().get_basis().size(K.variant());
 
149
    dis_idof.resize (nloc);
 
150
    for (size_type iloc = 0; iloc < nloc; iloc++) {
 
151
      dis_idof [iloc] = _constit.get_element().get_numbering().idof (
 
152
                get_geo().dis_size_by_dimension(),
 
153
                get_geo().dis_size_by_variant(),
 
154
                K,
 
155
                iloc);
171
156
    }
172
157
}
 
158
// ----------------------------------------------------------------------------
 
159
// stamp : e.g. "P1(square)", for field_expr<Expr> checks
 
160
// ----------------------------------------------------------------------------
 
161
template <class T, class M>
 
162
std::string
 
163
space_base_rep<T,M>::get_stamp() const
 
164
{
 
165
    if (get_geo().name() == "*nogeo*") return "";
 
166
    return get_element().name() + "(" + get_geo().name() + ")";
 
167
}
173
168
// ---------------------------------------------------------------
174
 
// space_seq_rep
 
169
// space_rep<seq>
175
170
// ---------------------------------------------------------------
176
171
template <class T>
177
 
space_seq_rep<T>::space_seq_rep (
 
172
space_rep<T,sequential>::space_rep (
178
173
    const geo_basic<T,sequential>& omega,
179
174
    std::string approx)
180
 
  : space_rep<T,sequential>::space_rep (omega, approx)
181
 
{
 
175
  : space_base_rep<T,sequential>::space_base_rep (omega, approx)
 
176
{
 
177
}
 
178
// ----------------------------------------------------------------------------
 
179
// u["left"] 
 
180
// => build here the requiresd temporary indirect array
 
181
// ----------------------------------------------------------------------------
 
182
/**
 
183
 Implementation note: there is two numbering styles
 
184
 
 
185
   idof     : from the current space
 
186
   dom_idof : from a space compacted to the domain "dom", as geo_domain does
 
187
 
 
188
  This function returns the renumbering array "dom_idof2dof".
 
189
  It is a temporary, used at the fly when using the u[dom] syntax.
 
190
 
 
191
  In that case, the geo_domain is not fully builded for efficiency reasons.
 
192
 
 
193
  Nevertheless, the renumbering procedure *MAY* be exactly the same as
 
194
  those of geo_domain, for compatibility purpose when:
 
195
 
 
196
        geo dom = omega["boundary"]);
 
197
        space Vh (omega);
 
198
        space Wh (dom);
 
199
        uh[dom] = gh;
 
200
 
 
201
  For a similar renumbering algorithm, see also geo_rep<T,M>::build_from_domain
 
202
  when numbering the geo_domain Wh.geo vertices.
 
203
 */
 
204
template <class T, class M>
 
205
array<typename space_base_rep<T,M>::size_type, M>
 
206
space_base_rep<T,M>::build_indirect_array (
 
207
    const domain_indirect_basic<M>& indirect) const
 
208
{
 
209
  const geo_basic<T,M>& omega = get_geo();
 
210
  size_type first_dis_idof = ownership().first_index();
 
211
  size_type  last_dis_idof = ownership().last_index();
 
212
  array<size_type> idof_is_on_domain (ownership(), 0); // logical, init to "false"
 
213
  size_type map_dim = indirect.map_dimension();
 
214
  std::vector<size_type> dis_idof;
 
215
  //
 
216
  // 1) loop on domain elements and mark used idofs
 
217
  //
 
218
  for (size_type ioige = 0, noige = indirect.size(); ioige < noige; ioige++) {
 
219
    size_type ige = indirect.oige (ioige).index();
 
220
    const geo_element& S = omega.get_geo_element (map_dim, ige);
 
221
    get_dis_idof (S, dis_idof);
 
222
    for (size_type iloc = 0, nloc = dis_idof.size(); iloc < nloc; iloc++) {
 
223
      idof_is_on_domain.dis_entry (dis_idof[iloc]) += 1;
 
224
    }
 
225
  }
 
226
  idof_is_on_domain.dis_entry_assembly (set_add_op<size_type,size_type>()); // logical "or"
 
227
  //
 
228
  // 2) count marked idofs
 
229
  //
 
230
  size_type dom_ndof = 0;
 
231
  for (size_type idof = 0, ndof = idof_is_on_domain.size(); idof < ndof; idof++) {
 
232
    if (idof_is_on_domain[idof] != 0) dom_ndof++;
 
233
  }
 
234
  //
 
235
  // 3) numbering dom_idof & permutation: dom_idof --> idof
 
236
  //
 
237
  communicator comm = ownership().comm();
 
238
  distributor dom_idof_ownership (distributor::decide, comm, dom_ndof);
 
239
  array<size_type, M> dom_idof2idof (dom_idof_ownership);
 
240
  for (size_type dom_idof = 0, idof = 0, ndof = idof_is_on_domain.size(); idof < ndof; idof++) {
 
241
    if (idof_is_on_domain[idof] == 0) continue;
 
242
    dom_idof2idof [dom_idof] = idof;
 
243
    dom_idof++;
 
244
  }
 
245
  return dom_idof2idof;
182
246
}
183
247
// ----------------------------------------------------------------------------
184
248
// instanciation in library
185
249
// ----------------------------------------------------------------------------
 
250
template class space_base_rep<Float,sequential>;
186
251
template class space_rep<Float,sequential>;
187
 
template class space_seq_rep<Float>;
188
252
 
189
253
#ifdef _RHEOLEF_HAVE_MPI
190
 
template class space_rep<Float,distributed>;
 
254
template class space_base_rep<Float,distributed>;
191
255
#endif // _RHEOLEF_HAVE_MPI
192
256
 
193
257
} // namespace rheolef