egs_union_geometry.cpp

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ union geometry
#  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
#                   Ernesto Mainegra-Hing
#                   Hubert Ho
#                   Marc Chamberland
#                   Hannah Gallop
#
###############################################################################
*/
 
 
#include "egs_union_geometry.h"
#include "egs_input.h"
#include "egs_functions.h"
 
using namespace std;
 
string EGS_UNIONG_LOCAL EGS_UnionGeometry::type = "EGS_UnionGeometry";
 
static bool EGS_UNIONG_LOCAL inputSet = false;
 
void EGS_UnionGeometry::setMedia(EGS_Input *,int,const int *) {
    egsWarning("EGS_UnionGeometry::setMedia: don't use this method. Use the\n"
               " setMedia() methods of the geometry objects that make up this geometry\n");
}
 
void EGS_UnionGeometry::setRelativeRho(int start, int end, EGS_Float rho) {
    setRelativeRho(0);
}
 
void EGS_UnionGeometry::setRelativeRho(EGS_Input *) {
    egsWarning("EGS_UnionGeometry::setRelativeRho(): don't use this method. "
               "Use the\n setRelativeRho() methods of the geometry objects that make "
               "up this geometry\n");
}
 
void EGS_UnionGeometry::setBScaling(int start, int end, EGS_Float bf) {
    setBScaling(0);
}
 
void EGS_UnionGeometry::setBScaling(EGS_Input *) {
    egsWarning("EGS_UnionGeometry::setBScaling(): don't use this method. "
               "Use the\n setBScaling() methods of the geometry objects that make "
               "up this geometry\n");
}
 
EGS_UnionGeometry::EGS_UnionGeometry(const vector<EGS_BaseGeometry *> &geoms,
                                     const int *priorities, const string &Name) :
    EGS_BaseGeometry(Name) {
    ng = geoms.size();
    if (ng <= 0) egsFatal("EGS_UnionGeometry::EGS_UnionGeometry: attempt "
                              " to construct a union geometry from zero geometries\n");
    is_convex = false;
    if (ng == 1) egsWarning("EGS_UnionGeometry::EGS_UnionGeometry: why "
                                "do you want to make a union out of a single geometry?\n");
    g = new EGS_BaseGeometry* [ng];
    nmax = 0;
    int j;
    int *order = new int [ng];
    if (priorities) {
        // user has definied priorities
        // order them using a very simplistic algorithm
        bool *is_used = new bool [ng];
        for (j=0; j<ng; j++) {
            is_used[j] = false;
        }
        for (j=0; j<ng; j++) {
            int imax;
            for (imax=0; imax<ng-1; imax++) if (!is_used[imax]) {
                    break;
                }
            int pmax = priorities[imax];
            for (int i=0; i<ng; i++) {
                if (!is_used[i] && priorities[i] > pmax) {
                    imax = i;
                    pmax = priorities[i];
                }
            }
            order[j] = imax;
            is_used[imax] = true;
        }
        delete [] is_used;
    }
    else {
        // user has not definied priorities
        for (j=0; j<ng; j++) {
            order[j] = j;
        }
    }
    has_rho_scaling = false;
    // now put the geometries into the array of geometries in
    // decreasing priority order.
    for (j=0; j<ng; j++) {
        int i = order[j];
        g[i] = geoms[j];
        g[i]->ref();
        int n = g[i]->regions();
        if (n > nmax) {
            nmax = n;
        }
        if (!has_rho_scaling) {
            has_rho_scaling = g[i]->hasRhoScaling();
        }
    }
    has_B_scaling = false;
    // now put the geometries into the array of geometries in
    // decreasing priority order.
    for (j=0; j<ng; j++) {
        int i = order[j];
        g[i] = geoms[j];
        g[i]->ref();
        int n = g[i]->regions();
        if (n > nmax) {
            nmax = n;
        }
        if (!has_B_scaling) {
            has_B_scaling = g[i]->hasBScaling();
        }
    }
    delete [] order;
    if (!nmax) egsFatal("EGS_UnionGeometry::EGS_UnionGeometry: all geometries"
                            " have zero regions?\n");
    nreg = nmax*ng;
}
 
EGS_UnionGeometry::~EGS_UnionGeometry() {
    for (int j=0; j<ng; j++) {
        if (!g[j]->deref()) {
            delete g[j];
        }
    }
    delete [] g;
}
 
