EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_envelope_geometry.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ envelope geometry
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: Frederic Tessier
27 # Ernesto Mainegra-Hing
28 #
29 ###############################################################################
30 */
31 
32 
38 #include "egs_envelope_geometry.h"
39 #include "egs_input.h"
40 #include "egs_functions.h"
41 
42 #include <cstdlib>
43 
44 using namespace std;
45 
46 string EGS_ENVELOPEG_LOCAL EGS_EnvelopeGeometry::type = "EGS_EnvelopeGeometry";
47 string EGS_ENVELOPEG_LOCAL EGS_FastEnvelope::type = "EGS_FastEnvelope";
48 
49 void EGS_EnvelopeGeometry::setMedia(EGS_Input *,int,const int *) {
50  egsWarning("EGS_EnvelopeGeometry::setMedia: don't use this method. Use the\n"
51  " setMedia() methods of the geometry objects that make up this geometry\n");
52 }
53 
54 void EGS_EnvelopeGeometry::setRelativeRho(int start, int end, EGS_Float rho) {
55  setRelativeRho(0);
56 }
57 
58 void EGS_EnvelopeGeometry::setRelativeRho(EGS_Input *) {
59  egsWarning("EGS_EnvelopeGeometry::setRelativeRho(): don't use this method."
60  " Use the\n setRelativeRho methods of the geometry objects that make up"
61  " this geometry\n");
62 }
63 
64 void EGS_EnvelopeGeometry::setBScaling(int start, int end, EGS_Float bf) {
65  setBScaling(0);
66 }
67 
68 void EGS_EnvelopeGeometry::setBScaling(EGS_Input *) {
69  egsWarning("EGS_EnvelopeGeometry::setBScaling(): don't use this method."
70  " Use the\n setBScaling methods of the geometry objects that make up"
71  " this geometry\n");
72 }
73 
74 void EGS_FastEnvelope::setMedia(EGS_Input *,int,const int *) {
75  egsWarning("EGS_FastEnvelope::setMedia: don't use this method. Use the\n"
76  " setMedia() methods of the geometry objects that make up this geometry\n");
77 }
78 
79 void EGS_FastEnvelope::setRelativeRho(int start, int end, EGS_Float rho) {
80  setRelativeRho(0);
81 }
82 
83 void EGS_FastEnvelope::setRelativeRho(EGS_Input *) {
84  egsWarning("EGS_FastEnvelope::setRelativeRho(): don't use this method."
85  " Use the\n setRelativeRho methods of the geometry objects that make up"
86  " this geometry\n");
87 }
88 
89 void EGS_FastEnvelope::setBScaling(int start, int end, EGS_Float bf) {
90  setBScaling(0);
91 }
92 
93 void EGS_FastEnvelope::setBScaling(EGS_Input *) {
94  egsWarning("EGS_FastEnvelope::setBScaling(): don't use this method."
95  " Use the\n setBScaling methods of the geometry objects that make up"
96  " this geometry\n");
97 }
98 
99 struct EGS_ENVELOPEG_LOCAL EnvelopeAux {
100  EGS_BaseGeometry *g;
101  int nreg;
102  int *regs;
103  EnvelopeAux() : nreg(0) {};
104  ~EnvelopeAux() {
105  if (nreg>0) {
106  delete [] regs;
107  }
108  };
109 };
110 
111 EGS_EnvelopeGeometry::EGS_EnvelopeGeometry(EGS_BaseGeometry *G,
112  const vector<EGS_BaseGeometry *> &geoms, const string &Name,
113  bool newindexing) :
114  EGS_BaseGeometry(Name), reg_to_inscr(0), local_start(0) {
115  if (!G) {
116  egsFatal("EGS_EnvelopeGeometry: base geometry must not be null\n");
117  }
118  g = G;
119  g->ref();
120  new_indexing = false;
121  nbase = g->regions();
122  n_in = geoms.size();
123  nmax = 1;
124  int nreg_inscribed = 0;
125  if (n_in > 0) {
126  geometries = new EGS_BaseGeometry * [n_in];
127  for (int j=0; j<n_in; j++) {
128  geometries[j] = geoms[j];
129  geometries[j]->ref();
130  int nj = geometries[j]->regions();
131  if (nj > nmax) {
132  nmax = nj;
133  }
134  nreg_inscribed += nj;
135  }
136  }
137  nreg = nbase + n_in*nmax;
138  is_convex = g->isConvex();
139 
140  has_rho_scaling = g->hasRhoScaling();
141  if (!has_rho_scaling) {
142  for (int j=0; j<n_in; j++) {
143  if (geometries[j]->hasRhoScaling()) {
144  has_rho_scaling = true;
145  break;
146  }
147  }
148  }
149  has_B_scaling = g->hasBScaling();
150  if (!has_B_scaling) {
151  for (int j=0; j<n_in; j++) {
152  if (geometries[j]->hasBScaling()) {
153  has_B_scaling = true;
154  break;
155  }
156  }
157  }
158 
159  if (!n_in) {
160  return;
161  }
162  if (newindexing) {
163  new_indexing = true;
164  local_start = new int [n_in];
165  reg_to_inscr = new int [nreg_inscribed];
166  int ir=0;
167  for (int j=0; j<n_in; j++) {
168  int n = geometries[j]->regions();
169  local_start[j] = nbase + ir;
170  for (int i=0; i<n; i++) {
171  reg_to_inscr[ir++] = j;
172  }
173  }
174  nreg = nbase + ir;
175  }
176 }
177 
178 EGS_FastEnvelope::EGS_FastEnvelope(EGS_BaseGeometry *G,
179  const vector<EnvelopeAux *> &fgeoms, const string &Name,
180  int newindexing) :
181  EGS_BaseGeometry(Name), reg_to_inscr(0), local_start(0) {
182  if (!G) {
183  egsFatal("EGS_FastEnvelope: base geometry must not be null\n");
184  }
185  g = G;
186  g->ref();
187  n_in = fgeoms.size();
188  nmax = 1;
189  new_indexing = false;
190  if (!n_in) {
191  egsFatal("EGS_FastEnvelope: no inscribed geometries!\n");
192  }
193  geometries = new EGS_BaseGeometry * [n_in];
194  nbase = g->regions();
195  int *iaux = new int [nbase];
196  int j;
197  for (j=0; j<nbase; j++) {
198  iaux[j] = 0;
199  }
200  int nlist = 0;
201  int nreg_inscribed = 0;
202  for (j=0; j<fgeoms.size(); j++) {
203  geometries[j] = fgeoms[j]->g;
204  geometries[j]->ref();
205  int nj = geometries[j]->regions();
206  nreg_inscribed += nj;
207  if (nj > nmax) {
208  nmax = nj;
209  }
210  for (int i=0; i<fgeoms[j]->nreg; i++) {
211  int k = fgeoms[j]->regs[i];
212  if (k >= 0 && k < nbase) {
213  ++iaux[k];
214  ++nlist;
215  }
216  }
217  }
218  n_start = new int [nbase+1];
219  glist = new int [nlist];
220  int ilist=0;
221  for (j=0; j<nbase; j++) {
222  n_start[j] = ilist;
223  if (iaux[j] > 0) {
224  for (int l=0; l<fgeoms.size(); l++) {
225  for (int i=0; i<fgeoms[l]->nreg; i++) {
226  int k = fgeoms[l]->regs[i];
227  if (k == j) {
228  glist[ilist++] = l;
229  }
230  }
231  }
232  }
233  }
234  n_start[nbase] = ilist;
235  nreg = nbase + n_in*nmax;
236  is_convex = g->isConvex();
237  has_rho_scaling = g->hasRhoScaling();
238  if (!has_rho_scaling) {
239  for (int j=0; j<n_in; j++) {
240  if (geometries[j]->hasRhoScaling()) {
241  has_rho_scaling = true;
242  break;
243  }
244  }
245  }
246 
247  has_B_scaling = g->hasBScaling();
248  if (!has_B_scaling) {
249  for (int j=0; j<n_in; j++) {
250  if (geometries[j]->hasBScaling()) {
251  has_B_scaling = true;
252  break;
253  }
254  }
255  }
256 
257  delete [] iaux;
258  if (newindexing) {
259  new_indexing = true;
260  local_start = new int [n_in];
261  reg_to_inscr = new int [nreg_inscribed];
262  int ir=0;
263  for (int j=0; j<n_in; j++) {
264  int n = geometries[j]->regions();
265  local_start[j] = nbase + ir;
266  for (int i=0; i<n; i++) {
267  reg_to_inscr[ir++] = j;
268  }
269  }
270  nreg = nbase + ir;
271  }
272 }
273 
274 EGS_EnvelopeGeometry::~EGS_EnvelopeGeometry() {
275  if (!g->deref()) {
276  delete g;
277  }
278  for (int j=0; j<n_in; j++) {
279  if (!geometries[j]->deref()) {
280  delete geometries[j];
281  }
282  }
283  delete [] geometries;
284  if (new_indexing) {
285  if (local_start) {
286  delete [] local_start;
287  }
288  if (reg_to_inscr) {
289  delete [] reg_to_inscr;
290  }
291  }
292 }
293 
294 EGS_FastEnvelope::~EGS_FastEnvelope() {
295  if (!g->deref()) {
296  delete g;
297  }
298  for (int j=0; j<n_in; j++) {
299  if (!geometries[j]->deref()) {
300  delete geometries[j];
301  }
302  }
303  delete [] geometries;
304  delete [] n_start;
305  delete [] glist;
306  if (new_indexing) {
307  if (local_start) {
308  delete [] local_start;
309  }
310  if (reg_to_inscr) {
311  delete [] reg_to_inscr;
312  }
313  }
314 }
315 
316 void EGS_EnvelopeGeometry::printInfo() const {
318  egsInformation(" base geometry = %s (type %s)\n",g->getName().c_str(),
319  g->getType().c_str());
320  egsInformation(" inscribed geometries:\n");
321  for (int j=0; j<n_in; j++) egsInformation(" %s (type %s)\n",
322  geometries[j]->getName().c_str(),geometries[j]->getType().c_str());
324  "=======================================================\n");
325 }
326 
327 void EGS_FastEnvelope::printInfo() const {
329  egsInformation(" base geometry = %s (type %s)\n",g->getName().c_str(),
330  g->getType().c_str());
331  egsInformation(" inscribed geometries:\n");
332  for (int j=0; j<n_in; j++) egsInformation(" %s (type %s)\n",
333  geometries[j]->getName().c_str(),geometries[j]->getType().c_str());
335  "=======================================================\n");
336 }
337 
338 
339 static char EGS_ENVELOPEG_LOCAL eeg_message1[] =
340  "createGeometry(envelope geometry): %s\n";
341 static char EGS_ENVELOPEG_LOCAL eeg_message2[] =
342  "null input?";
343 static char EGS_ENVELOPEG_LOCAL eeg_message3[] =
344  "no 'base geometry' input?";
345 static char EGS_ENVELOPEG_LOCAL eeg_message4[] =
346  "incorrect base geometry definition";
347 static char EGS_ENVELOPEG_LOCAL eeg_message5[] =
348  "missing/incorrect 'base geometry' input";
349 static char EGS_ENVELOPEG_LOCAL eeg_message6[] =
350  "createGeometry(envelope geometry): no geometry with name %s defined\n";
351 static char EGS_ENVELOPEG_LOCAL eeg_message7[] =
352  "no inscirebed geometries defined?. I hope you know what you are doing";
353 static char EGS_ENVELOPEG_LOCAL eeg_message8[] =
354  "an error occured while constructing inscibed geometries";
355 
356 static char EGS_ENVELOPEG_LOCAL eeg_keyword1[] = "base geometry";
357 static char EGS_ENVELOPEG_LOCAL eeg_keyword2[] = "geometry";
358 static char EGS_ENVELOPEG_LOCAL eeg_keyword3[] = "inscribed geometries";
359 
360 extern "C" {
361 
362  EGS_ENVELOPEG_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
363  if (!input) {
364  egsWarning(eeg_message1,eeg_message2);
365  return 0;
366  }
367  //
368  // *** new indexing style
369  //
370  int indexing = 0;
371  input->getInput("new indexing style",indexing);
372  //
373  // *** Base geometry
374  //
375  EGS_Input *i = input->takeInputItem(eeg_keyword1);
376  if (!i) {
377  egsWarning(eeg_message1,eeg_message3);
378  return 0;
379  }
380  EGS_Input *ig = i->takeInputItem(eeg_keyword2);
381  EGS_BaseGeometry *g;
382  if (ig) { // defined inline
384  delete ig;
385  if (!g) {
386  egsWarning(eeg_message1,eeg_message4);
387  delete i;
388  return 0;
389  }
390  }
391  else { // defined via a name of a previously defined geometry
392  string bgname;
393  int err = i->getInput(eeg_keyword1,bgname);
394  delete i;
395  if (err) {
396  egsWarning(eeg_message1,eeg_message5);
397  return 0;
398  }
399  g = EGS_BaseGeometry::getGeometry(bgname);
400  if (!g) {
401  egsWarning(eeg_message6,bgname.c_str());
402  return 0;
403  }
404  }
405  vector<EnvelopeAux *> fgeoms;
406  EGS_Input *ix;
407  while ((ix = input->takeInputItem("inscribe in regions")) != 0) {
408  vector<string> values;
409  ix->getInput("inscribe in regions",values);
410  if (values.size() < 2) egsWarning("createGeometry(envelope geometry):"
411  " %d inputs for 'inscribe in regions'? 2 or more are needed\n",values.size());
412  else {
414  if (!gj) {
415  egsWarning(eeg_message6,values[0].c_str());
416  }
417  else {
418  EnvelopeAux *aux = new EnvelopeAux;
419  aux->g = gj;
420  aux->nreg = values.size()-1;
421  aux->regs = new int [aux->nreg];
422  for (int j=0; j<aux->nreg; j++) {
423  aux->regs[j] = atoi(values[j+1].c_str());
424  }
425  fgeoms.push_back(aux);
426  }
427  }
428  delete ix;
429  }
430  if (fgeoms.size() > 0) {
431  EGS_BaseGeometry *result = new EGS_FastEnvelope(g,fgeoms,"",indexing);
432  result->setName(input);
433  result->setBoundaryTolerance(input);
434  for (int j=0; j<fgeoms.size(); j++) {
435  delete fgeoms[j];
436  }
437  return result;
438  }
439 
440  //
441  // *** Inscribed geometries
442  //
443  i = input->takeInputItem(eeg_keyword3);
444  vector<EGS_BaseGeometry *> geoms;
445  if (!i) {
446  if (fgeoms.size()) {
447  egsWarning(eeg_message1,eeg_message7);
448  }
449  }
450  else {
451  // first try for inscribed geometries defined inline
452  while ((ig = i->takeInputItem(eeg_keyword2)) != 0) {
454  delete ig;
455  if (!gj) {
456  egsWarning(eeg_message1,eeg_message8);
457  }
458  else {
459  geoms.push_back(gj);
460  }
461  }
462  // if geoms.size() is 0, check if inscribed geometries are defined
463  // via names of previously defined geometries.
464  if (!geoms.size()) {
465  vector<string> igeoms;
466  int err = i->getInput(eeg_keyword3,igeoms);
467  if (err || !igeoms.size()) {
468  egsWarning(eeg_message1,eeg_message7);
469  }
470  else {
471  for (unsigned int j=0; j<igeoms.size(); j++) {
472  EGS_BaseGeometry *gj =
474  if (!gj) {
475  egsWarning(eeg_message6,igeoms[j].c_str());
476  }
477  else {
478  geoms.push_back(gj);
479  }
480  }
481  }
482  }
483  delete i;
484  }
485  /*
486  EGS_BaseGeometry *result = fgeoms.size() ?
487  new EGS_EnvelopeGeometry(g,fgeoms,geoms) :
488  new EGS_EnvelopeGeometry(g,geoms);
489  */
490  EGS_BaseGeometry *result = new EGS_EnvelopeGeometry(g,geoms,"",indexing);
491  result->setName(input);
492  result->setLabels(input);
493  return result;
494 
495  }
496 
497 
498  void EGS_EnvelopeGeometry::getLabelRegions(const string &str, vector<int> &regs) {
499 
500  // label defined in the envelope geometry
501  g->getLabelRegions(str, regs);
502 
503  // label defined in the inscribed geometries
504  vector<int> gregs;
505  int shift=0;
506  for (int i=0; i<n_in; i++) {
507 
508  // add regions from set geometries
509  gregs.clear();
510  if (geometries[i]) {
511  geometries[i]->getLabelRegions(str, gregs);
512  }
513 
514  // shift region numbers according to indexing style
515  if (new_indexing) {
516  shift = local_start[i];
517  }
518  else {
519  shift = nbase+i*nmax;
520  }
521  for (int j=0; j<gregs.size(); j++) {
522  gregs[j] += shift;
523  }
524 
525  // add regions to the list
526  regs.insert(regs.end(), gregs.begin(), gregs.end());
527 
528  }
529 
530  // label defined in self (envelope geometry input block)
532 
533  }
534 
535  void EGS_FastEnvelope::getLabelRegions(const string &str, vector<int> &regs) {
536 
537  // label defined in the envelope geometry
538  g->getLabelRegions(str, regs);
539 
540  // label defined in the inscribed geometries
541  vector<int> gregs;
542  int shift=0;
543  for (int i=0; i<n_in; i++) {
544 
545  // add regions from set geometries
546  gregs.clear();
547  if (geometries[i]) {
548  geometries[i]->getLabelRegions(str, gregs);
549  }
550 
551  // shift region numbers according to indexing style
552  if (new_indexing) {
553  shift = local_start[i];
554  }
555  else {
556  shift = nmax;
557  }
558  for (int j=0; j<gregs.size(); j++) {
559  gregs[j] += shift;
560  }
561 
562  // add regions to the list
563  regs.insert(regs.end(), gregs.begin(), gregs.end());
564 
565  }
566 
567  // label defined in self (envelope geometry input block)
569 
570  }
571 
572 }
EGS_BaseGeometry * g
The envelope geometry.
virtual const string & getType() const =0
Get the geometry type.
EGS_BaseGeometry ** geometries
The inscribed geometries.
int * local_start
First region for each inscribed geometry.
void setMedia(EGS_Input *, int, const int *)
Don&#39;t set media for an envelope geometry.
int regions() const
Returns the number of local regions in this geometry.
int deref()
Decrease the reference count to this geometry.
virtual void printInfo() const
Print information about this geometry.
An envelope geometry class.
EGS_Input class header file.
static string type
Geometry type.
int setLabels(EGS_Input *input)
Set the labels from an input block.
bool hasRhoScaling() const
Does this geometry object have a mass density scaling feature?
Global egspp functions header file.
EGS_BaseGeometry * g
The envelope geometry.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
int n_in
Number of inscribed geometries.
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
void setMedia(EGS_Input *, int, const int *)
Don&#39;t set media for an envelope geometry.
bool isConvex() const
Is the geometry convex?
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
int * reg_to_inscr
Region to inscribed geometry conversion.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
bool new_indexing
If true, use new indexing style.
An envelope geometry: header.
bool hasBScaling() const
Does this geometry object have a B field scaling feature?
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
int ref()
Increase the reference count to this geometry.
static EGS_BaseGeometry * createSingleGeometry(EGS_Input *inp)
Create a single geometry from the input inp.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
const string & getName() const
Get the name of this geometry.
An envelope geometry class.
bool new_indexing
If true, use new indexing style.
int n_in
Number of inscribed geometries.
int * local_start
First region for each inscribed geometry.
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
int * reg_to_inscr
Region to inscribed geometry conversion.
EGS_BaseGeometry ** geometries
The inscribed geometries.
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
static string type
Geometry type.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.