BoundaryGeometry.h

Go to the documentation of this file.

//
// Implementation of the BoundaryGeometry class
//
// Copyright (c) 200x - 2020, Achim Gsell,
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL.  If not, see <https://www.gnu.org/licenses/>.
//
 
#ifndef _OPAL_BOUNDARY_GEOMETRY_H
#define _OPAL_BOUNDARY_GEOMETRY_H
 
class OpalBeamline;
class ElementBase;
 
#include "AbstractObjects/Definition.h"
#include "Attributes/Attributes.h"
#include "Utilities/Util.h"
#include "Utility/IpplTimings.h"
#include "Utility/PAssert.h"
 
#include "Utilities/Random.h"
 
#include <array>
#include <unordered_map>
#include <unordered_set>
#include <vector>
 
enum class Topology : unsigned short { RECTANGULAR, BOXCORNER, ELLIPTIC };
 
class BoundaryGeometry : public Definition {
public:
    BoundaryGeometry();
    virtual ~BoundaryGeometry();
 
    virtual bool canReplaceBy(Object* object);
 
    virtual BoundaryGeometry* clone(const std::string& name);
 
    // Check the GEOMETRY data.
    virtual void execute();
 
    // Find named GEOMETRY.
    static BoundaryGeometry* find(const std::string& name);
 
    // Update the GEOMETRY data.
    virtual void update();
 
    void updateElement(ElementBase* element);
 
    void initialize();
 
    int partInside(
            const Vector_t<double, 3>& r, const Vector_t<double, 3>& v, const double dt,
            Vector_t<double, 3>& intecoords, int& triId);
 
    Inform& printInfo(Inform& os) const;
 
    void writeGeomToVtk(std::string fn);
 
    inline std::string getFilename() const { return Attributes::getString(itsAttr[FGEOM]); }
 
    inline Topology getTopology() const {
        static const std::unordered_map<std::string, Topology> stringTopology_s = {
                {"RECTANGULAR", Topology::RECTANGULAR},
                {"BOXCORNER", Topology::BOXCORNER},
                {"ELLIPTIC", Topology::ELLIPTIC}};
        Topology topo = stringTopology_s.at(Attributes::getString(itsAttr[TOPO]));
        return topo;
    }
 
    inline double getA() { return (double)Attributes::getReal(itsAttr[A]); }
 
    inline double getB() { return (double)Attributes::getReal(itsAttr[B]); }
 
    inline double getC() { return (double)Attributes::getReal(itsAttr[C]); }
 
    inline double getS() { return (double)Attributes::getReal(itsAttr[S]); }
 
    inline double getLength() { return (double)Attributes::getReal(itsAttr[LENGTH]); }
 
    inline double getL1() { return (double)Attributes::getReal(itsAttr[L1]); }
 
    inline double getL2() { return (double)Attributes::getReal(itsAttr[L2]); }
 
    inline size_t getNumBFaces() { return Triangles_m.size(); }
 
    inline Vector_t<double, 3> gethr() { return voxelMesh_m.sizeOfVoxel; }
    inline Vector_t<int, 3> getnr() { return voxelMesh_m.nr_m; }
 
    inline Vector_t<double, 3> getmincoords() { return minExtent_m; }
    inline Vector_t<double, 3> getmaxcoords() { return maxExtent_m; }
 
    inline bool getInsidePoint(Vector_t<double, 3>& pt) {
        if (haveInsidePoint_m == false) {
            return false;
        }
        pt = insidePoint_m;
        return true;
    }
 
    bool findInsidePoint(void);
 
    int intersectRayBoundary(
            const Vector_t<double, 3>& P, const Vector_t<double, 3>& v, Vector_t<double, 3>& I);
 
    int fastIsInside(
            const Vector_t<double, 3>& reference_pt,  // [in] a reference point
            const Vector_t<double, 3>& P              // [in] point to test
    );
 
    enum DebugFlags {
        debug_isInside                         = 0x0001,
        debug_fastIsInside                     = 0x0002,
        debug_intersectRayBoundary             = 0x0004,
        debug_intersectLineSegmentBoundary     = 0x0008,
        debug_intersectTinyLineSegmentBoundary = 0x0010,
        debug_PartInside                       = 0x0020,
    };
 
    inline void enableDebug(enum DebugFlags flags) { debugFlags_m |= flags; }
 
