egs_base_geometry.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ base geometry headers
#  Copyright (C) 2015 National Research Council Canada
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Iwan Kawrakow, 2005
#
#  Contributors:    Frederic Tessier
#                   Blake Walters
#                   Marc Chamberland
#                   Reid Townson
#                   Ernesto Mainegra-Hing
#                   Hugo Bouchard
#                   Martin Martinov
#                   Alexandre Demelo
#
###############################################################################
*/
 
 
#ifndef EGS_BASE_GEOMETRY_
#define EGS_BASE_GEOMETRY_
 
#include "egs_vector.h"
#include "egs_rndm.h"
#include "egs_input_struct.h"
 
#include <string>
#include <vector>
#include <iostream>
#include <memory>
 
using std::string;
using std::vector;
 
class EGS_Application; // forward declaration
class EGS_Input;
struct EGS_GeometryIntersections;
 
#ifdef BPROPERTY64
    typedef EGS_I64 EGS_BPType;
#elif defined BPROPERTY32
    typedef unsigned int EGS_BPType;
#elif defined BPROPERTY16
    typedef unsigned short EGS_BPType;
#else
    typedef unsigned char EGS_BPType;
#endif
 
class EGS_Label {
public:
    string      name;
    vector<int> regions;
};
 
static shared_ptr<EGS_BlockInput> geomBlockInput = make_shared<EGS_BlockInput>("geometry");
inline void setBaseGeometryInputs(bool includeMediaBlock = true) {
    geomBlockInput->addSingleInput("library", true, "The type of geometry, loaded by shared library in egs++/dso.");
    geomBlockInput->addSingleInput("name", true, "The user-declared unique name of this geometry. This is the name you may refer to elsewhere in the input file");
    geomBlockInput->addSingleInput("set label", false, "A name for the label, followed by a list of local region numbers (found by viewing only this geometry). Then use the label name elsewhere in the input file to refer to those regions. E.g. 'set label = myLabel 0 1'");
 
    if (includeMediaBlock) {
        shared_ptr<EGS_BlockInput> mediaBlock = geomBlockInput->addBlockInput("media input");
        mediaBlock->addSingleInput("media", true, "A list of media that are used in this geometry");
        mediaBlock->addSingleInput("set medium", false, "2, 3 or 4 integers defining the medium for a region or range of regions.\nFor 2: region #, medium index from the media list for this geometry (starts at 0). For 3: start region, stop region, medium index. For 4: Same as 3, plus a step size for the region range.\nNeglect this input for a homogeneous geometry of the first medium in the media list. Repeat this input to specify each medium.");
    }
}
 
class EGS_EXPORT EGS_BaseGeometry {
 
public:
 
    EGS_BaseGeometry(const string &Name);
 
    virtual ~EGS_BaseGeometry();
 
    inline bool isConvex() const {
        return is_convex;
    };
 
    virtual int inside(const EGS_Vector &x) = 0;
 
    virtual bool isInside(const EGS_Vector &x) = 0;
 
    virtual int isWhere(const EGS_Vector &x) = 0;
 
    /* getNextGeom is the equivalent of getNextParticle but for the simulation
     * object. Its goal is to determine the next state of the geometry, either
     * by synchronizing itself to the source time parameter, or by sampling its
     * own time parameter and updating itself accordingly if the source has
     * provided no time index.
     *
     * This function has a non-empty implementation in 2 cases.
     *
     * 1) it is reimplemented in any composite geometry, where it will call next
     *    geom on all of its components
     *
     * 2) it is reimplemented in the dynamic geometry class. This is where the
     *    code will find the current (non static) state of the geometry. */
    virtual void getNextGeom(EGS_RandomGenerator *rndm) {
        (void)rndm;
    };
 
    static int findRegion(EGS_Float xp, int np, const EGS_Float *p) {
        int ml = 0, mu = np;
        while (mu - ml > 1) {
            int mav = (ml+mu)/2;
            if (xp <= p[mav]) {
                mu = mav;
            }
            else {
                ml = mav;
            }
        }
        return  mu - 1;
    };
 
    virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
                       EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) = 0;
 
    virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
                                      const EGS_Vector &u);
 
    virtual EGS_Float hownear(int ireg, const EGS_Vector &x) = 0;
 
    virtual EGS_Float getVolume(int ireg) {
        (void)ireg;
        return 1.0;
    }
 
    virtual EGS_Float getBound(int idir, int ind) {
        (void)idir;
        (void)ind;
        return 0.0;
    }
 
    virtual int getNRegDir(int idir) {
        (void)idir;
        return 0;
    }
 
    int regions() const {
        return nreg;
    };
 
    virtual bool isRealRegion(int ireg) const {
        return (ireg >= 0 && ireg < nreg);
    };
 
    virtual int medium(int ireg) const {
        return region_media ? region_media[ireg] : med;
    };
 
    virtual int getMaxStep() const {
        return nreg+1;
    };
 
    virtual int computeIntersections(int ireg, int n, const EGS_Vector &x,
                                     const EGS_Vector &u, EGS_GeometryIntersections *isections);
 
    void setMedium(const string &Name);
 
    void setMedium(int start, int end, const string &Name, int delta=1);
 
    void setMedium(int imed) {
        med = imed;
    };
 
    void setMedium(int istart, int iend, int imed, int delta=1);
 
    void setMedia(EGS_Input *inp);
 
    static int nMedia();
 
    static const char *getMediumName(int ind);
 
    static int addMedium(const string &medname);
 
    static int getMediumIndex(const string &medname);
 
    virtual bool hasRhoScaling() {
        return has_rho_scaling;
    };
 
    virtual EGS_Float getRelativeRho(int ireg) const {
        return rhor && ireg >= 0 && ireg < nreg ? rhor[ireg] : 1;
    };
 
    virtual void setRelativeRho(int start, int end, EGS_Float rho);
 
    virtual void setRelativeRho(EGS_Input *);
 
    EGS_Float getMediumRho(int ind) const;
 
    virtual void setApplication(EGS_Application *app);
 
    inline bool hasBScaling() const {
        return (has_B_scaling || has_Ref_rho);
    };
 
    virtual EGS_Float getBScaling(int ireg) const {
        if (has_Ref_rho && has_B_scaling) {
            if (bfactor && ireg >= 0 && ireg < nreg) {
                return getMediumRho(medium(ireg))/rhoRef*bfactor[ireg];
            }
            else {
                return 1.0;
            }
        }
        else if (has_Ref_rho && !has_B_scaling) {
            if (ireg >= 0 && ireg < nreg) {
                return  getMediumRho(medium(ireg))/rhoRef;
            }
            else {
                return 1.0;
            }
        }
        else if (!has_Ref_rho && has_B_scaling) {
            if (bfactor && ireg >= 0 && ireg < nreg) {
                return bfactor[ireg];
 
            }
            else {
                return 1.0;
            }
        }
        else {
            return 1.0;
        }
    }
 
    virtual void setBScaling(int start, int end, EGS_Float bf);
 
    virtual void setBScaling(EGS_Input *);
 
    const string &getName() const {
        return name;
    };
 
    virtual const string &getType() const = 0;
 
    static EGS_BaseGeometry *createGeometry(EGS_Input *);
 
    static EGS_BaseGeometry *createSingleGeometry(EGS_Input *inp);
 
    static void clearGeometries();
 
    void setDebug(bool deb) {
        debug = deb;
    };
 
    static EGS_BaseGeometry *getGeometry(const string &Name);
 
    static EGS_BaseGeometry **getGeometries();
 
    static int getNGeometries();
 
    static string getUniqueName();
 
    void   setName(EGS_Input *inp);
 
    void    setBoundaryTolerance(EGS_Input *inp);
 
    void    setBoundaryTolerance(EGS_Float tol) {
        boundaryTolerance = tol;
        halfBoundaryTolerance = tol/2.;
    }
 
    virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const {
        if (!bp_array) {
            return (prop & bproperty);
        }
        return ireg >= 0 && ireg < nreg ? prop & bp_array[ireg] : false;
    };
 
    virtual void setBooleanProperty(EGS_BPType prop);
 
    virtual void addBooleanProperty(int bit);
 
    virtual void setBooleanProperty(EGS_BPType prop, int start, int end,
                                    int step=1);
 
    virtual void addBooleanProperty(int bit, int start, int end, int step=1);
 
    virtual void printInfo() const;
 
    static void describeGeometries();
 
    inline int ref() {
        return ++nref;
    };
 
 
    inline int deref() {
        return --nref;
    };
 
    static void setActiveGeometryList(int list);
 
    static int getLastError() {
        return error_flag;
    };
 
    static void resetErrorFlag() {
        error_flag = 0;
    };
 
    EGS_Float getBoundaryTolerance() {
        return boundaryTolerance;
    };
 
    virtual int getGlobalRegionOffset(const string geomName);
 
    virtual void getNumberRegions(const string &str, vector<int> &regs);
 
    virtual void getLabelRegions(const string &str, vector<int> &regs, bool sanitize=true);
 
    virtual const string &getLabelName(const int i) {
        return labels[i].name;
    }
 
    virtual int getLabelCount() {
        return labels.size();
    }
 
    int setLabels(EGS_Input *input);
 
    int setLabels(const string &inp);
 
    virtual void updatePosition(EGS_Float time) {
        (void)time;
    };
 
    bool validateRegions(const std::vector<int> &regions);
 
    virtual void finishInitialization() { };
 
    /* This method is essentially used to determine whether the simulation
     * geometry contains a dynamic geometry. Like getNextGeom(), the only
     * non-empty implementations of this function are in composite geometries
     * (where it simply calls containsDynamic on its components), and in the
     * dynamic geometry, where it will update the boolean reference to true and
     * call on its base geometry. This function was conceived to be used in the
     * view/viewcontrol (to determine whether time index objects are visible or
     * hidden), and track scoring */
    virtual void containsDynamic(bool &hasdynamic) {
        (void)hasdynamic;
    };
 
