← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/space.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/space.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
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.
 
10
///
 
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.
 
15
///
 
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
 
19
/// 
 
20
/// =========================================================================
 
21
 
 
22
#include "rheolef/space.h"
 
23
#include "rheolef/space_mult.h"
 
24
#include "rheolef/piola.h"
 
25
 
 
26
namespace rheolef {
 
27
 
 
28
// ======================================================================
 
29
// space_base_rep
 
30
// ======================================================================
 
31
template <class T, class M>
 
32
space_base_rep<T,M>::space_base_rep (
 
33
    const geo_basic<T,M>& omega_in,
 
34
    std::string approx,
 
35
    std::string valued)
 
36
  : _constit(omega_in,approx,valued),
 
37
    _xdof(),
 
38
    _have_freezed(false),
 
39
    _idof2blk_iub(),
 
40
    _iu_ownership(),
 
41
    _ib_ownership()
 
42
{
 
43
  // omega_in is compressed by space_scalar_constitution() allocator
 
44
  // when it is a domain: then it becomes a geo_domain, with compressed numbering
 
45
  // so, omega_in can be different from omega:
 
46
  if (approx == "") return; // empty element => default cstor
 
47
  _idof2blk_iub.resize (_constit.ownership());
 
48
  init_xdof();
 
49
}
 
50
template <class T, class M>
 
51
space_base_rep<T,M>::space_base_rep (const space_constitution<T,M>& constit)
 
52
  : _constit(constit),
 
53
    _xdof(),
 
54
    _have_freezed(false),
 
55
    _idof2blk_iub(),
 
56
    _iu_ownership(),
 
57
    _ib_ownership()
 
58
{
 
59
  distributor idof_ownership = _constit.ownership();
 
60
  _idof2blk_iub.resize (idof_ownership);
 
61
  init_xdof();
 
62
}
 
63
template <class T, class M>
 
64
space_base_rep<T,M>::space_base_rep (const space_mult_list<T,M>& expr)
 
65
  : _constit(expr),
 
66
    _xdof(),
 
67
    _have_freezed(false),
 
68
    _idof2blk_iub(),
 
69
    _iu_ownership(),
 
70
    _ib_ownership()
 
71
{
 
72
  _idof2blk_iub.resize (_constit.ownership());
 
73
  init_xdof();
 
74
}
 
75
template <class T, class M>
 
76
void
 
77
space_base_rep<T,M>::init_xdof()
 
78
{
 
79
  if (_constit.valued_tag() == space_constant::mixed) {
 
80
    trace_macro ("init_xdof and mixed valued space: not yet");
 
81
    return;
 
82
  }
 
83
  // precompute one time for all all dof nodes, for an easy interpolation
 
84
  // apply only to nodal fem basis (eg. Lagrange) only
 
85
  // TODO: in a lasy way: at first call to xdof(idof), call this initialization
 
86
  const geo_basic<T,M>& omega = get_geo();
 
87
  const numbering<T,M>& fem   = get_numbering();
 
88
  bool is_isoparametric = (fem.name() == omega.get_piola_basis().name());
 
89
  if (is_isoparametric) {
 
90
    // 1.a) geometric nodes coincides to fem dof nodes:  
 
91
    _xdof = omega.get_nodes();
 
92
    return;
 
93
  }
 
94
  // 1.b) here is the non-isoparametric case: eg. P0 fem space on a P1 geo
 
95
  distributor xdof_ownership;
 
96
  if (_constit.valued_tag() == space_constant::scalar) {
 
97
    xdof_ownership = ownership();
 
98
  } else {
 
99
    // get xdof sizes from any component:
 
100
    xdof_ownership = _constit[0].ownership();
 
101
  }
 
102
  _xdof.resize (xdof_ownership);
 
103
  std::vector<bool> marked (ndof(), false);
 
104
  size_type count_mark = 0;
 
105
  std::vector<size_type> dis_idof;
 
106
  std::vector<size_type> dis_inod;
 
107
  size_type first_dis_idof = xdof_ownership.first_index();
 
108
  basis_on_nodal_basis b (fem.get_basis(), omega.get_piola_basis());
 
109
  for (typename geo_basic<T,M>::const_iterator iter = omega.begin(), last = omega.end(); iter != last; iter++) {
 
110
    const geo_element& K = *iter;
 
111
    fem.dis_idof (omega.sizes(), K, dis_idof);
 
112
    omega.dis_inod (K, dis_inod);
 
113
    for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
 
114
      assert_macro (dis_idof[loc_idof] < xdof_ownership.dis_size(), 
 
115
        "dis_idof " << dis_idof[loc_idof] << " out of range [0:" << xdof_ownership.dis_size() << "[");
 
116
      if (! xdof_ownership.is_owned (dis_idof[loc_idof])) continue;
 
117
      size_type idof = dis_idof[loc_idof] - first_dis_idof;
 
118
      if (marked [idof]) continue;
 
119
      marked [idof] = true;
 
120
      count_mark++;
 
121
      _xdof[idof] = piola_transformation (omega, b, K.variant(), dis_inod, loc_idof);
 
122
    }
 
123
  }
 
124
  size_type dis_n_missing = marked.size() - count_mark;
 
125
#ifdef _RHEOLEF_HAVE_MPI
 
126
  dis_n_missing = mpi::all_reduce (comm(), dis_n_missing, std::plus<size_type>());
 
127
#endif // _RHEOLEF_HAVE_MPI
 
128
  if (dis_n_missing == 0) {
 
129
    return;
 
130
  }
 
131
  //
 
132
  // 2) on some rare 3d bdry geo_domains with distributed nproc > 1, there are missing xdofs...
 
133
  // => second pass: assembly only external xdofs, an performs some comms
 
134
  for (typename geo_basic<T,M>::const_iterator iter = omega.begin(), last = omega.end(); iter != last; iter++) {
 
135
    const geo_element& K = *iter;
 
136
    fem.dis_idof (omega.sizes(), K, dis_idof);
 
137
    omega.dis_inod (K, dis_inod);
 
138
    for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
 
139
      if (xdof_ownership.is_owned (dis_idof[loc_idof])) continue;
 
140
      _xdof.dis_entry (dis_idof[loc_idof]) = piola_transformation (omega, b, K.variant(), dis_inod, loc_idof);
 
141
    }
 
142
  }
 