    inline void disableDebug(enum DebugFlags flags) { debugFlags_m &= ~flags; }
 
private:
    bool isInside(
            const Vector_t<double, 3>& P  // [in] point to test
    );
 
    int intersectTriangleVoxel(const int triangle_id, const int i, const int j, const int k);
 
    int intersectTinyLineSegmentBoundary(
            const Vector_t<double, 3>&, const Vector_t<double, 3>&, Vector_t<double, 3>&, int&);
 
    int intersectLineSegmentBoundary(
            const Vector_t<double, 3>& P0, const Vector_t<double, 3>& P1,
            Vector_t<double, 3>& intersection_pt, int& triangle_id);
 
    std::string h5FileName_m;  // H5hut filename
 
    std::vector<Vector_t<double, 3>> Points_m;  // geometry point coordinates
    std::vector<std::array<unsigned int, 4>>
            Triangles_m;  // boundary faces defined via point IDs
                          // please note: 4 is correct, historical reasons!
 
    std::vector<Vector_t<double, 3>> TriNormals_m;  // oriented normal vector of triangles
    std::vector<double> TriAreas_m;                 // area of triangles
 
    Vector_t<double, 3> minExtent_m;  // minimum of geometry coordinate.
    Vector_t<double, 3> maxExtent_m;  // maximum of geometry coordinate.
 
    struct {
        Vector_t<double, 3> minExtent;
        Vector_t<double, 3> maxExtent;
        Vector_t<double, 3> sizeOfVoxel;
        Vector_t<int, 3> nr_m;  // number of intervals of geometry in X,Y,Z direction
        std::unordered_map<
                int,  // map voxel IDs ->
                std::unordered_set<int>>
                ids;  // intersecting triangles
 
    } voxelMesh_m;
 
    int debugFlags_m;
 
    bool haveInsidePoint_m;
    Vector_t<double, 3> insidePoint_m;  // attribute INSIDEPOINT
 
    gsl_rng* randGen_m;  //
 
    IpplTimings::TimerRef Tinitialize_m;  // initialize geometry
    IpplTimings::TimerRef TisInside_m;
    IpplTimings::TimerRef TfastIsInside_m;
    IpplTimings::TimerRef TRayTrace_m;    // ray tracing
    IpplTimings::TimerRef TPartInside_m;  // particle inside
 
    BoundaryGeometry(const BoundaryGeometry&);
    void operator=(const BoundaryGeometry&);
 
    // Clone constructor.
    BoundaryGeometry(const std::string& name, BoundaryGeometry* parent);
 
    inline const Vector_t<double, 3>& getPoint(const int triangle_id, const int vertex_id) {
        PAssert(1 <= vertex_id && vertex_id <= 3);
        return Points_m[Triangles_m[triangle_id][vertex_id]];
    }
 
    enum INTERSECTION_TESTS { SEGMENT, RAY, LINE };
 
    int intersectLineTriangle(
            const enum INTERSECTION_TESTS kind, const Vector_t<double, 3>& P0,
            const Vector_t<double, 3>& P1, const int triangle_id, Vector_t<double, 3>& I);
 
    inline int mapVoxelIndices2ID(const int i, const int j, const int k);
    inline Vector_t<double, 3> mapIndices2Voxel(const int, const int, const int);
    inline Vector_t<double, 3> mapPoint2Voxel(const Vector_t<double, 3>&);
    inline void computeMeshVoxelization(void);
 
    enum {
        FGEOM,     // file holding the geometry
        LENGTH,    // length of elliptic tube or boxcorner
        S,         // start of the geometry
        L1,        // in case of BOXCORNER first part of geometry with height B
        L2,        // in case of BOXCORNER second part of geometry with height B-C
        A,         // major semi-axis of elliptic tube
        B,         // minor semi-axis of ellitpic tube
        C,         // in case of BOXCORNER height of corner
        TOPO,      // RECTANGULAR, BOXCORNER, ELLIPTIC if FGEOM is selected TOPO is over-written
        ZSHIFT,    // Shift in z direction
        XYZSCALE,  // Multiplicative scaling factor for coordinates
        XSCALE,    // Multiplicative scaling factor for x-coordinates
        YSCALE,    // Multiplicative scaling factor for y-coordinates
        ZSCALE,    // Multiplicative scaling factor for z-coordinates
        INSIDEPOINT,
        SIZE
    };
};
 
inline Inform& operator<<(Inform& os, const BoundaryGeometry& b) { return b.printInfo(os); }
#endif