Rheolef  7.1
an efficient C++ finite element environment
space.cc
Go to the documentation of this file.
1 
22 #include "rheolef/space.h"
23 #include "rheolef/space_numbering.h"
24 #include "rheolef/space_mult.h"
25 #include "rheolef/piola_util.h"
26 #include "rheolef/piola_on_pointset.h"
27 
28 namespace rheolef {
29 
30 // ======================================================================
31 // space_base_rep
32 // ======================================================================
33 template <class T, class M>
35  const geo_basic<T,M>& omega_in,
36  std::string approx,
37  std::string valued)
38  : _constit(omega_in,
39  (approx == "" || valued == "scalar") ? approx : valued+"("+approx+")"),
40  _xdof(),
41  _have_freezed(false),
42  _idof2blk_dis_iub(),
43  _has_nt_basis(),
44  _normal(),
45  _iu_ownership(),
46  _ib_ownership()
47 {
48  // omega_in is compressed by space_scalar_constitution() allocator
49  // when it is a domain: then it becomes a geo_domain, with compressed numbering
50  // so, omega_in can be different from omega:
51  if (approx == "") return; // empty element => default cstor
53  _has_nt_basis.resize (_constit.ownership()); // TODO: _constit.scalar_ownership()
54  _normal.resize (_constit.ownership());
55  init_xdof();
56 }
57 template <class T, class M>
59  const geo_basic<T,M>& omega_in,
60  const basis_basic<T>& b)
61  : _constit(omega_in, b.name()),
62  _xdof(),
63  _have_freezed(false),
64  _idof2blk_dis_iub(),
65  _has_nt_basis(),
66  _normal(),
67  _iu_ownership(),
68  _ib_ownership()
69 {
70  // omega_in is compressed by space_scalar_constitution() allocator
71  // when it is a domain: then it becomes a geo_domain, with compressed numbering
72  // so, omega_in can be different from omega:
73  if (! b.is_initialized()) return; // empty element => default cstor
75  _has_nt_basis.resize (_constit.ownership());
76  _normal.resize (_constit.ownership());
77  init_xdof();
78 }
79 template <class T, class M>
81  : _constit(constit),
82  _xdof(),
83  _have_freezed(false),
84  _idof2blk_dis_iub(),
85  _has_nt_basis(),
86  _normal(),
87  _iu_ownership(),
88  _ib_ownership()
89 {
90  distributor idof_ownership = _constit.ownership();
91  _idof2blk_dis_iub.resize (idof_ownership);
92  _has_nt_basis.resize (idof_ownership);
93  _normal.resize (idof_ownership);
94  init_xdof();
95 }
96 template <class T, class M>
98  : _constit(expr),
99  _xdof(),
100  _have_freezed(false),
101  _idof2blk_dis_iub(),
102  _has_nt_basis(),
103  _normal(),
104  _iu_ownership(),
105  _ib_ownership()
106 {
108  _has_nt_basis.resize (_constit.ownership());
109  _normal.resize (_constit.ownership());
110  init_xdof();
111 }
112 template <class T, class M>
113 void
115 {
116  if (_constit.is_hierarchical()) {
117  trace_macro ("init_xdof: hierarchical space: xdof undefined");
118  return;
119  }
120  const geo_basic<T,M>& omega = get_geo();
121  const basis_basic<T>& b = get_basis();
122  size_type dis_nnod = space_numbering::dis_nnod (b, omega.sizes(), omega.map_dimension());
123  size_type nnod = space_numbering:: nnod (b, omega.sizes(), omega.map_dimension());
124  communicator comm = ownership().comm();
125  distributor nod_ownership (dis_nnod, comm, nnod);
126  _xdof.resize (nod_ownership);
128  pops.initialize (omega.get_piola_basis(), b, integrate_option());
129  std::vector<size_type> dis_inod_tab;
130  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
131  const geo_element& K = omega [ie];
132  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
133  space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
134  for (size_type loc_inod = 0, loc_nnod = piola.size(); loc_inod < loc_nnod; ++loc_inod) {
135  _xdof.dis_entry (dis_inod_tab[loc_inod]) = piola [loc_inod].F;
136  }
137  }
138  _xdof.dis_entry_assembly();
139  // add external nodes on current element:
140  if (is_distributed<M>::value && comm.size() > 1) {
141  index_set dis_inod_set;
142  for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
143  const geo_element& K = omega [ie];
144  space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
145  for (size_type loc_inod = 0, loc_nnod = dis_inod_tab.size(); loc_inod < loc_nnod; ++loc_inod) {
146  size_type dis_inod = dis_inod_tab[loc_inod];
147  if (nod_ownership.is_owned(dis_inod)) continue;
148  check_macro (dis_inod < dis_nnod, "invalid dis_inod="<<dis_inod<<" out of range [0:"<<dis_nnod<<"[");
149  dis_inod_set.insert (dis_inod);
150  }
151  }
152  _xdof.set_dis_indexes (dis_inod_set);
153  }
154 }
155 template <class T, class M>
156 void
158 {
159  _constit.neighbour_guard(); // geo could uses S.master(i), a neighbour connnectivity
160  // -----------------------------------------------------------------------
161  // 1) loop on domains: mark blocked dofs
162  // -----------------------------------------------------------------------
163  disarray<size_type,M> blocked_flag = _constit.build_blocked_flag();
164 
165  // copy the blocked_flag into the dis_iub disarray, as the "blocked" bit:
166  for (size_type idof = 0, ndof = blocked_flag.size(); idof < ndof; idof++) {
167  _idof2blk_dis_iub [idof].set_blocked (blocked_flag[idof]);
168  }
169  // -----------------------------------------------------------------------
170  // 2) init numbering
171  // -----------------------------------------------------------------------
172  size_type n_unknown = 0;
173  size_type n_blocked = 0;
174  for (size_type idof = 0, ndof = _idof2blk_dis_iub.size(); idof < ndof; idof++) {
175  bool blk = _idof2blk_dis_iub [idof].is_blocked();
176  if (! blk) {
177  _idof2blk_dis_iub[idof].set_dis_iub (n_unknown);
178  n_unknown++;
179  } else {
180  _idof2blk_dis_iub[idof].set_dis_iub (n_blocked);
181  n_blocked++;
182  }
183  }
184  size_type dis_n_unknown = n_unknown;
185  size_type dis_n_blocked = n_blocked;
186 #ifdef _RHEOLEF_HAVE_MPI
188  dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
189  dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
190  }
191 #endif // // _RHEOLEF_HAVE_MPI
192  _iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
193  _ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
194 }
195 // ----------------------------------------------------------------------------
196 // still used by geo_subdivide.cc and internally in space.cc
197 // ----------------------------------------------------------------------------
198 template <class T, class M>
199 void
200 space_base_rep<T,M>::dis_idof (const geo_element& K, std::vector<size_type>& dis_idof) const
201 {
202  freeze_guard();
203  _constit.assembly_dis_idof (get_geo(), K, dis_idof);
204 }
205 // ----------------------------------------------------------------------------
206 // name : e.g. "P1(square)", for field_expr<Expr> checks
207 // ----------------------------------------------------------------------------
208 template <class T, class M>
209 std::string
211 {
212  return _constit.name();
213 }
214 // ----------------------------------------------------------------------------
215 // u["left"]
216 // => build here the requiresd temporary indirect disarray
217 // ----------------------------------------------------------------------------
227 template <class T, class M>
230  const space_base_rep<T,M>& Wh, // Wh = space on domain
231  const std::string& dom_name // redundant: contained in Wh.get_geo()
232  ) const
233 {
234  const geo_basic<T,M>& bgd_gamma = get_geo()[dom_name];
235  return build_indirect_array (Wh, bgd_gamma);
236 }
237 template <class T, class M>
240  const space_base_rep<T,M>& Wh, // Wh = space on domain
241  const geo_basic<T,M>& bgd_gamma2
242  ) const
243 {
244  // TODO: move it to the call ?
245  const geo_basic<T,M>& bgd_gamma = bgd_gamma2.get_background_domain();
246 
247  // TODO: verifier que le domaine bgd_gamma est compatible:
248  // => il doit y avoir un meme maillage de base a bgd_gamma & Wh.get_geo
249  const space_base_rep<T,M>& Vh = *this; // Vh = space on background mesh
250  Vh.freeze_guard();
252  const geo_basic<T,M>& dom_gamma = Wh.get_geo();
253  size_type map_dim = dom_gamma.map_dimension();
254  check_macro (dom_gamma.size() == bgd_gamma.size(), "incompatible domains");
255  distributor dom_ownership = Wh.ownership();
256  distributor bgd_ownership = Vh.ownership();
257 trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_idof2bgd_idof.size=Wh.size="<<dom_ownership.size()<<"...");
258  size_type first_dom_dis_idof = dom_ownership.first_index();
259  size_type first_bgd_dis_idof = bgd_ownership.first_index();
260  std::vector<size_type> dom_dis_idofs, bgd_dis_idofs;
261  disarray<size_type, M> dom_idof2bgd_idof (dom_ownership, std::numeric_limits<size_type>::max());
262  for (size_type ige = 0, nge = dom_gamma.size(); ige < nge; ige++) {
263  const geo_element& dom_S = dom_gamma[ige];
264  const geo_element& bgd_S = bgd_gamma[ige];
265  Wh.dis_idof (dom_S, dom_dis_idofs);
266  Vh.dis_idof (bgd_S, bgd_dis_idofs);
267 trace_macro("ige="<<ige<<": dom_S="<<dom_S.name()<<dom_S.dis_ie()<<", bgd_S="<<bgd_S.name()<<bgd_S.dis_ie()
268  << ": dom_dis_idofs.size="<<dom_dis_idofs.size()
269  << ", bgd_dis_idofs.size="<<bgd_dis_idofs.size());
270  for (size_type loc_idof = 0, loc_ndof = dom_dis_idofs.size(); loc_idof < loc_ndof; loc_idof++) {
271  size_type dom_dis_idof = dom_dis_idofs [loc_idof];
272  size_type bgd_dis_idof = bgd_dis_idofs [loc_idof];
273  dom_idof2bgd_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
274 trace_macro("ige="<<ige<<": dom_idof2bgd_idof["<<dom_dis_idof<<"] = "<<bgd_dis_idof);
275  }
276  }
277  dom_idof2bgd_idof.dis_entry_assembly();
278  // move to local numbering:
279  for (size_type dom_idof = 0, dom_ndof = dom_idof2bgd_idof.size(); dom_idof < dom_ndof; dom_idof++) {
280  size_type bgd_dis_idof = dom_idof2bgd_idof [dom_idof];
281  size_type dom_dis_idof = dom_idof + first_dom_dis_idof;
282  check_macro (bgd_ownership.is_owned (bgd_dis_idof), "bgd_dis_idof="<<bgd_dis_idof<<" is out of range ["<<first_bgd_dis_idof<<":"<< bgd_ownership.last_index()<<"[");
283  size_type bgd_idof = bgd_dis_idof - first_bgd_dis_idof;
284  dom_idof2bgd_idof [dom_idof] = bgd_idof;
285  }
286 trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_idof2bgd_idof.size=Wh.size="<<dom_ownership.size()<<" done");
287  return dom_idof2bgd_idof;
288 }
289 #define _RHEOLEF_space_real(M) \
290 template <class T> \
291 space_basic<T,M> \
292 space_basic<T,M>::real() \
293 { \
294  return space_basic<T,M> (geo_basic<T,M>(details::zero_dimension()), "P1"); \
295 }
297 #ifdef _RHEOLEF_HAVE_MPI
299 #endif // _RHEOLEF_HAVE_MPI
300 #undef _RHEOLEF_space_real
301 // ----------------------------------------------------------------------------
302 // instanciation in library
303 // ----------------------------------------------------------------------------
304 #define _RHEOLEF_instanciation(T,M) \
305 template class space_base_rep<T,M>; \
306 template space_basic<T,M> space_basic<T,M>::real();
307 
309 #ifdef _RHEOLEF_HAVE_MPI
311 #endif // _RHEOLEF_HAVE_MPI
312 #undef _RHEOLEF_instanciation
313 
314 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:425
see the Float page for the full documentation
see the distributor page for the full documentation
Definition: distributor.h:62
size_type last_index(size_type iproc) const
Definition: distributor.h:157
size_type size(size_type iproc) const
Definition: distributor.h:163
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:213
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type dis_ie() const
Definition: geo_element.h:163
char name() const
Definition: geo_element.h:169
void insert(size_type dis_i)
see the integrate_option page for the full documentation
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
const Eigen::Matrix< piola< T >, Eigen::Dynamic, 1 > & get_piola(const geo_basic< T, M > &omega, const geo_element &K) const
disarray< point_basic< T >, M > _normal
Definition: space.h:259
void dis_idof(const geo_element &K, std::vector< size_type > &dis_idof) const
Definition: space.cc:200
void base_freeze_body() const
Definition: space.cc:157
std::string name() const
Definition: space.cc:210
disarray< size_type, M > build_indirect_array(const space_base_rep< T, M > &Wh, const std::string &dom_name) const
Definition: space.cc:229
space_constitution< T, M > _constit
Definition: space.h:254
const geo_basic< T, M > & get_geo() const
Definition: space.h:189
disarray< int, M > _has_nt_basis
Definition: space.h:258
const distributor & ownership() const
Definition: space.h:183
space_pair_type::size_type size_type
Definition: space.h:164
void freeze_guard() const
Definition: space.h:243
disarray< space_pair_type, M > _idof2blk_dis_iub
Definition: space.h:257
const distributor & ownership() const
#define trace_macro(message)
Definition: dis_macros.h:111
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type dis_nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
_RHEOLEF_space_real(sequential) _RHEOLEF_space_real(distributed) _RHEOLEF_instanciation(Float
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
static const bool value
Definition: distributed.h:37
point_basic< T > F
Definition: piola.h:79
Expr1::memory_type M
Definition: vec_expr_v2.h:416