void EGS_UnionGeometry::printInfo() const {
    EGS_BaseGeometry::printInfo();
    egsInformation(" geometries:\n");
    for (int j=0; j<ng; j++) egsInformation("   %s (type %s)\n",
                                                g[j]->getName().c_str(),g[j]->getType().c_str());
    egsInformation(
        "=======================================================\n");
}
 
extern "C" {
    static void setInputs() {
        inputSet = true;
 
        setBaseGeometryInputs(false);
 
        geomBlockInput->getSingleInput("library")->setValues({"egs_gunion"});
 
        // Format: name, isRequired, description, vector string of allowed values
        geomBlockInput->addSingleInput("geometries", true, "A list of names of previously defined geometries");
        geomBlockInput->addSingleInput("priorities", false, "A list of integers defining the geometry priorities. If neglected, the priority decreases when moving from the first to the last of the geometry list.");
    }
 
    EGS_UNIONG_EXPORT string getExample() {
        string example;
        example = {
            R"(
    # Example of egs_gunion
    #:start geometry:
        name = my_union
        library = egs_union
        geometries = my_box my_sphere
    :stop geometry:
)"};
        return example;
    }
 
    EGS_UNIONG_EXPORT shared_ptr<EGS_BlockInput> getInputs() {
        if(!inputSet) {
            setInputs();
        }
        return geomBlockInput;
    }
 
    EGS_UNIONG_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
        if (!input) {
            egsWarning("createGeometry(union): null input?\n");
            return 0;
        }
        vector<EGS_BaseGeometry *> geoms;
        vector<string> gnames;
        int err = input->getInput("geometries",gnames);
        if (err || gnames.size() < 1) {
            egsWarning("createGeometry(union): missing/wrong 'geometries' input\n");
            return 0;
        }
        for (unsigned int j=0; j<gnames.size(); j++) {
            EGS_BaseGeometry *gj = EGS_BaseGeometry::getGeometry(gnames[j]);
            if (!gj) egsWarning("createGeometry(union): no geometry named %s "
                                    "defined\n",gnames[j].c_str());
            else {
                geoms.push_back(gj);
            }
        }
        if (geoms.size() < 1) {
            egsWarning("createGeometry(union): must have at least one geometry\n");
            return 0;
        }
        vector<int> pri;
        err = input->getInput("priorities",pri);
        int *p = 0;
        if (!err) {
            if (pri.size() == geoms.size()) {
                p = new int [pri.size()];
                for (int i=0; i<pri.size(); i++) {
                    p[i] = pri[i];
                }
            }
            else egsWarning("createGeometry(union): the number of priorities (%d)"
            " is not the same as the number of geometries (%d) => ignoring\n",
            pri.size(),geoms.size());
        }
        EGS_BaseGeometry *result = new EGS_UnionGeometry(geoms,p);
        result->setName(input);
        result->setBoundaryTolerance(input);
        result->setLabels(input);
        if (p) {
            delete [] p;
        }
        return result;
    }
 
    int EGS_UnionGeometry::getGlobalRegionOffset(const string geomName) {
        // Look for the named geometry in the inscribed geometries
        for (int i=0; i<ng; i++) {
            if (g[i] && g[i]->getName() == geomName) {
                return i*nmax;
            }
        }
 
        // If it's not found above, search through the inscribed geometries in case they are composite geometries
        for (int i=0; i<ng; i++) {
            int shift = g[i]->getGlobalRegionOffset(geomName);
            if (shift >= 0) {
                shift += i*nmax;
                return shift;
            }
        }
 
        // Return -1 for not found
        return -1;
    }
 
    void EGS_UnionGeometry::getLabelRegions(const string &str, vector<int> &regs, bool sanitize) {
 
        // label defined in the sub-geometries
        vector<int> gregs;
        for (int i=0; i<ng; i++) {
 
            // add regions from set geometries
            gregs.clear();
            if (g[i]) {
                g[i]->getLabelRegions(str, gregs, sanitize);
            }
 
            // shift region numbers according to indexing style
            for (int j=0; j<gregs.size(); j++) {
                gregs[j] += i*nmax;
            }
 
            // add regions to the list
            regs.insert(regs.end(), gregs.begin(), gregs.end());
 
        }
 
        // label defined in self (union geometry input block)
        EGS_BaseGeometry::getLabelRegions(str, regs, sanitize);
 
    }
 
}