BoundaryGeometry.h 9.26 KB
Newer Older
gsell's avatar
gsell committed
1
//
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
// 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/>.
gsell's avatar
gsell committed
17 18
//

gsell's avatar
gsell committed
19 20 21 22 23 24
/**
   \brief class BoundaryGeometry

   A GEOMETRY definition is used by most physics commands to define the
   particle charge and the reference momentum, together with some other
   data.
25

gsell's avatar
gsell committed
26 27 28 29 30
   i.e:
   G1: Geometry, FILE="input.h5"
   G2: Geometry, L=1.0, A=0.0025, B=0.0001

   :TODO: update above section
gsell's avatar
gsell committed
31 32 33 34
 */

#ifndef _OPAL_BOUNDARY_GEOMETRY_H
#define _OPAL_BOUNDARY_GEOMETRY_H
gsell's avatar
gsell committed
35

36 37 38
class OpalBeamline;
class ElementBase;

39
#include <cassert>
40 41
#include <unordered_map>
#include <unordered_set>
42
#include <array>
gsell's avatar
gsell committed
43

gsell's avatar
gsell committed
44
#include "AbstractObjects/Definition.h"
45 46
#include "Attributes/Attributes.h"
#include "Structure/SecondaryEmissionPhysics.h"
47

frey_m's avatar
frey_m committed
48 49
#include "Utilities/Util.h"

50 51
#include <gsl/gsl_rng.h>

52 53
extern Inform* gmsg;

54
class BoundaryGeometry : public Definition {
gsell's avatar
gsell committed
55

56
public:
gsell's avatar
gsell committed
57 58 59
    BoundaryGeometry();
    virtual ~BoundaryGeometry();

60 61
    virtual bool canReplaceBy (
        Object* object);
gsell's avatar
gsell committed
62

63 64
    virtual BoundaryGeometry* clone (
        const std::string& name);
gsell's avatar
gsell committed
65

66 67
    // Check the GEOMETRY data.
    virtual void execute ();
gsell's avatar
gsell committed
68

69
    // Find named GEOMETRY.
70 71
    static BoundaryGeometry* find (
        const std::string& name);
gsell's avatar
gsell committed
72

73 74
    // Update the GEOMETRY data.
    virtual void update ();
gsell's avatar
gsell committed
75

76 77
    void updateElement (
        ElementBase* element);
gsell's avatar
gsell committed
78

79
    void initialize ();
80

81
    int partInside (
gsell's avatar
gsell committed
82 83
        const Vector_t& r,
        const Vector_t& v,
84 85
        const double dt,
        Vector_t& intecoords,
86
        int& triId);
87 88 89 90

    Inform& printInfo (
        Inform& os) const;

91
    void writeGeomToVtk (std::string fn);
92

frey_m's avatar
frey_m committed
93 94
    inline std::string getFilename() const {
        return Attributes::getString(itsAttr[FGEOM]);
95
    }
gsell's avatar
gsell committed
96

frey_m's avatar
frey_m committed
97
    inline std::string getTopology() const {
98
        return Util::toUpper(Attributes::getString(itsAttr[TOPO]));
99
    }
100

101 102 103
    inline double getA() {
        return (double)Attributes::getReal(itsAttr[A]);
    }
gsell's avatar
gsell committed
104

105 106 107 108 109 110 111 112
    inline double getB() {
        return (double)Attributes::getReal(itsAttr[B]);
    }

    inline double getC() {
        return (double)Attributes::getReal(itsAttr[C]);
    }

113 114 115 116
    inline double getS() {
        return (double)Attributes::getReal(itsAttr[S]);
    }

117 118
    inline double getLength() {
        return (double)Attributes::getReal(itsAttr[LENGTH]);
119 120 121 122 123 124 125 126 127
    }

    inline double getL1() {
        return (double)Attributes::getReal(itsAttr[L1]);
    }

    inline double getL2() {
        return (double)Attributes::getReal(itsAttr[L2]);
    }
gsell's avatar
gsell committed
128 129

