← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/lib/hazel_1d.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/lib/hazel_1d.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
#include "rheolef/hazel_1d.h"
 
22
#include "rheolef/georep.h"
 
23
using namespace std;
 
24
void
 
25
hazel_1d::insert (pair<size_type,size_type>& v2e_ix, size_type K_idx)
 
26
{
 
27
  static const size_type none = numeric_limits<size_type>::max();
 
28
  if (v2e_ix.first == none) { v2e_ix.first = K_idx; return; }
 
29
  if (v2e_ix.second == none) { v2e_ix.second = K_idx; return; }
 
30
  error_macro ("unexpected 1d mesh connectivity arround element " << K_idx);
 
31
}
 
32
hazel_1d::hazel_1d (const georep& Ih)
 
33
 : vector<vector<size_type> >(),
 
34
   _interval()
 
35
{
 
36
  initialize(Ih);
 
37
}
 
38
void
 
39
hazel_1d::initialize (const georep& Ih)
 
40
{
 
41
  check_macro (Ih.map_dimension() == 1, "only 1D mesh here");
 
42
  // ------------------------------------
 
43
  // 1) initialize vertex2edge table
 
44
  // ------------------------------------
 
45
  size_type dimension = Ih.dimension(); // number of physical coordinates
 
46
  static const size_type none = numeric_limits<size_type>::max();
 
47
  pair<size_type,size_type> none_pair = make_pair(none,none);
 
48
  vector<pair<size_type,size_type> > vertex2edge (Ih.n_vertex(), none_pair);
 
49
  georep::const_iterator K = Ih.begin();
 
50
  georep::const_iterator_vertex x = Ih.begin_vertex();
 
51
  for (size_type k = 0; k < Ih.size(); k++) {
 
52
    insert (vertex2edge[K[k][0]], k);
 
53
    insert (vertex2edge[K[k][1]], k);
 
54
  }
 
55
  // sanity check
 
56
  for (size_type i = 0; i < vertex2edge.size(); i++) {
 
57
    check_macro (vertex2edge[i].first != none, "vertex " << i
 
58
        << " has no elements");
 
59
  }
 
60
  // -------------------------------------
 
61
  // 2) build list of connected components
 
62
  // -------------------------------------
 
63
  list<list<size_type> > path_list;
 
64
  vector<bool> mark (Ih.n_vertex(), false);
 
65
  list<pair<Float,Float> > interval_list;
 
66
  for (size_type i = 0; i < Ih.n_vertex(); i++) {
 
67
    // search for an extremity of a connected component
 
68
    // i.e a boundary none that is included in only one element
 
69
    if (mark[i]) {
 
70
      continue;
 
71
    }
 
72
    if (dimension == 1 && vertex2edge[i].second != none) {
 
73
      continue;
 
74
    }
 
75
    path_list.push_back(list<size_type>());
 
76
    list<size_type>& path = *(path_list.rbegin());
 
77
    mark[i] = true;
 
78
    size_type k = vertex2edge[i].first;
 
79
    size_type j = i;
 
80
    do {
 
81
      size_type j_prec = j;
 
82
      j = (K[k][0] != j_prec) ? K[k][0] : K[k][1];
 
83
      path.push_back(k);
 
84
      if (mark[j]) {
 
85
        // closed surface
 
86
        break;
 
87
      }
 
88
      mark[j] = true;
 
89
      size_type k_prev = k;
 
90
      k = (vertex2edge[j].first != k_prev) ? vertex2edge[j].first : vertex2edge[j].second;
 
91
    } while (k != none);
 
92
    // path from i to j : ordered by increasing coordinate ?
 
93
    if (dimension == 1) {
 
94
      pair<Float,Float> boundary;
 
95
      if (x[i][0] < x[j][0]) {
 
96
        boundary = make_pair(x[i][0], x[j][0]);
 
97
      } else {
 
98
        boundary = make_pair(x[j][0], x[i][0]);
 
99
        path.reverse();
 
100
      }
 
101
      interval_list.push_back(boundary);
 
102
    }
 
103
  }
 
104
  // copy list<list> into hard-sized vector<vector> : save memory
 
105
  size_type n_component = path_list.size();
 
106
  vector<vector<size_type> >::resize(n_component);
 
107
  size_type i_component = 0;
 
108
  for (list<list<size_type> >::const_iterator p = path_list.begin(); p != path_list.end(); p++, i_component++) {
 
109
    const list<size_type>& path = *p;
 
110
    vector<size_type>& new_path = vector<vector<size_type> >::operator[](i_component);
 
111
    new_path.resize(path.size());
 
112
    copy (path.begin(), path.end(), new_path.begin());
 
113
  }
 
114
  // copy also interval list into hard-sized version (when dimension = 1)
 
115
  if (dimension == 1) {
 
116
    _interval.resize(n_component);
 
117
    copy (interval_list.begin(), interval_list.end(), _interval.begin());
 
118
  }
 
119
}
 
