Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_trace_move.cc
Go to the documentation of this file.
1
21// 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
37namespace rheolef {
38
39// --------------------------------------------------------------------------
40// 1) one point x
41// --------------------------------------------------------------------------
42// TODO: add dis_ie_guest
43template <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}
81template <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// --------------------------------------------------------------------------
109template <class T>
110void
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
126template <class T>
127void
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// ----------------------------------------------------------------------------
194template class geo_base_rep<Float,sequential>;
195template class geo_rep<Float,sequential>;
196#ifdef _RHEOLEF_HAVE_MPI
198template 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:497
see the distributor page for the full documentation
Definition distributor.h:69
size_type dis_size() const
global and local sizes
static const size_type decide
Definition distributor.h:83
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
base::size_type size_type
Definition geo.h:791
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.