143
  _xdof.dis_entry_assembly();
 
144
}
 
145
template <class T, class M>
 
146
void
 
147
space_base_rep<T,M>::base_freeze_body () const
 
148
{
 
149
  // -----------------------------------------------------------------------
 
150
  // 1) loop on domains: mark blocked dofs
 
151
  // -----------------------------------------------------------------------
 
152
  array<size_type,M> blocked_flag = _constit.build_blocked_flag();
 
153
 
 
154
  // copy the blocked_flag into the iub array, as the "blocked" bit:
 
155
  for (size_type idof = 0, ndof = blocked_flag.size(); idof < ndof; idof++) {
 
156
    _idof2blk_iub [idof].set_blocked (blocked_flag[idof]);
 
157
  }
 
158
  // -----------------------------------------------------------------------
 
159
  // 2) init numbering
 
160
  // -----------------------------------------------------------------------
 
161
  size_type n_unknown = 0;
 
162
  size_type n_blocked = 0;
 
163
  for (size_type idof = 0, ndof = _idof2blk_iub.size(); idof < ndof; idof++) {
 
164
    bool blk = _idof2blk_iub [idof].is_blocked();
 
165
    if (! blk) {
 
166
      _idof2blk_iub[idof].set_iub (n_unknown);
 
167
      n_unknown++;
 
168
    } else {
 
169
      _idof2blk_iub[idof].set_iub (n_blocked);
 
170
      n_blocked++;
 
171
    }
 
172
  }
 
173
  size_type dis_n_unknown = n_unknown;
 
174
  size_type dis_n_blocked = n_blocked;
 
175
#ifdef _RHEOLEF_HAVE_MPI
 
176
  if (is_distributed<M>::value) {
 
177
    dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
 
178
    dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
 
179
  }
 
180
#endif // // _RHEOLEF_HAVE_MPI
 
181
  _iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
 
182
  _ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
 
183
}
 
184
// ----------------------------------------------------------------------------
 
185
// used by field & form_assembly
 
186
// ----------------------------------------------------------------------------
 
187
template <class T, class M>
 
188
void
 
189
space_base_rep<T,M>::dis_idof (const geo_element& K, std::vector<size_type>& dis_idof) const
 
190
{
 
191
  freeze_guard();
 
192
  _constit.dis_idof (K, dis_idof);
 
193
}
 
194
// ----------------------------------------------------------------------------
 
195
// stamp : e.g. "P1(square)", for field_expr<Expr> checks
 
196
// ----------------------------------------------------------------------------
 