    /**
130 131
       Return number of boundary faces.
    */
132 133
    inline size_t getNumBFaces () {
        return Triangles_m.size();
134
    }
gsell's avatar
gsell committed
135 136

    /**
137
       Return the hr_m.
138
    */
139
    inline Vector_t gethr () {
gsell's avatar
gsell committed
140
        return voxelMesh_m.sizeOfVoxel;
141 142 143 144 145
    }
    /**
       Return the nr_m.
     */
    inline Vektor<int, 3> getnr () {
gsell's avatar
gsell committed
146
        return voxelMesh_m.nr_m;
147
    }
148

gsell's avatar
gsell committed
149
    /**
150 151 152
       Return the mincoords_m.
     */
    inline Vector_t getmincoords () {
gsell's avatar
gsell committed
153
        return minExtent_m;
154 155 156
    }
    /**
       Return the maxcoords_m.
157
    */
158
    inline Vector_t getmaxcoords () {
gsell's avatar
gsell committed
159
        return maxExtent_m;
160
    }
gsell's avatar
gsell committed
161

162 163 164 165 166 167 168 169 170 171
    inline bool getInsidePoint (Vector_t& pt) {
        if (haveInsidePoint_m == false) {
            return false;
        }
        pt = insidePoint_m;
        return true;
    }

    bool findInsidePoint (void);

adelmann's avatar
adelmann committed
172 173
    inline bool isOutsideApperture(Vector_t x) {
        if (hasApperture()) {
174
            for (size_t i = 0; i < apert_m.size(); i += 3) {
adelmann's avatar
adelmann committed
175 176 177
                if ((apert_m[i] <= x(2)) && (x(2) < apert_m[i+1])) {
                    // yes we are inside the interval
                    const double r = apert_m[i+2] * apert_m[i+2];
adelmann's avatar
adelmann committed
178
                    return ((x(0)*x(0)) + (x(1)*x(1))) > r;
adelmann's avatar
adelmann committed
179 180 181 182 183
                }
            }
        }
        return false;
    }
184

185 186 187 188 189
    int intersectRayBoundary (
        const Vector_t& P,
        const Vector_t& v,
        Vector_t& I);

190
    int fastIsInside (
191 192
        const Vector_t& reference_pt,        // [in] a reference point
        const Vector_t& P                    // [in] point to test
193 194
        );

195 196 197 198 199 200 201 202 203 204 205
    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;
206 207
    }

208 209
    inline void disableDebug(enum DebugFlags flags) {
        debugFlags_m &= ~flags;
210 211
    }

212
private:
213 214 215 216
    bool isInside (
        const Vector_t& P                    // [in] point to test
        );

217 218
    int intersectTriangleVoxel (
        const int triangle_id,
219 220 221
        const int i,
        const int j,
        const int k);
222

223
    int intersectTinyLineSegmentBoundary (
gsell's avatar
gsell committed
224 225
        const Vector_t&,
        const Vector_t&,
226 227 228
        Vector_t&,
        int&
        );
229

230
    int intersectLineSegmentBoundary (
gsell's avatar
gsell committed
231 232
        const Vector_t& P0,
        const Vector_t& P1,
233 234 235 236
        Vector_t& intersection_pt,
        int& triangle_id
        );

237
    std::string h5FileName_m;           // H5hut filename
gsell's avatar
gsell committed
238

239
    std::vector<Vector_t> Points_m;     // geometry point coordinates
240
    std::vector<std::array<unsigned int,4>> Triangles_m;   // boundary faces defined via point IDs
241
                                        // please note: 4 is correct, historical reasons!
gsell's avatar
gsell committed
242

gsell's avatar
gsell committed
243 244
    std::vector<Vector_t> TriNormals_m; // oriented normal vector of triangles
    std::vector<double> TriAreas_m;     // area of triangles
245

gsell's avatar
gsell committed
246 247
    Vector_t minExtent_m;               // minimum of geometry coordinate.
    Vector_t maxExtent_m;               // maximum of geometry coordinate.
