← Back to branch summary

~sophie-middleton08/maus/devel

« back to all changes in this revision

Viewing changes to src/legacy/Interface/Mesh.hh

  • Committer: Chris Rogers
  • Date: 2012-10-03 07:19:33 UTC
  • mfrom: (659.1.40 release-candidate)
  • Revision ID: chris.rogers@stfc.ac.uk-20121003071933-kgrhvl1ec6w2jmug
Tags: MAUS-v0.3, MAUS-v0.3.3
MAUS-v0.3.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/legacy/Interface/Mesh.hh
5
5
#include <vector>
6
6
#include <algorithm>
7
7
#include <math.h>
8
 
 
9
 
#include "Maths/Vector.hh"
10
 
 
11
 
//
12
 
//Mesh is a set of classes to describe a grid of points in various dimensions etc
13
 
//Initially this is for interpolation type algorithms, but one might also use these
14
 
//things for e.g. numerical integration etc
15
 
//
16
 
//So I define Mesh which is an abstract data type (i.e. can't be imlpemented, only inherited from)
17
 
//which has an Iterator object that behaves like a normal stl iterator
18
 
//
19
 
//Concrete classes are TwoDGrid,  ThreeDGrid and NDGrid, rectangular grids with either constant 
20
 
//spacing or dynamically allocated spacing
21
 
//
22
 
//See also TriangularMesh for an example of a (Delaunay) 2D triangular meshing algorithm
23
 
//
24
 
 
25
 
class VectorMap;
26
 
 
27
 
class Mesh
28
 
{
29
 
public:
30
 
    //Iterable sub-class
31
 
    class Iterator;
32
 
    //Constructor
33
 
    Mesh();
34
 
    //Return the first point on the mesh
35
 
    virtual Mesh::Iterator Begin() const = 0;
36
 
    //Return the last point+1 on the mesh
37
 
    virtual Mesh::Iterator End()   const = 0;
38
 
    //Return a copy of child object - copy constructor will only copy the parent object
39
 
    virtual Mesh* Clone() = 0;
40
 
    //Return the "Dual" of the mesh
41
 
    //Dual is a polyhedron that has centre of each face as a point on the mesh
42
 
    virtual Mesh* Dual () = 0;
43
 
    virtual ~Mesh();
44
 
 
45
 
    //Overload these functions for Mesh::Iterator to work with your base class
46
 
    //Return the position of a point in the mesh
47
 
    virtual void Position(const Mesh::Iterator& it, double * position) const      = 0;
48
 
    //Return the dimension of the mesh
49
 
    virtual int  PositionDimension() const = 0;
50
 
    //Map from iterator to an integer used to index the iterator
51
 
    virtual int  ToInteger(const Mesh::Iterator& lhs)      const = 0;
52
 
    //Find the point on the mesh nearest to some point
53
 
    virtual Mesh::Iterator Nearest(const double* position) const = 0;
54
 
 
55
 
protected:
56
 
    //Change the position of an iterator
57
 
    //This is the minimal set of functions you need to define for the iterator to work
58
 
    virtual Mesh::Iterator& AddEquals(Mesh::Iterator& lhs, int difference) const = 0;
59
 
    virtual Mesh::Iterator& SubEquals(Mesh::Iterator& lhs, int difference) const = 0;
60
 
    virtual Mesh::Iterator& AddEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const = 0;
61
 
    virtual Mesh::Iterator& SubEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const = 0;
62
 
    virtual Mesh::Iterator& AddOne   (Mesh::Iterator& lhs) const = 0;
63
 
    virtual Mesh::Iterator& SubOne   (Mesh::Iterator& lhs) const = 0;
64
 
    //Check relative position
65
 
    virtual bool IsGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const = 0;
66
 
 
67
 
public:
68
 
    //Iterator needs to know about internal workings of the Mesh to work properly
69
 
    friend Mesh::Iterator  operator++(Mesh::Iterator& lhs, int);
70
 
    friend Mesh::Iterator  operator--(Mesh::Iterator& lhs, int);
71
 
    friend Mesh::Iterator& operator++(Mesh::Iterator& lhs);
72
 
    friend Mesh::Iterator& operator--(Mesh::Iterator& lhs);
73
 
    friend Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
74
 
    friend Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
75
 
    friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
76
 
    friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
77
 
    friend Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const int& rhs);
78
 
    friend Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const int& rhs);
79
 
    friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs,       const int& rhs);
80
 
    friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs,       const int& rhs);
81
 
 
82
 
    friend bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
83
 
    friend bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
84
 
    friend bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
85
 
    friend bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
86
 
    friend bool operator< (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
87
 
    friend bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
88
 
 
89
 
};
90
 
 
91
 
 
92
 
//Iterator object
93
 
//Used to loop over some, or all, points in the mesh, as in stl
94
 
//Enables e,g, generic file I/O operations on a field map without knowing
95
 
//the details of what is in the mesh or what shape it has...
96
 
class Mesh::Iterator
97
 
{
98
 
public:
99
 
    Iterator();
100
 
    Iterator(const Mesh::Iterator& in);
101
 
    Iterator(std::vector<int> state, const Mesh* mesh);
102
 
    virtual ~Iterator();
103
 
    const Mesh::Iterator&  operator=(const Mesh::Iterator& rhs);
104
 
    //return position referenced by the iterator
105
 
    virtual void                Position(double* point) const;
106
 
    virtual std::vector<double> Position()              const;
107
 
 
108
 
    int                         ToInteger()             const {return _mesh->ToInteger(*this);}
109
 
    std::vector<int>            State()                 const {return _state;}
110
 
    int&       operator[](int i)       {return _state[i]; }
111
 
    const int& operator[](int i) const {return _state[i]; }
112
 
    const Mesh* GetMesh()                  const {return _mesh;}
113
 
 
114
 
    friend class Mesh;
115
 
    friend class TwoDGrid;
116
 
    friend class ThreeDGrid;
117
 
    friend class NDGrid;
118
 
    friend class TriangularMesh;
119
 
 
120
 
    friend Mesh::Iterator  operator++(Mesh::Iterator& lhs, int);
121
 
    friend Mesh::Iterator  operator--(Mesh::Iterator& lhs, int);
122
 
    friend Mesh::Iterator& operator++(Mesh::Iterator& lhs);
123
 
    friend Mesh::Iterator& operator--(Mesh::Iterator& lhs);
124
 
    friend Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const int& difference);
125
 
    friend Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const int& difference);
126
 
    friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs,       const int& difference);
127
 
    friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs,       const int& difference);
128
 
    friend Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
129
 
    friend Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
130
 
    friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
131
 
    friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
132
 
 
133
 
    friend bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
134
 
    friend bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
135
 
    friend bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
136
 
    friend bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
137
 
    friend bool operator< (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
138
 
    friend bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
139
 
 
140
 
 
141
 
private:
142
 
    const Mesh*      _mesh;
143
 
    std::vector<int> _state;
144
 
};
145
 
 
146
 
Mesh::Iterator& operator++(Mesh::Iterator& lhs); //no int means prefix operator ++it
147
 
Mesh::Iterator& operator--(Mesh::Iterator& lhs);
148
 
Mesh::Iterator  operator++(Mesh::Iterator& lhs, int); //int means postfix operator it++
149
 
Mesh::Iterator  operator--(Mesh::Iterator& lhs, int);
150
 
Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
151
 
Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
152
 
Mesh::Iterator& operator-=(Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
153
 
Mesh::Iterator& operator+=(Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
154
 
Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const int& rhs);
155
 
Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const int& rhs);
156
 
Mesh::Iterator& operator-=(Mesh::Iterator& lhs, const int& rhs);
157
 
Mesh::Iterator& operator+=(Mesh::Iterator& lhs, const int& rhs);
158
 
 
159
 
bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
160
 
bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
161
 
bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
162
 
bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
163
 
bool operator< (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
164
 
bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
165
 
 
166
 
std::ostream& operator<<(std::ostream& out, const Mesh::Iterator& it);
 
8
#include <ostream>
 
9
 
 
10
#include "src/legacy/Interface/Meshing/Mesh.hh"
 
11
#include "src/legacy/Interface/Meshing/ThreeDGrid.hh"
167
12
 
168
13
//Holds grid info for 2dTo1d interpolation algorithm; 
169
14
//also controls memory usage (essentially a boost::shared_pointer but not as elegent);
271
116
 
272
117
};
273
118
 
274
 
//Holds grid info for 2dTo1d interpolation algorithm; 
275
 