protected:
 
    int  nreg;
 
    string name;
 
    short *region_media;
 
    int med;
 
 
    bool has_rho_scaling;
 
    EGS_Float *rhor;
 
    bool has_B_scaling, has_Ref_rho;
 
    EGS_Float *bfactor;
 
    EGS_Float rhoRef;
 
    int nref;
 
    virtual void setMedia(EGS_Input *inp, int nmed, const int *med_ind);
 
    bool debug;
 
    bool is_convex;
 
    EGS_BPType   bproperty;
 
    EGS_BPType   *bp_array;
 
    EGS_Float boundaryTolerance, halfBoundaryTolerance;
 
    static int       error_flag;
 
    vector<EGS_Label> labels;
 
    EGS_Application *app;
 
private:
 
    static int active_glist;
 
};
 
struct EGS_GeometryIntersections {
    EGS_Float t;     
    EGS_Float rhof;  
    EGS_I32   ireg;  
    EGS_I32   imed;  
};
 
/* \example geometry/example1/geometry_example1.cpp
 
  Suppose that you frequently use the same complex geometry and you
  are tired of always having to include the long definition of this geometry
  into your input file. In this case it may be worthwhile to implement your
  own geometry DSO that constructs the geometry of interest with very little
  input (or even no input at all). This example illustrates how to do this
  programatically.
 
  As this is just an example, the geometry we will be constructing is not very
  complex: it consists of a cylinder inscribed in a sphere and the sphere
  inscribed in a box as shown in the figure ???. We want to be able to define
  the media of the cylinder, the sphere and the box in the input file
  (so that we can easily use different PEGS data sets, for example). All other
  dimensions of our example geometry are fixed.
 
  To create our new geometry, we make a new subdirectory in the egspp
  geometries subdirectory called \c example1. To be able to use the
  predefined rules for building geometry libraries we create a dummy
  header file \c %geometry_example1.h that does nothing else except to
  define the \c EGS_EXAMPLE1_DLL
  macro:
  \include geometry/example1/geometry_example1.h
  and put the implementation in \c %geometry_example1.cpp.
  \dontinclude geometry/example1/geometry_example1.cpp
  In this file we include the just created header file and the EGS_Input
  class header file to get access to its definition,
  needed in the \c %createGeometry() function (see below):
  \until egs_input
  We also need the declarations of the various geometry classes that we
  will use to construct our geometry:
  \until egs_envelope_geometry.h
  We then define the dimensions of the various geometrical structures
  \until box_size
 
  The remaining task is to implement a C-style \c %createGeometry() function
  exported by the DSO
  \until EGS_EXAMPLE1_EXPORT
  We check that the input object is valid and return \c null if it is not
  \until }
  We then look for a definition of the cylinder medium and print a warning if
  such a definition is missing:
  \skipline cyl_medium
  \until }
  We repeat basically the same code to obtain the sphere and the box
  media
  \until }
  \until }
  and return \c null if any of the media definition was missing:
  \until if
  Now we are ready to construct our geometry. We first create a set of
  parallel z-planes
  \until EGS_ZProjector
  As the planes will only be used by our geometry object internally, we don't
  care how the plane set is named and therefore give it a name using
  the static function EGS_BaseGeometry::getUniqueName().
  In a similar way we create a set of z-cylinders consisting of a single
  cylinder
  \until ZProjector
  and make a closed cylinder modeled as a 2D-geometry from the planes and the
  cylinder surface:
  \until getUniqueName
  We then set the medium of the cylinder
  \until setMedium
  The next step is to create the sphere as a set of spheres consisting of
  a single sphere
  \until getUniqueName
  and to fill it with the sphere medium defined in the input
  \until setMedium
  We then use an envelope geometry to inscribe the cylinder into the just
  created sphere:
  \until getUniqueName
  The next step is to create the outer box and fill it with the box medium
  defined in the input
  \until setMedium
  (as our geometry will always be centered about the origin, we pass a
   \c null affine transformation to the constructor of the box).
  We now inscribe the \c the_sphere geometry object (which is the sphere
  with the cylinder inscribed in it) into the just created box using
  again an envelope geometry,
  \until the_box
  set the name of \c the_box from the input using the inherited
  \link EGS_BaseGeometry::setName(EGS_Input*) setName() \endlink method,
  \until setName
  and return a pointer to the just created geometry
  \until }
  \until }
  That's all.
 
  In the Makefile we include the egspp make config files,
  \dontinclude geometry/example1/Makefile
  \skipline EGS_CONFIG
  \until my_machine
  set the C-preprocessor defines to the standard set of defines
  for our egspp configuration,
  \until DEFS
  define the name of our DSO,
  \until library
  and and the set of files needed to build the DSO
  \until lib_files
  Because we are using planes, cylinders, spheres, a ND-geometry, a box
  geometry and an envelope geometry to construct our geometry, we must
  link against the geometry libraries containing these geometries.
  We accomplish this by seting the \c link2_libs make variable to the
  list of needed libraries
  \until egs_genvelope
  We then include the standard rules for building a DSO
  \until include
  and set the dependencies of our object files
  \until make_depend
  We must also add the dependencies on the various geometry header files
  that we are including
  \until egs_envelope_geometry.h
  We can now build our geometry DSO by typing \c make and use
  it by including input like this in the input file
  \verbatim
  :start geometry:
      library = geometry_example1
      name = some_name
      cylinder medium = H2O521ICRU
      sphere medium = AL521ICRU
      box medium = AIR521ICRU
  :stop geometry:
  \endverbatim
  It is worth noting that we could have constructed this geometry as the
  \link EGS_UnionGeometry union \endlink of the cylinder pointed to
  by \c the_cyl, the sphere pointed to by \c sphere and the
  box pointed to by \c box. However, such an implementation would be
  much slower at run time because for each invocation of a geometry
  method the union will have to interogate the geometry methods of all
  3 objects. In contrast, the envelope geometry will be checking only
  against the box for points outside of the box, the box and the sphere for
  points inside the box but outside the sphere, the sphere and the cylinder
  for points inside the sphere but outside the cylinder and the cylinder
  only for points inside the cylinder.
 
  The complete source code of this example geometry is below
  \include geometry/example1/Makefile
  \include geometry/example1/geometry_example1.h
 
*/
 
