← Back to branch summary

~christopher-hunt08/maus/beam_selection_development

« back to all changes in this revision

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

  • Committer: Chris Rogers
  • Date: 2012-10-13 18:43:00 UTC
  • mfrom: (663.6.137 merge)
  • mto: (663.6.204 merge)
  • mto: This revision was merged to the branch mainline in revision 680.
  • Revision ID: chris.rogers@stfc.ac.uk-20121013184300-ry9q81m45dmtgejr
Bring control room branch into line with trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/legacy/Interface/Mesh.cc
3
3
#include "Interface/Squeal.hh"
4
4
#include <iomanip>
5
5
#include <math.h>
6
 
///////// Mesh ////////
7
 
Mesh::Mesh() 
8
 
{}
9
 
 
10
 
Mesh::~Mesh() 
11
 
{}
12
 
 
13
 
/////// End of Mesh //////
14
6
 
15
7
/////// TwoDGrid ///////
16
8
TwoDGrid::TwoDGrid() : _x(2,0), _y(2,0), _xSize(2), _ySize(2), _maps(), _constantSpacing(false)
161
153
 
162
154
///// End of TwoDGrid /////
163
155
 
164
 
///// Mesh::Iterator /////
165
 
 
166
 
Mesh::Iterator::Iterator() {;}
167
 
Mesh::Iterator::Iterator(const Mesh::Iterator& in) : _mesh(in._mesh), _state(in._state)   {;}
168
 
Mesh::Iterator::Iterator(std::vector<int> state, const Mesh* mesh) : _mesh(mesh), _state(state) {;} 
169
 
Mesh::Iterator::~Iterator() {;}
170
 
 
171
 
const Mesh::Iterator&  Mesh::Iterator::operator= (const Mesh::Iterator& rhs)
172
 
{
173
 
    _mesh  = rhs._mesh;
174
 
    _state = rhs._state;
175
 
    return *this;
176
 
}
177
 
 
178
 
void Mesh::Iterator::Position(double* point) const
179
 
{
180
 
    _mesh->Position(*this, point);
181
 
}
182
 
 
183
 
std::vector<double> Mesh::Iterator::Position()              const
184
 
{
185
 
    std::vector<double> PointV(_mesh->PositionDimension());
186
 
    _mesh->Position(*this, &PointV[0]);
187
 
    return PointV;
188
 
}
189
 
 
190
 
Mesh::Iterator& operator++(Mesh::Iterator& lhs)      {return lhs._mesh->AddOne(lhs);}
191
 
Mesh::Iterator& operator--(Mesh::Iterator& lhs)      {return lhs._mesh->SubOne(lhs);}
192
 
Mesh::Iterator  operator++(Mesh::Iterator& lhs, int) {Mesh::Iterator copy = lhs; lhs._mesh->AddOne(lhs); return copy;}
193
 
Mesh::Iterator  operator--(Mesh::Iterator& lhs, int) {Mesh::Iterator copy = lhs; lhs._mesh->SubOne(lhs); return copy;}
194
 
 
195
 
Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) 
196
 
{
197
 
    Mesh::Iterator lhsCopy(lhs); 
198
 
    Mesh::Iterator lhsNew(lhs._mesh->SubEquals(lhsCopy, rhs)); 
199
 
    return lhsNew;
200
 
}
201
 
 
202
 
Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const int& difference) 
203
 
{
204
 
    Mesh::Iterator lhsCopy(lhs); 
205
 
    Mesh::Iterator lhsNew(lhs._mesh->SubEquals(lhsCopy, difference)); 
206
 
    return lhsNew;
207
 
}
208
 
 
209
 
Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) 
210
 
{
211
 
    Mesh::Iterator lhsCopy(lhs); 
212
 
    Mesh::Iterator lhsNew(lhs._mesh->AddEquals(lhsCopy, rhs)); 
213
 
    return lhsNew;
214
 
}
215
 
 
216
 
Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const int& difference) 
217
 
{
218
 
    Mesh::Iterator lhsCopy(lhs); 
219
 
    Mesh::Iterator lhsNew(lhs._mesh->AddEquals(lhsCopy, difference)); 
220
 
    return lhsNew;
221
 
}
222
 
 
223
 
 
224
 
Mesh::Iterator& operator+= (Mesh::Iterator& lhs, const Mesh::Iterator& rhs) { return lhs._mesh->AddEquals(lhs, rhs); }
225
 
Mesh::Iterator& operator-= (Mesh::Iterator& lhs, const Mesh::Iterator& rhs) { return lhs._mesh->SubEquals(lhs, rhs); }
226
 
Mesh::Iterator& operator+= (Mesh::Iterator& lhs, const int& rhs) { return lhs._mesh->AddEquals(lhs, rhs); }
227
 
Mesh::Iterator& operator-= (Mesh::Iterator& lhs, const int& rhs) { return lhs._mesh->SubEquals(lhs, rhs); }
228
 
 
229
 
 
230
 
bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs)
231
 
{
232
 
    if(lhs._mesh->IsGreater(lhs, rhs) || lhs==rhs) return true;
233
 
    return false;
234
 
}
235
 
 
236
 
bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs)
237
 
{
238
 
    if(lhs._mesh->IsGreater(lhs, rhs)) return false;
239
 
    return true;
240
 
}
241
 
 
242
 
bool operator>(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs)
243
 
{
244
 
    return !(lhs<=rhs);
245
 
}
246
 
 
247
 
bool operator<(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs)
248
 
{
249
 
    return !(lhs>=rhs);
250
 
}
251
 
 
252
 
bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) 
253
 
{
254
 
    if(lhs._state==rhs._state) 
255
 
        return true; 
256
 
    return false;
257
 
}
258
 
 
259
 
bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs)  
260
 
{
261
 
    if(lhs._mesh==rhs._mesh && lhs._state==rhs._state) 
262
 
        return false; 
263
 
    return true;
264
 
}
265
 
 
266
 
std::ostream& operator<<(std::ostream& out, const Mesh::Iterator& it)
267
 
{
268
 
    out << std::setw(5) << it.ToInteger() << " ** ";
269
 
    for(unsigned int i=0; i<it.State().size(); i++)    out << std::setw(5) << it[i] << " ";
270
 
    out << "** ";
271
 
    for(unsigned int i=0; i<it.Position().size(); i++) out << std::scientific << std::setprecision(3) << std::setw(12) << it.Position()[i] << " ";
272
 
    return out;
273
 
}
274
 
 
275
 
 
276
 
////// Mesh::Iterator end /////
277
 
 
278
 
////// ThreeDGrid //////
279
 
 
280
 
ThreeDGrid::ThreeDGrid() : _x(2, 0), _y(2, 0), _z(2, 0), _xSize(0), _ySize(0), _zSize(0), _maps(0), _constantSpacing(false) 
281
 
{
282
 
  SetConstantSpacing();
283
 
  _x[1] = _y[1] = _z[1] = 1.;
284
 
}
285
 
 
286
 
ThreeDGrid::ThreeDGrid(int xSize, const double *x, int ySize, const double *y, int zSize, const double *z) 
287
 
             : _x (x, x+xSize), _y(y, y+ySize), _z(z, z+zSize), _xSize(xSize), _ySize(ySize), _zSize(zSize), _maps(), _constantSpacing(false) 
288
 
{
289
 
  SetConstantSpacing();
290
 
  if(_xSize < 2 || _ySize < 2 || _zSize <2) throw(Squeal(Squeal::recoverable, "3D Grid must be at least 2x2x2 grid", "ThreeDGrid::ThreeDGrid(...)"));
291
 
}
292
 
 
293
 
ThreeDGrid::ThreeDGrid(std::vector<double> x, std::vector<double> y, std::vector<double> z)
294
 
             : _x (x), _y(y), _z(z), _xSize(x.size()), _ySize(y.size()), _zSize(z.size()), _maps(), _constantSpacing(false) 
295
 
{
296
 
  SetConstantSpacing();
297
 
  if(_xSize < 2 || _ySize < 2 || _zSize <2) throw(Squeal(Squeal::recoverable, "3D Grid must be at least 2x2x2 grid", "ThreeDGrid::ThreeDGrid(...)"));
298
 
}
299
 
 
300
 
ThreeDGrid::ThreeDGrid(double dX, double dY, double dZ, double minX, double minY, double minZ, 
301
 
                       int numberOfXCoords, int numberOfYCoords, int numberOfZCoords)
302
 
        :  _x(numberOfXCoords), _y(numberOfYCoords), _z(numberOfZCoords), _xSize(numberOfXCoords), _ySize(numberOfYCoords), _zSize(numberOfZCoords), _maps(), _constantSpacing(true)