//also controls memory usage (essentially a boost::shared_pointer but not as elegant);
276
 
class ThreeDGrid : public Mesh
277
 
{
278
 
public:
279
 
    class Iterator;
280
 
    Mesh * Clone() {return new ThreeDGrid(*this);}
281
 
    Mesh * Dual () {return NULL;}
282
 
    //get value
283
 
    //ERROR SHOULD NOT ALLOW MESH WITH 1 GRID POINT (causes error in LowerBound)
284
 
    ThreeDGrid();
285
 
    ThreeDGrid(double dX, double dY, double dZ, double minX, double minY, double minZ, int numberOfXCoords, int numberOfYCoords, int numberOfZCoords);
286
 
    ThreeDGrid(int xSize, const double *x, int ySize, const double *y, int zSize, const double *z);
287
 
    ThreeDGrid(std::vector<double> x, std::vector<double> y, std::vector<double> z);
288
 
    ~ThreeDGrid();
289
 
 
290
 
    //get coordinate; round bracket indexing starts at 1 goes to nCoords
291
 
    inline double& x    (const int& i) {return _x[i-1];}
292
 
    inline double& y    (const int& j) {return _y[j-1];}
293
 
    inline double& z    (const int& k) {return _z[k-1];}
294
 
    inline const double& x    (const int& i) const {return _x[i-1];}
295
 
    inline const double& y    (const int& j) const {return _y[j-1];}
296
 
    inline const double& z    (const int& j) const {return _z[j-1];}
297
 
    inline int     xSize() const       {return int(_x.size());}
298
 
    inline int     ySize() const       {return int(_y.size());}
299
 
    inline int     zSize() const       {return int(_z.size());}
300
 
 
301
 
    std::vector<double> xVector()      {return std::vector<double>(_x);}
302
 
    std::vector<double> yVector()      {return std::vector<double>(_y);}
303
 
    std::vector<double> zVector()      {return std::vector<double>(_z);}
304
 
 
305
 
    double* newXArray()                {double *x = new double[_x.size()]; for(unsigned int i=0; i<_x.size(); i++) x[i]=_x[i]; return x;}
306
 
    double* newYArray()                {double *y = new double[_y.size()]; for(unsigned int i=0; i<_y.size(); i++) y[i]=_y[i]; return y;}
307
 
    double* newZArray()                {double *z = new double[_z.size()]; for(unsigned int i=0; i<_z.size(); i++) z[i]=_z[i]; return z;}
308
 
 
309
 
    //if you are sure the grid has constant spacing ConstSpacing is quicker; VarSpacing in either case
310
 
    //indexing starts at 0 and goes to nCoords+1
311
 
    inline void xLowerBound            (const double& x, int& xIndex) const 
312
 
    {if(_constantSpacing) xIndex = static_cast<int>(floor( (x - _x[0])/(_x[1]-_x[0]) )); else xIndex = std::lower_bound(_x.begin(), _x.end(), x)-_x.begin()-1;}
313
 
    inline void yLowerBound            (const double& y, int& yIndex) const
314
 
    {if(_constantSpacing) yIndex = static_cast<int>(floor( (y - _y[0])/(_y[1]-_y[0]) )); else yIndex = std::lower_bound(_y.begin(), _y.end(), y)-_y.begin()-1;}
315
 
    inline void zLowerBound            (const double& z, int& zIndex) const
316
 
    {if(_constantSpacing) zIndex = static_cast<int>(floor( (z - _z[0])/(_z[1]-_z[0]) )); else zIndex = std::lower_bound(_z.begin(), _z.end(), z)-_z.begin()-1;}
317
 
    inline void LowerBound             (const double& x, int& xIndex, const double& y, int& yIndex, const double& z, int& zIndex) const
318
 
    {xLowerBound(x, xIndex); yLowerBound(y, yIndex); zLowerBound(z, zIndex);}
319
 
    inline void LowerBound             (const double& x, const double& y, const double& z, Mesh::Iterator& it) const
320
 
    {xLowerBound(x, it[0]); yLowerBound(y, it[1]); zLowerBound(z, it[2]); it[0]++; it[1]++; it[2]++;}
321
 
    inline double    MinX()            const {return _x[0];}
322
 
    inline double    MaxX()            const {return _x[_xSize-1];}
323
 
    inline double    MinY()            const {return _y[0];}
324
 
    inline double    MaxY()            const {return _y[_ySize-1];}
325
 
    inline double    MinZ()            const {return _z[0];}
326
 
    inline double    MaxZ()            const {return _z[_zSize-1];}
327
 
 
328
 
    void    SetX(int nXCoords, double * x) {_x = std::vector<double>(x,x+nXCoords);}
329
 
    void    SetY(int nYCoords, double * y) {_y = std::vector<double>(y,y+nYCoords);}
330
 
    void    SetZ(int nZCoords, double * z) {_z = std::vector<double>(z,z+nZCoords);}
331
 
 
332
 
    void    Add(VectorMap* map); //add *map if it has not already been added
333
 
    void    Remove(VectorMap* map); //remove *map if it exists; delete this if there are no more VectorMaps
334
 
 
335
 
    Mesh::Iterator Begin() const;
336
 
    Mesh::Iterator End()   const;
337
 
    virtual void   Position(const Mesh::Iterator& it, double * position) const;
338
 
    int            PositionDimension() const {return 3;}
339
 
    int            ToInteger(const Mesh::Iterator& lhs) const  {return (lhs.State()[0]-1)*_zSize*_ySize+(lhs.State()[1]-1)*_zSize+(lhs.State()[2]-1);}
340
 
 
341
 
    void           SetConstantSpacing(bool spacing) {_constantSpacing = spacing;}
342
 
    void           SetConstantSpacing();
343
 
    bool           GetConstantSpacing() const       {return _constantSpacing; }
344
 
    Mesh::Iterator Nearest(const double* position) const;
345
 
 
346
 
protected:
347
 
 
348
 
    //Change position
349
 
    virtual Mesh::Iterator& AddEquals(Mesh::Iterator& lhs, int difference) const;
350
 
    virtual Mesh::Iterator& SubEquals(Mesh::Iterator& lhs, int difference) const;
351
 
    virtual Mesh::Iterator& AddEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
352
 
    virtual Mesh::Iterator& SubEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
353
 
    virtual Mesh::Iterator& AddOne   (Mesh::Iterator& lhs) const;
354
 
    virtual Mesh::Iterator& SubOne   (Mesh::Iterator& lhs) const;
355
 
    //Check relative position
356
 
    virtual bool IsGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
357
 
 
358
 
private:
359
 
    std::vector<double>     _x;
360
 
    std::vector<double>     _y;
361
 
    std::vector<double>     _z;
362
 
    int                     _xSize;
363
 
    int                     _ySize;
364
 
    int                     _zSize;
365
 
    std::vector<VectorMap*> _maps;
366
 
    bool                    _constantSpacing;
367
 
 
368
 
    friend Mesh::Iterator  operator++(Mesh::Iterator& lhs, int);
369
 
    friend Mesh::Iterator  operator--(Mesh::Iterator& lhs, int);
370
 
    friend Mesh::Iterator& operator++(Mesh::Iterator& lhs);
371
 
    friend Mesh::Iterator& operator--(Mesh::Iterator& lhs);
372
 
    friend Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
373
 
    friend Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
374
 
    friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
375
 
    friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
376
 
 
377
 
    friend bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
378
 
    friend bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
379
 
    friend bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
380
 
    friend bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
381
 
    friend bool operator< (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
382
 
    friend bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
383
 
};
384
 
 
385
 
 
386
119
//Holds grid info for some vector mapping; 
387
120
//also controls memory usage (essentially a boost::shared_pointer but not as elegent);
388
121
class NDGrid : public Mesh
416
149
 
417
150
    inline void    LowerBound      (const std::vector<double>& pos, std::vector<int>& xIndex) const
418
151
    {xIndex = std::vector<int>(pos.size()); for(unsigned int i=0; i<pos.size(); i++) coordLowerBound(pos[i], i, xIndex[i]);}
419
 
    inline void    LowerBound      (const MAUS::Vector<double>& pos,    std::vector<int>& xIndex) const
420
 
    {xIndex = std::vector<int>(pos.size()); for(size_t i=0; i<pos.size(); i++) coordLowerBound(pos(i+1), i, xIndex[i]);}
421
152
    inline double  Min             (const int& dimension) const {return _coord[dimension][0];}
422
153
    inline double  Max             (const int& dimension) const {return _coord[dimension][_coord[dimension].size()-1];}
423
154