248

249
    struct {
gsell's avatar
gsell committed
250 251 252 253 254 255 256
        Vector_t minExtent;
        Vector_t maxExtent;
        Vector_t sizeOfVoxel;
        Vektor<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

257
    } voxelMesh_m;
gsell's avatar
gsell committed
258

259
    int debugFlags_m;
260

261 262 263
    bool haveInsidePoint_m;
    Vector_t insidePoint_m;             // attribute INSIDEPOINT

264
    /*
adelmann's avatar
adelmann committed
265 266 267 268 269 270
       An additional structure to hold apperture information
       to prevent that particles go past the geometry. The user
       can specify n trippel with the form: (zmin, zmax, r)
    */
    std::vector<double> apert_m;

271 272
    gsl_rng *randGen_m;         //

gsell's avatar
gsell committed
273 274 275 276
    IpplTimings::TimerRef Tinitialize_m; // initialize geometry
    IpplTimings::TimerRef TisInside_m;
    IpplTimings::TimerRef TfastIsInside_m;
    IpplTimings::TimerRef TRayTrace_m;   // ray tracing
277
    IpplTimings::TimerRef TPartInside_m; // particle inside
gsell's avatar
gsell committed
278

279 280
    BoundaryGeometry(const BoundaryGeometry&);
    void operator= (const BoundaryGeometry&);
gsell's avatar
gsell committed
281 282

    // Clone constructor.
283
    BoundaryGeometry(const std::string& name, BoundaryGeometry* parent);
284

adelmann's avatar
adelmann committed
285 286 287
    inline bool hasApperture() {
        return (apert_m.size() != 0);
    }
288

gsell's avatar
gsell committed
289
    inline const Vector_t& getPoint (const int triangle_id, const int vertex_id) {
gsell's avatar
gsell committed
290
        assert (1 <= vertex_id && vertex_id <=3);
291
        return Points_m[Triangles_m[triangle_id][vertex_id]];
gsell's avatar
gsell committed
292 293
    }

294 295 296 297 298 299
    enum INTERSECTION_TESTS {
        SEGMENT,
        RAY,
        LINE
    };

300
    int intersectLineTriangle (
301 302 303 304 305
        const enum INTERSECTION_TESTS kind,
        const Vector_t& P0,
        const Vector_t& P1,
        const int triangle_id,
        Vector_t& I);
306

307
    inline int mapVoxelIndices2ID (const int i, const int j, const int k);
308 309
    inline Vector_t mapIndices2Voxel (const int, const int, const int);
    inline Vector_t mapPoint2Voxel (const Vector_t&);
310
    inline void computeMeshVoxelization (void);
311

312 313
    enum {
        FGEOM,    // file holding the geometry
314
        LENGTH,   // length of elliptic tube or boxcorner
315 316 317 318 319 320 321 322
        S,        // start of the geometry
        L1,       // in case of BOXCORNER first part of geometry with hight B
        L2,       // in case of BOXCORNER second part of geometry with hight B-C
        A,        // major semi-axis of elliptic tube
        B,        // minor semi-axis of ellitpic tube
        C,        // in case of BOXCORNER hight of corner
        TOPO,     // BOX, BOXCORNER, ELLIPTIC if FGEOM is selected topo is over-written
        ZSHIFT,   // Shift in z direction
323 324 325 326
        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
327
        APERTURE,    // in addition to the geometry
328
        INSIDEPOINT,
329 330
        SIZE
    };
gsell's avatar
gsell committed
331 332
};

333 334
inline Inform &operator<< (Inform& os, const BoundaryGeometry& b) {
    return b.printInfo (os);
gsell's avatar
gsell committed
335
}
336
#endif
gsell's avatar
gsell committed
337
// vi: set et ts=4 sw=4 sts=4:
338 339 340
// Local Variables:
// mode:c
// c-basic-offset: 4
341 342
// indent-tabs-mode: nil
// require-final-newline: nil
343
// End: