EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_union_geometry.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ union geometry headers
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Iwan Kawrakow, 2005
25 #
26 # Contributors: Ernesto Mainegra-Hing
27 # Frederic Tessier
28 # Reid Townson
29 #
30 ###############################################################################
31 */
32 
33 
39 #ifndef EGS_UNION_GEOMETRY_
40 #define EGS_UNION_GEOMETRY_
41 
42 #ifdef WIN32
43 
44  #ifdef BUILD_UNIONG_DLL
45  #define EGS_UNIONG_EXPORT __declspec(dllexport)
46  #else
47  #define EGS_UNIONG_EXPORT __declspec(dllimport)
48  #endif
49  #define EGS_UNIONG_LOCAL
50 
51 #else
52 
53  #ifdef HAVE_VISIBILITY
54  #define EGS_UNIONG_EXPORT __attribute__ ((visibility ("default")))
55  #define EGS_UNIONG_LOCAL __attribute__ ((visibility ("hidden")))
56  #else
57  #define EGS_UNIONG_EXPORT
58  #define EGS_UNIONG_LOCAL
59  #endif
60 
61 #endif
62 
63 
64 #include "egs_base_geometry.h"
65 
66 #include<vector>
67 using std::vector;
68 
133 class EGS_UNIONG_EXPORT EGS_UnionGeometry : public EGS_BaseGeometry {
134 
135 public:
136 
140  EGS_UnionGeometry(const vector<EGS_BaseGeometry *> &geoms,
141  const int *priorities = 0, const string &Name = "");
142 
144 
145  bool isRealRegion(int ireg) const {
146  if (ireg < 0 || ireg >= nreg) {
147  return false;
148  }
149  int j = ireg/nmax;
150  return g[j]->isRealRegion(ireg-j*nmax);
151  };
152 
153  bool isInside(const EGS_Vector &x) {
154  for (int j=0; j<ng; j++) if (g[j]->isInside(x)) {
155  return true;
156  }
157  return false;
158  };
159 
160  int isWhere(const EGS_Vector &x) {
161  for (int j=0; j<ng; j++) {
162  int ij = g[j]->isWhere(x);
163  if (ij >= 0) {
164  return ij + j*nmax;
165  }
166  }
167  return -1;
168  };
169 
170  int inside(const EGS_Vector &x) {
171  return isWhere(x);
172  };
173 
174  int medium(int ireg) const {
175  if (ireg < 0 || ireg >= nreg) {
176  return -1;
177  }
178  int j = ireg/nmax;
179  return g[j]->medium(ireg-j*nmax);
180  };
181 
182  int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
183  EGS_Float &t,int *newmed=0, EGS_Vector *normal=0) {
184  if (ireg >= 0) {
185  // we are inside, set current geometry
186  int jg = ireg/nmax;
187  // distance to boundary and new region in this geometry.
188  int inew = g[jg]->howfar(ireg-jg*nmax,x,u,t,newmed,normal);
189  // jgnew is the geometry after the step
190  int jgnew = inew >= 0 ? jg : -1;
191  // geometries are now ordered in decreasing priorities.
192  // Then,
193  // - if the particle remains in the current geometry,
194  // we only need to check geometries up to jg-1
195  // For these geometries the particle must be outside,
196  // otherwise it would have been in one of them
197  // - if the particle exits the current geometry, then
198  // we must also check jg+1...ng-1
199  for (int j=0; j<jg; j++) {
200  int ii = g[j]->howfar(-1,x,u,t,newmed,normal);
201  if (ii >= 0) {
202  jgnew = j;
203  inew = ii;
204  }
205  }
206  if (inew < 0) {
207  // the particle didn't enter any of the higher priority
208  // geometries but exits the current one.
209  // => we need to check if the particle is in one
210  // of the lower priority geometries at the exit point.
211  EGS_Vector xnew(x+u*t);
212  for (int j=jg+1; j<ng; j++) {
213  int ii = g[j]->isWhere(xnew);
214  if (ii >= 0) {
215  // when exiting jg, particle is in region ii of geometry j
216  // we don't need to check other geometries because they
217  // have a lower priority.
218  jgnew = j;
219  inew = ii;
220  break;
221  }
222  }
223  }
224  if (inew < 0) {
225  return inew;
226  }
227  if (newmed) {
228  *newmed = g[jgnew]->medium(inew);
229  }
230  return inew + jgnew*nmax;
231  }
232  // if here, we are currently outside of all geometries in the union.
233  int jg, inew=-1;
234  for (int j=0; j<ng; j++) {
235  int ii = g[j]->howfar(-1,x,u,t,newmed,normal);
236  if (ii >= 0) {
237  jg = j;
238  inew = ii;
239  }
240  }
241  return inew < 0 ? -1 : inew + jg*nmax;
242  };
243 
244  EGS_Float hownear(int ireg, const EGS_Vector &x) {
245  if (ireg >= 0) {
246  int jg = ireg/nmax;
247  EGS_Float tmin = g[jg]->hownear(ireg-jg*nmax,x);
248  if (tmin <= 0) {
249  return 0;
250  }
251  // tmin is now the perpendicular distance to a boundary
252  // in the current geometry.
253  // we need to check all geometries with higher priorities
254  // i.e., all geometries between 0 and jg-1.
255  // as their priorities are higher, we know that we are
256  // outside of such geometries.
257  for (int j=jg-1; j>=0; --j) {
258  EGS_Float t = g[j]->hownear(-1,x);
259  if (t < tmin) {
260  tmin = t;
261  if (tmin <= 0) {
262  return 0;
263  }
264  }
265  }
266  return tmin;
267  }
268  // if here, we are outside of all geomtries in the union.
269  EGS_Float tmin = veryFar;
270  for (int j=ng-1; j>=0; j--) {
271  EGS_Float t = g[j]->hownear(-1,x);
272  if (t < tmin) {
273  tmin = t;
274  }
275  if (tmin <= 0) {
276  return 0;
277  }
278  }
279  return tmin;
280  };
281 
282  int getMaxStep() const {
283  int nstep = 1;
284  for (int j=0; j<ng; ++j) {
285  nstep += g[j]->getMaxStep();
286  }
287  return nstep;
288  };
289 
290  const string &getType() const {
291  return type;
292  };
293 
294  void printInfo() const;
295 
296  EGS_Float getRelativeRho(int ireg) const {
297  if (ireg < 0 || ireg >= nreg) {
298  return 1;
299  }
300  int jg = ireg/nmax;
301  return g[jg]->getRelativeRho(ireg-jg*nmax);
302  };
303  void setRelativeRho(int start, int end, EGS_Float rho);
304  void setRelativeRho(EGS_Input *);
305 
306  void setBScaling(int start, int end, EGS_Float rho);
307  void setBScaling(EGS_Input *);
308  EGS_Float getBScaling(int ireg) const {
309  if (ireg < 0 || ireg >= nreg) {
310  return 1;
311  }
312  int jg = ireg/nmax;
313  return g[jg]->getBScaling(ireg-jg*nmax);
314  };
315 
316  virtual void getLabelRegions(const string &str, vector<int> &regs);
317 
318 protected:
319 
321  int ng;
322  int nmax;
323  static string type;
324 
331  void setMedia(EGS_Input *,int,const int *);
332 
333 };
334 
335 #endif
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
A class representing 3D vectors.
Definition: egs_vector.h:56
EGS_BaseGeometry ** g
the geometries that form the union.
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0)=0
Calculate the distance to a boundary from x along the direction u.
const EGS_Float veryFar
A very large float.
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
virtual EGS_Float getBScaling(int ireg) const
Get the B field scaling factor in region ireg.
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
A geometry constructed as the union of other geometries.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
static string type
the geometry type
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
int nmax
max. number of regions in all of the geoms.
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
EGS_BaseGeometry class header file.
int ng
number of geometries.