303
 
{
304
 
    for(int i=0; i<numberOfXCoords; i++) _x[i] = minX+i*dX;
305
 
    for(int j=0; j<numberOfYCoords; j++) _y[j] = minY+j*dY;
306
 
    for(int k=0; k<numberOfZCoords; k++) _z[k] = minZ+k*dZ;
307
 
 
308
 
    SetConstantSpacing();
309
 
}
310
 
 
311
 
ThreeDGrid::~ThreeDGrid() 
312
 
{
313
 
}
314
 
 
315
 
//state starts at 1,1,1
316
 
Mesh::Iterator ThreeDGrid::Begin() const 
317
 
{
318
 
    return Mesh::Iterator(std::vector<int>(3,1), this);
319
 
}
320
 
 
321
 
Mesh::Iterator ThreeDGrid::End()   const 
322
 
{
323
 
    std::vector<int> end(3, 1); end[0] = _xSize+1; return Mesh::Iterator(end, this);
324
 
}
325
 
 
326
 
void ThreeDGrid::Position(const Mesh::Iterator& it, double * position) const
327
 
{
328
 
    position[0] = x(it._state[0]);
329
 
    position[1] = y(it._state[1]);
330
 
    position[2] = z(it._state[2]);
331
 
}
332
 
 
333
 
 
334
 
Mesh::Iterator& ThreeDGrid::AddEquals(Mesh::Iterator& lhs, int difference) const
335
 
{
336
 
    if(difference < 0) return SubEquals(lhs, -difference);
337
 
 
338
 
    lhs._state[0] +=  difference/(_ySize*_zSize);
339
 
    difference     = difference%(_ySize*_zSize);
340
 
    lhs._state[1] += difference/(_zSize);
341
 
    lhs._state[2] += difference%(_zSize);
342
 
 
343
 
    if(lhs._state[1] > _ySize) {lhs._state[0]++; lhs._state[1] -= _ySize;}
344
 
    if(lhs._state[2] > _zSize) {lhs._state[1]++; lhs._state[2] -= _zSize;}
345
 
 
346
 
    return lhs;
347
 
}
348
 
/*
349
 
Mesh::Iterator& TwoDGrid::SubEquals(Mesh::Iterator& lhs, int difference) const
350
 
{
351
 
    if(difference < 0) return AddEquals(lhs, -difference);
352
 
    lhs._state[0] -= difference/(_ySize);
353
 
    lhs._state[1] -= difference%(_ySize);
354
 
    if(lhs._state[1] < 1) {lhs._state[1] += _ySize; lhs._state[0]--;}
355
 
    return lhs;
356
 
}
357
 
*/
358
 
Mesh::Iterator& ThreeDGrid::SubEquals(Mesh::Iterator& lhs, int difference) const
359
 
{
360
 
    if(difference < 0) return AddEquals(lhs, -difference);
361
 
 
362
 
    lhs._state[0] -= difference/(_ySize*_zSize);
363
 
    difference     = difference%(_ySize*_zSize);
364
 
    lhs._state[1] -= difference/(_zSize);
365
 
    lhs._state[2] -= difference%(_zSize);
366
 
 
367
 
    while(lhs._state[2] < 1) {lhs._state[1]--; lhs._state[2] += _zSize;}
368
 
    while(lhs._state[1] < 1) {lhs._state[0]--; lhs._state[1] += _ySize;}
369
 
 
370
 
    return lhs;
371
 
}
372
 
 
373
 
 
374
 
Mesh::Iterator& ThreeDGrid::AddEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const
375
 
{
376
 
    return AddEquals(lhs, rhs.ToInteger());
377
 
}
378
 
 
379
 
Mesh::Iterator& ThreeDGrid::SubEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const
380
 
{
381
 
    return SubEquals(lhs, rhs.ToInteger());
382
 
}
383
 
 
384
 
Mesh::Iterator& ThreeDGrid::AddOne   (Mesh::Iterator& lhs) const
385
 
{
386
 
    if(lhs._state[1] == _ySize && lhs._state[2] == _zSize) {++lhs._state[0]; lhs._state[1] = lhs._state[2] = 1;}
387
 
    else if(lhs._state[2] == _zSize) {++lhs._state[1]; lhs._state[2] = 1;}
388
 
    else ++lhs._state[2];
389
 
    return lhs;
390
 
}
391
 
 
392
 