120
hazel_1d::size_type
 
121
hazel_1d::localize (const georep& Ih, const Float& x0) const
 
122
{
 
123
  check_macro (_interval.size() != 0, "cannot localize in a " << Ih.dimension() 
 
124
        << "D line path");
 
125
  // step 1 : find connex component containing x0
 
126
  static const size_type none = numeric_limits<size_type>::max();
 
127
  size_type i0 = none;
 
128
  for (size_type i = 0; i < _interval.size(); i++) {
 
129
    if (_interval[i].first <= x0 && _interval[i].second >= x0) {
 
130
      i0 = i;
 
131
      break;
 
132
    }
 
133
  }
 
134
  if (i0 == none) return none;
 
135
  const vector<size_type>& component = operator[](i0);
 
136
  size_type n = component.size();
 
137
  // step 2 : find element containing x0 in component
 
138
  // by using dichotomy on elements
 
139
  georep::const_iterator K = Ih.begin();
 
140
  georep::const_iterator_vertex x = Ih.begin_vertex();
 
141
  size_type kmin = 0;
 
142
  size_type kmax = n;
 
143
  size_type k = 0;
 
144
  while (kmin+1 < kmax) {
 
145
    k = size_type((kmin + kmax)/2.0);
 
146
    Float xk = min(x[K[k][0]][0],x[K[k][1]][0]);
 
147
    if (x0 >= xk) kmin = k;
 
148
    else          kmax = k;
 
149
  }
 
150
  k = kmin;
 
151
  // K[k] contains x
 
152
#undef PARANO
 
153
#ifdef PARANO
 
154
  Float x1 = min(x[K[k][0]][0],x[K[k][1]][0]);
 
155
  Float x2 = max(x[K[k][0]][0],x[K[k][1]][0]);
 
156
  check_macro (x1 <= x0 && x0 <= x2, "locate error");
 
157
#endif // PARANO
 
158
  return k;
 
159
}
 
160
bool
 
161
hazel_1d::localize_nearest (const georep& Ih, const Float& x0,
 
162
        Float& y0, size_type& k) const
 
163
{
 
164
  static const size_type none = numeric_limits<size_type>::max();
 
165
  k = localize (Ih, x0);
 
166
  if (k != none) { y0 = x0; return true; }
 
167
  // loop on intervals and search closest boundary
 
168
  size_type i0 = none;
 
169
  Float dist_min = numeric_limits<Float>::max();
 
170
  size_t extremity = 0;
 
171
  for (size_type i = 0; i < _interval.size(); i++) {
 
172
    Float d1 = fabs(x0 - _interval[i].first);
 
173
    if (d1 < dist_min) {
 
174
      y0 = _interval[i].first;
 
175
      i0 = i;
 
176
      extremity = 0;
 
177
      dist_min = d1;
 
178
    }
 
179
    Float d2 = fabs(x0 - _interval[i].second);
 
180
    if (d2 < dist_min) {
 
181
      y0 = _interval[i].second;
 
182
      i0 = i;
 
183
      extremity = 1;
 
184
      dist_min = d2;
 
185
    }
 
186
  }
 
187
  const vector<size_type>& component = operator[](i0);
 
188
  if (extremity == 0) {
 
189
    k = component[0];
 
190
  } else {
 
191
    k = component[component.size()-1];
 
192
  }
 
193
  return false;
 
194
}
 
195
// trace x0+v :
 
196
// return true if x0+v is inside the mesh
 
197
// otherwise return nearest element containing x,
 
198
// also t such that x=x0+t*v (for the method of characteristic)
 
199
bool
 
200
hazel_1d::trace (const georep& Ih, const Float& x0, const Float& v, Float& x, Float& t, size_type& k) const 
 
201
{
 
202
  static const size_type none = numeric_limits<size_type>::max();
 
203
  bool is_inside = localize_nearest (Ih, x0+v, x, k);
 
204
  if (is_inside) {
 
205
    t = 1;
 
206
  } else {
 
207
    t = (1+v != 1) ? (x-x0)/v : 0;
 
208
  }
 
209
  return is_inside;
 
210
}