Rheolef  7.1
an efficient C++ finite element environment
geo_trace_move.cc
Go to the documentation of this file.
1 // y = trace_move(x,v);
22 // given x and v, search K in mesh such that y=x+v in K
23 // and when x+v goes outside, ray trace the y point on the boundary
24 //
25 // author: Pierre.Saramito@imag.fr
26 //
27 // date: 15 march 2012
28 //
29 // TODO: avoid the 2nd locate in trace_move ?
30 // en utilisant une table de liens S -> K (face-element: meme procs)
31 // => trace_ray renvoie un numero d'element tout de suite
32 // les liens S -> sont dans le maillage des la lecture : geo::get
33 // ils sont aussi utiles pour Galerkin-discontinu
34 //
35 #include "rheolef/geo.h"
36 
37 namespace rheolef {
38 
39 // --------------------------------------------------------------------------
40 // 1) one point x
41 // --------------------------------------------------------------------------
42 // TODO: add dis_ie_guest
43 template <class T, class M>
46  const point_basic<T>& x,
47  const point_basic<T>& v,
48  point_basic<T>& y) const
49 {
50  y = x+v;
51  size_type dis_ie = _locator.seq_locate (*this,y);
52  if (dis_ie != std::numeric_limits<size_type>::max()) {
53  // y = x+v falls inside Omega
54  return dis_ie;
55  }
56  // y = x+v falls outside Omega: compute y = ray(x,v) intersection with the Omega boundary
57  bool hit = _tracer_ray_boundary.seq_trace_ray_boundary (*this, x, v, y);
58  if (hit) {
59  // TODO: hit knows the face boundary: could get the element without locator
60  dis_ie = _locator.seq_locate (*this,y);
61  check_macro (dis_ie != std::numeric_limits<size_type>::max(), "invalid ray computation");
62  return dis_ie;
63  }
64  // x is on the boundary and also y (tangential velocity)
65  dis_ie = _locator.seq_locate (*this,y);
66  if (dis_ie != std::numeric_limits<size_type>::max()) {
67  return dis_ie;
68  }
69  // x was outside Omega ? most of the time, y is exactly on the boundary:
70  // boundary projection in the multi-connected case with: y = omega.nearest(x+v)
71  point_basic<T> y_nearest;
72  dis_ie = seq_nearest (y, y_nearest);
73  if (dis_ie != std::numeric_limits<size_type>::max()) {
74  y = y_nearest;
75  return dis_ie;
76  }
77  // machine precision problem ?
78  error_macro ("invalid seq_trace_move");
79  return std::numeric_limits<size_type>::max();
80 }
81 template <class T, class M>
84  const point_basic<T>& x,
85  const point_basic<T>& v,
86  point_basic<T>& y) const
87 {
88  y = x+v;
89  size_type dis_ie = _locator.dis_locate (*this,y);
90  if (dis_ie != std::numeric_limits<size_type>::max()) {
91  // y = x+v falls inside Omega
92  return dis_ie;
93  }
94  // y = x+v falls outside Omega: compute y = ray(x,v) intersection with the Omega boundary
95  bool hit = _tracer_ray_boundary.dis_trace_ray_boundary (*this, x, v, y);
96  if (hit) {
97  // TODO: hit knows the face boundary: could get the element without locator
98  dis_ie = _locator.dis_locate (*this,y);
99  check_macro (dis_ie != std::numeric_limits<size_type>::max(), "invalid ray computation");
100  return dis_ie;
101  }
102  // x was outside Omega ?
103  error_macro ("invalid dis_trace_move");
104  return std::numeric_limits<size_type>::max();
105 }
106 // --------------------------------------------------------------------------
107 // 2) disarrays x & v : groups comms in the distributed case
108 // --------------------------------------------------------------------------
109 template <class T>
110 void
116 {
117  check_macro (x.ownership() == v.ownership(), "invalid ranges for x and v argumenst");
118  dis_ie.resize (x.ownership());
119  y.resize (x.ownership());
120  for (size_type i = 0, n = x.size(); i < n; i++) {
121  // TODO: use dis_ie[i] as a guest
122  dis_ie[i] = base::seq_trace_move (x[i], v[i], y[i]);
123  }
124 }
125 #ifdef _RHEOLEF_HAVE_MPI
126 template <class T>
127 void
131  disarray<size_type, distributed>& dis_ie, // in: guest; out: localized
133 {
134  trace_macro("trace_move...");
135  check_macro (x.ownership() == v.ownership(), "invalid ranges for x and v argumenst");
136  if (dis_ie.ownership() != x.ownership()) dis_ie.resize (x.ownership());
137  y.resize (x.ownership());
138  // ---------------------------------------------------------------
139  // 1) locate: solve all x(i)+y(i) that remains inside the domain:
140  // ---------------------------------------------------------------
141  for (size_type i = 0, n = x.size(); i < n; i++) {
142  y[i] = x[i] + v[i];
143  }
144  bool do_stop_when_failed = false;
145  trace_macro("trace_move/locate...");
146  locate (y, dis_ie, do_stop_when_failed);
147  // -----------------------------------------------------------------
148  // 2) for all x+v that goes outside
149  // compact all failed points and run trace_ray_boundary (fld_x,fld_v)
150  // ---------------------------------------------------------------
151  std::list<size_type> failed;
152  for (size_type i = 0, n = x.size(); i < n; i++) {
153  if (dis_ie[i] != std::numeric_limits<size_type>::max()) continue;
154  // here x[i)+v[i] cross the boundary:
155  failed.push_back (i);
156  }
157  distributor fld_ownership (distributor::decide, base::comm(), failed.size());
158  if (fld_ownership.dis_size() == 0) {
159  // all are solved: nothing more to do
160  trace_macro("trace_move done(1)");
161  return;
162  }
163  disarray<point_basic<T> > fld_x (fld_ownership);
164  disarray<point_basic<T> > fld_v (fld_ownership);
165  disarray<point_basic<T> > fld_y (fld_ownership);
166  disarray<size_type> fld_dis_ie (fld_ownership, std::numeric_limits<size_type>::max());
167  typename std::list<size_type>::const_iterator iter = failed.begin();
168  for (size_type fld_i = 0, fld_n = fld_x.size(); fld_i < fld_n; ++fld_i, ++iter) {
169  size_type i = *iter;
170  fld_x [fld_i] = x[i];
171  fld_v [fld_i] = v[i];
172  }
173  do_stop_when_failed = true;
174  trace_macro("trace_move/ray_boundary...");
175  trace_ray_boundary (fld_x, fld_v, fld_dis_ie, fld_y, do_stop_when_failed);
176  trace_macro("trace_move/ray_boundary = "
177  << fld_x.dis_size() << "/" << x.dis_size() << " = "
178  << 1.0*fld_x.dis_size()/x.dis_size());
179  // -----------------------------------------------------------------
180  // 3) unpack the result
181  // ---------------------------------------------------------------
182  iter = failed.begin();
183  for (size_type fld_i = 0, fld_n = fld_x.size(); fld_i < fld_n; ++fld_i, ++iter) {
184  size_type i = *iter;
185  dis_ie[i] = fld_dis_ie [fld_i];
186  y[i] = fld_y [fld_i];
187  }
188  trace_macro("trace_move done(2)");
189 }
190 #endif // _RHEOLEF_HAVE_MPI
191 // ----------------------------------------------------------------------------
192 // instanciation in library
193 // ----------------------------------------------------------------------------
194 template class geo_base_rep<Float,sequential>;
195 template class geo_rep<Float,sequential>;
196 #ifdef _RHEOLEF_HAVE_MPI
197 template class geo_base_rep<Float,distributed>;
198 template class geo_rep<Float,distributed>;
199 #endif // _RHEOLEF_HAVE_MPI
200 
201 } // namespace rheolef
see the disarray page for the full documentation
Definition: disarray.h:459
see the distributor page for the full documentation
Definition: distributor.h:62
size_type dis_size() const
global and local sizes
Definition: distributor.h:207
static const size_type decide
Definition: distributor.h:76
geo_element_hack::size_type size_type
Definition: geo.h:260
base class for M=sequential or distributed meshes representations
Definition: geo.h:528
base::size_type size_type
Definition: geo.h:533
size_type seq_trace_move(const point_basic< T > &x, const point_basic< T > &v, point_basic< T > &y) const
size_type dis_trace_move(const point_basic< T > &x, const point_basic< T > &v, point_basic< T > &y) const
base::size_type size_type
Definition: geo.h:934
sequential mesh representation
Definition: geo.h:778
#define trace_macro(message)
Definition: dis_macros.h:111
#define error_macro(message)
Definition: dis_macros.h:49
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.