197
template <class T, class M>
 
198
std::string
 
199
space_base_rep<T,M>::stamp() const
 
200
{
 
201
    if (get_geo().name() == "*nogeo*") return "";
 
202
    return get_numbering().name() + "(" + get_geo().name() + ")";
 
203
}
 
204
// ----------------------------------------------------------------------------
 
205
// u["left"] 
 
206
// => build here the requiresd temporary indirect array
 
207
// ----------------------------------------------------------------------------
 
208
/**
 
209
 Implementation note: there is two numbering styles
 
210
 
 
211
   bgd_idof : from the current space
 
212
   dom_idof : from a space compacted to the domain "dom", as geo_domain does
 
213
 
 
214
  This function returns the renumbering array "dom_idof2bgd_idof".
 
215
  It is a temporary, used at the fly when using the u[dom] syntax.
 
216
 */
 
217
template <class T, class M>
 
218
array<typename space_base_rep<T,M>::size_type, M>
 
219
space_base_rep<T,M>::build_indirect_array (
 
220
    const space_base_rep<T,M>& Wh,      // Wh = space on domain
 
221
    const std::string&       dom_name   // redundant: contained in Wh.get_geo()
 
222
  ) const
 
223
{
 
224
  const space_base_rep<T,M>& Vh = *this;  // Vh = space on background mesh
 
225
  Vh.freeze_guard();
 
226
  Wh.freeze_guard();
 
227
  const geo_basic<T,M>& dom_gamma = Wh.get_geo();
 
228
  const geo_basic<T,M>& bgd_gamma = Vh.get_geo()[dom_name];
 
229
  size_type map_dim = dom_gamma.map_dimension();
 
230
  check_macro (dom_gamma.size() == bgd_gamma.size(), "incompatible domains");
 
231
  distributor dom_ownership = Wh.ownership();
 
232
  distributor bgd_ownership = Vh.ownership();
 
233
  size_type first_dom_dis_idof = dom_ownership.first_index();
 
234
  size_type first_bgd_dis_idof = bgd_ownership.first_index();
 
235
  std::vector<size_type> dom_dis_idofs, bgd_dis_idofs;
 
236
  array<size_type, M> dom_idof2bgd_idof (dom_ownership, std::numeric_limits<size_type>::max());
 
237
  for (size_type ige = 0, nge = dom_gamma.size(); ige < nge; ige++) {
 
238
    const geo_element& dom_S = dom_gamma[ige];
 
239
    const geo_element& bgd_S = bgd_gamma[ige];
 
240
    Wh.dis_idof (dom_S, dom_dis_idofs);
 
241
    Vh.dis_idof (bgd_S, bgd_dis_idofs);
 
242
    for (size_type loc_idof = 0, loc_ndof = dom_dis_idofs.size(); loc_idof < loc_ndof; loc_idof++) {
 
243
      size_type dom_dis_idof = dom_dis_idofs [loc_idof];
 
244
      size_type bgd_dis_idof = bgd_dis_idofs [loc_idof];
 
245
      dom_idof2bgd_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
 
246
    }
 
247
  }
 
248
  dom_idof2bgd_idof.dis_entry_assembly();
 
249
  // move to local numbering:
 
250
  for (size_type dom_idof = 0, dom_ndof = dom_idof2bgd_idof.size(); dom_idof < dom_ndof; dom_idof++) {
 
251
      size_type bgd_dis_idof = dom_idof2bgd_idof [dom_idof];
 
252
      size_type dom_dis_idof = dom_idof + first_dom_dis_idof;
 
253
      check_macro (bgd_ownership.is_owned (bgd_dis_idof), "unexpected external dom2bgd("<<dom_dis_idof<<")="<<bgd_dis_idof);
 
254
      size_type bgd_idof = bgd_dis_idof - first_bgd_dis_idof;
 
255
      dom_idof2bgd_idof [dom_idof] = bgd_idof;
 
256
  }
 
257
  return dom_idof2bgd_idof;
 
258
}
 
259
// ----------------------------------------------------------------------------
 
260
// instanciation in library
 
261
// ----------------------------------------------------------------------------
 
262
template class space_base_rep<Float,sequential>;
 
263
 
 
264
#ifdef _RHEOLEF_HAVE_MPI
 
265
template class space_base_rep<Float,distributed>;
 
266
#endif // _RHEOLEF_HAVE_MPI
 
267
 
 
268
} // namespace rheolef