Mesh::Iterator& ThreeDGrid::SubOne   (Mesh::Iterator& lhs) const
393
 
{
394
 
    if(lhs._state[1] == 1 && lhs._state[2] == 1) {--lhs._state[0]; lhs._state[1] = _ySize; lhs._state[2] = _zSize;}
395
 
    else if(lhs._state[2] == 1) {--lhs._state[1]; lhs._state[2] = _zSize;}
396
 
    else --lhs._state[2];
397
 
    return lhs;
398
 
}
399
 
//Check relative position
400
 
bool ThreeDGrid::IsGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const
401
 
{
402
 
    if(lhs._state[0] > rhs._state[0]) return true;
403
 
    else if (lhs._state[0] == rhs._state[0] && lhs._state[1] > rhs._state[1]) return true;
404
 
    else if (lhs._state[0] == rhs._state[0] && lhs._state[1] == rhs._state[1] && lhs._state[2] > rhs._state[2]) return true;
405
 
    return false;
406
 
}
407
 
 
408
 
void ThreeDGrid::Remove(VectorMap* map) //remove *map if it exists; delete this if there are no more VectorMaps
409
 
{
410
 
    std::vector<VectorMap*>::iterator it = std::find(_maps.begin(), _maps.end(), map); 
411
 
    if(it<_maps.end()) { _maps.erase(it);  }
412
 
    if(_maps.begin() == _maps.end()) { delete this; }
413
 
} 
414
 
 
415
 
void ThreeDGrid::Add(VectorMap* map) //add *map if it has not already been added
416
 
{
417
 
    std::vector<VectorMap*>::iterator it = std::find(_maps.begin(), _maps.end(), map); 
418
 
    if(it==_maps.end()) { _maps.push_back(map);    }
419
 
}
420
 
 
421
 
void ThreeDGrid::SetConstantSpacing()
422
 
{
423
 
    _constantSpacing = true;
424
 
    for(unsigned int i=0; i<_x.size()-1; i++) if( fabs(1-(_x[i+1]-_x[i])/(_x[1]-_x[0])) > 1e-9 ) {_constantSpacing = false; return;}
425
 
    for(unsigned int i=0; i<_y.size()-1; i++) if( fabs(1-(_y[i+1]-_y[i])/(_y[1]-_y[0])) > 1e-9 ) {_constantSpacing = false; return;}
426
 
    for(unsigned int i=0; i<_z.size()-1; i++) if( fabs(1-(_z[i+1]-_z[i])/(_z[1]-_z[0])) > 1e-9 ) {_constantSpacing = false; return;}
427
 
}
428
 
 
429
 
Mesh::Iterator ThreeDGrid::Nearest(const double* position) const
430
 
{
431
 
    std::vector<int> index(3);
432
 
    LowerBound(position[0], index[0], position[1], index[1], position[2], index[2]);
433
 
    if(index[0] < _xSize-1 && index[0]>=0) index[0] += (2*(position[0] - _x[index[0]])  > _x[index[0]+1]-_x[index[0]] ? 2 : 1);
434
 
    else index[0]++;
435
 
    if(index[1] < _ySize-1 && index[1]>=0) index[1] += (2*(position[1] - _y[index[1]])  > _y[index[1]+1]-_y[index[1]] ? 2 : 1);
436
 
    else index[1]++;
437
 
    if(index[2] < _zSize-1 && index[2]>=0) index[2] += (2*(position[2] - _z[index[2]])  > _z[index[2]+1]-_z[index[2]] ? 2 : 1);
438
 
    else index[2]++;
439
 
    if(index[0] < 1) index[0] = 1; 
440
 
    if(index[1] < 1) index[1] = 1; 
441
 
    if(index[2] < 1) index[2] = 1; 
442
 
    if(index[0] > _xSize) index[0] = _xSize; 
443
 
    if(index[1] > _ySize) index[1] = _ySize; 
444
 
    if(index[2] > _zSize) index[2] = _zSize; 
445
 
    return Mesh::Iterator(index, this);
446
 
}
447
 
 
448
 
/// ////// ThreeDGrid END //////
449
156
 
450
157
/// ////// NDGrid ///////
451
158
NDGrid::NDGrid() : _coord(), _maps(), _constantSpacing(false)