← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/geo.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/geo.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 sequential 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
#include "rheolef/geo.h"
 
22
#include "rheolef/geo_domain.h"
 
23
#include "rheolef/rheostream.h"
 
24
#include "rheolef/iorheo.h"
 
25
 
 
26
#include "geo_header.h"
 
27
 
 
28
namespace rheolef {
 
29
 
 
30
// =========================================================================
 
31
// geo_base_rep members
 
32
// =========================================================================
 
33
template <class T, class M>
 
34
geo_base_rep<T,M>::~geo_base_rep()
 
35
{
 
36
#ifdef TO_CLEAN
 
37
  trace_macro ("dstor geo_base_rep: name="<<name()<<", dim="<<dimension()<<", map_dim="<<map_dimension()<<" size="<<size());
 
38
#endif // TO_CLEAN
 
39
}
 
40
template <class T, class M>
 
41
std::string
 
42
geo_base_rep<T,M>::name() const
 
43
{
 
44
  if (_serial_number == 0) return _name;
 
45
  return _name + "-" + itos(_serial_number);
 
46
}
 
47
template <class T, class M>
 
48
bool
 
49
geo_base_rep<T,M>::have_domain_indirect (const std::string& name) const
 
50
{
 
51
  for (typename std::vector<domain_indirect_basic<M> >::const_iterator
 
52
        iter = _domains.begin(), last = _domains.end(); iter != last; iter++) {
 
53
    const domain_indirect_basic<M>& dom = *iter;
 
54
    if (name == dom.name()) return true;
 
55
  }
 
56
  return false;
 
57
}
 
58
template <class T, class M>
 
59
const domain_indirect_basic<M>&
 
60
geo_base_rep<T,M>::get_domain_indirect (const std::string& name) const
 
61
{
 
62
  for (typename std::vector<domain_indirect_basic<M> >::const_iterator
 
63
        iter = _domains.begin(), last = _domains.end(); iter != last; iter++) {
 
64
    const domain_indirect_basic<M>& dom = *iter;
 
65
    if (name == dom.name()) return dom;
 
66
  }
 
67
  error_macro ("undefined domain \""<<name<<"\" in mesh \"" << _name << "\"");
 
68
  return *(_domains.begin()); // not reached
 
69
}
 
70
template <class T, class M>
 
71
void
 
72
geo_base_rep<T,M>::insert_domain_indirect (const domain_indirect_basic<M>& dom) const
 
73
{
 
74
  for (typename std::vector<domain_indirect_basic<M> >::iterator
 
75
        iter = _domains.begin(), last = _domains.end(); iter != last; iter++) {
 
76
    domain_indirect_basic<M>& curr_dom = *iter;
 
77
    if (dom.name() == curr_dom.name()) {
 
78
      // a domain with the same name already exists: replace it
 
79
      curr_dom = dom;
 
80
      return;
 
81
    }
 
82
  }
 
83
  // insert a new domain:
 
84
  _domains.push_back (dom);
 
85
}
 
86
template <class T, class M>
 
87
typename geo_base_rep<T,M>::size_type
 
88
geo_base_rep<T,M>::size (size_type dim) const
 
89
{
 
90
  size_type sz = 0;
 
91
  for (size_type variant = reference_element::first_variant_by_dimension(dim);
 
92
                 variant < reference_element:: last_variant_by_dimension(dim); variant++) {
 
93
    sz += _geo_element [variant].size();
 
94
  }
 
95
  return sz;
 
96
}
 
97
template <class T, class M>
 
98
typename geo_base_rep<T,M>::size_type
 
99
geo_base_rep<T,M>::dis_size (size_type dim) const
 
100
{
 
101
  size_type sz = 0;
 
102
  for (size_type variant = reference_element::first_variant_by_dimension(dim);
 
103
                 variant < reference_element:: last_variant_by_dimension(dim); variant++) {
 
104
    sz += _geo_element [variant].dis_size();
 
105
  }
 
106
  return sz;
 
107
}
 
108
template <class T, class M>
 
109
void
 
110
geo_base_rep<T,M>::compute_bbox()
 
111
{
 
112
  // first, compute bbox sequentialy:
 
113
  for (size_type j = 0; j < _dimension; j++) {
 
114
    _xmin[j] =  std::numeric_limits<T>::max();
 
115
    _xmax[j] = -std::numeric_limits<T>::max();
 
116
  }
 
117
  for (size_type j = _dimension+1; j < 3; j++) {
 
118
    _xmin [j] = _xmax [j] = 0;
 
119
  }
 
120
  for (size_type inod = 0, nnod = _node.size(); inod < nnod; inod++) {
 
121
    const point_basic<T>& x = _node [inod];
 
122
    for (size_type j = 0 ; j < _dimension; j++) {
 
123
      _xmin[j] = min(x[j], _xmin[j]);
 
124
      _xmax[j] = max(x[j], _xmax[j]);
 
125
    }
 
126
  }
 
127
  // distributed case: min & max are grouped accross procs 
 
128
#ifdef _RHEOLEF_HAVE_MPI
 
129
  if (is_distributed<M>::value) {
 
130
    for (size_type j = 0 ; j < _dimension; j++) {
 
131
      _xmin[j] = mpi::all_reduce (comm(), _xmin[j], mpi::minimum<T>());
 
132
      _xmax[j] = mpi::all_reduce (comm(), _xmax[j], mpi::maximum<T>());
 
133
    }
 
134
  }
 
135
#endif // _RHEOLEF_HAVE_MPI
 
136
}
 
137
template <class T, class M>
 
138
typename geo_base_rep<T,M>::const_reference
 
139
geo_base_rep<T,M>::get_geo_element (size_type dim, size_type ige) const
 
140
{
 
141
  size_type first_ige = 0, last_ige = 0;
 
142
  for (size_type variant = reference_element::first_variant_by_dimension(dim);
 
143
                 variant < reference_element:: last_variant_by_dimension(dim); variant++) {
 
144
    last_ige += _geo_element [variant].size();
 
145
    if (ige < last_ige) return _geo_element [variant] [ige - first_ige];
 
146
    first_ige = last_ige;
 
147
  }
 
148
  error_macro ("geo_element index " << ige << " out of range [0:"<< last_ige << "[");
 
149
  return _geo_element [0][0]; // not reached
 
150
}
 
151
template <class T, class M>
 
152
typename geo_base_rep<T,M>::reference
 
153
geo_base_rep<T,M>::get_geo_element (size_type dim, size_type ige)
 
154
{
 
155
  size_type first_ige = 0, last_ige = 0;
 
156
  for (size_type variant = reference_element::first_variant_by_dimension(dim);
 
157
                 variant < reference_element:: last_variant_by_dimension(dim); variant++) {
 
158
    last_ige += _geo_element [variant].size();
 
159
    if (ige < last_ige) return _geo_element [variant] [ige - first_ige];
 
160
    first_ige = last_ige;
 
161
  }
 
162
  error_macro ("geo_element index " << ige << " out of range [0:"<< last_ige << "[");
 
163
  return _geo_element [0][0]; // not reached
 
164
}
 
165
template <class T, class M>
 
166
typename geo_base_rep<T,M>::node_type
 
167
geo_base_rep<T,M>::piola (const geo_element& K, const node_type& hat_x) const
 
168
{
 
169
  // TODO: use basis_on_lattice [K.variant] instead of values[]
 
170
  std::vector<T> values;
 
171
  reference_element hat_K (K.variant());
 
172
  _numbering.get_basis().eval (hat_K, hat_x, values);
 
173
 
 
174
  std::vector<size_type> dis_inod;
 
175
  _numbering.dis_idof (_gs, K, dis_inod);
 
176
  node_type x;
 
177
  for (size_type loc_inod = 0, loc_nnod = values.size(); loc_inod < loc_nnod; loc_inod++) {
 
178
    const node_type& xnod = dis_node(dis_inod[loc_inod]);
 
179
    const T& coeff = values[loc_inod];
 
180
    for (size_type i = 0, d = dimension(); i < d; i++) {
 
181
      x[i] += coeff*xnod[i];
 
182
    }
 
183
  }
 
184
  return x;
 
185
}
 
186
// ----------------------------------------------------------------------------
 
187
// set order & resize node array ; internal node are not computed
 
188
// since they depend on the curved boundary information (cad data)
 
189
// that is not aivailable here
 
190
// ----------------------------------------------------------------------------
 
191
#define _RHEOLEF_reset_order(M)                                         \
 
192
template <class T>                                                      \
 
193
void                                                                    \
 
194
geo_rep<T,M>::reset_order (size_type new_order)                         \
 
195
{                                                                       \
 
196
  if (new_order == base::_numbering.degree()) return;                   \
 
197
  base::_numbering.set_degree (new_order);                              \
 
198
  size_type dis_nnod = base::_numbering.dis_ndof (base::_gs, base::_map_dimension); \
 
199
  size_type     nnod = base::_numbering.ndof     (base::_gs, base::_map_dimension); \
 
200
  base::_gs.node_ownership = distributor (nnod, base::comm(), nnod);    \
 
201
  array<point_basic<T>, M> new_node (base::_gs.node_ownership);         \
 
202
  for (size_type iv = 0, nv = base::_gs.ownership_by_dimension[0].size(); iv < nv; iv++) { \
 
203
    new_node [iv] = base::_node [iv];                                   \
 
204
  }                                                                     \
 
205
  base::_node = new_node;                                               \
 
206
  build_external_entities ();                                           \
 
207
}
 
208
_RHEOLEF_reset_order(sequential)
 
209
#ifdef _RHEOLEF_HAVE_MPI
 
210
_RHEOLEF_reset_order(distributed)
 
211
#endif // _RHEOLEF_HAVE_MPI
 
212
#undef _RHEOLEF_reset_order
 
213
// ----------------------------------------------------------------------------
 
214
// instanciation in library
 
215
// ----------------------------------------------------------------------------
 
216
template class geo_base_rep<Float,sequential>;
 
217
template class geo_rep<Float,sequential>;
 
218
#ifdef _RHEOLEF_HAVE_MPI
 
219
template class geo_base_rep<Float,distributed>;
 
220
template class geo_rep<Float,distributed>;
 
221
#endif // _RHEOLEF_HAVE_MPI
 
222
 
 
223
template geo_basic<Float,sequential> compact (const geo_basic<Float,sequential>&);
 
224
 
 
225
} // namespace rheolef