/* \example geometry/example2/geometry_example2.cpp
 
  This example illustrates how to create the same geometry as
  described in
  <a href="geometry_2example1_2geometry__example1_8cpp-example.html">
  this example</a> but using EGS_BaseGeometry::createSingleGeometry
  to create the geometry from definitions stored in an EGS_Input objects.
  The header file is similar as in the previous
  <a href="geometry_2example1_2geometry__example1_8cpp-example.html">
  example</a>, but unlike in the previous example we don't need the
  header files of the various geometries being used.
  Our implementation consists of a single C-style function
  \c %createGeometry()
  \dontinclude geometry/example2/geometry_example2.cpp
  \skipline createGeometry
  As with the previous
  <a href="geometry_2example1_2geometry__example1_8cpp-example.html">
  example</a>, we check for valid input and extract the
  media of the cylinder, sphere and the box from the input
  \until !ok
  We will create the various geometry objects needed in our composite
  geometry by passing EGS_Input objects containing geometry definitions
  to the EGS_BaseGeometry::createSingleGeometry function.
  We start with the set of planes needed to close the cylinder.
  We first create an EGS_Input object named \c geometry that we will
  be passing to the geometry creation function:
  \until plane_input
  We then create properties for the library name,
  \until plane_library
  the plane-set type,
  \until plane_type
  the name of the set of planes,
  \until plane_name
  and the plane positions
  \until plane_positions
  and add them to the geometry definition
  \until plane_positions
  We then create the set of planes
  \until createSingleGeometry
  The definition and creation of the set of cylinders is very similar
  \until createSingleGeometry
  We now make a closed cylinder as a 2D geometry
  \until createSingleGeometry
  and set its medium
  \until setMedium
  In a similar fashion we create the sphere, use an envelope geometry
  to inscribe the cylinder into the sphere, create a box, and use again
  an envelope to inscribe the sphere into the box. The code is not shown
  here as it is lengthy and boring but can be found at the end of this page.
  Finally we set the name of the newly created geometry \c the_box and
  return it:
  \skipline the_box->setName
  \until }
  \until }
 
  The Makefile is now simpler as we don't need to link against the
  various geometry DSOs as they are loaded dynamically at run time
  by the EGS_BaseGeometry::createSingleGeometry function and
  our implementation also does not depend on the header files of the
  various geometry classes being used:
  \include geometry/example2/Makefile
 
  For completeness here is the header file:
  \include geometry/example2/geometry_example2.h
  and the complete implementation:
*/
 
 
#endif