ESyS-Particle  2.3.2
ParticleFitter.h
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 #ifndef ESYS_LSMPARTICLEFITTER_H
14 #define ESYS_LSMPARTICLEFITTER_H
15 
17 #include "Geometry/Sphere3d.h"
18 #include "Geometry/Sphere2d.h"
19 
20 namespace esys
21 {
22  namespace lsm
23  {
25  {
26  public:
28 
30  : m_pGenerator(&blockGenerator),
32  m_getFitCount(0),
34  {
35  }
36 
37  virtual ~ParticleFitter() {}
38 
40  const SimpleParticle &particle,
41  const ParticleVector &neighbours,
42  const Plane3D &plane
43  ) = 0;
44 
45  static const SimpleParticle INVALID;
46 
47  void incrGetFit()
48  {
49  m_getFitCount++;
50  }
51 
53  {
55  }
56 
58  {
60  }
61 
62  virtual std::string getName() const = 0;
63 
64  void write(std::ostream &oStream) const
65  {
66  oStream
67  << (getName() + ": ")
68  << "fit count = " << m_getFitCount
69  << ", success = " << m_successfulFitCount
70  << ", fail = " << m_failedFitCount;
71  }
72 
73  std::string toString() const
74  {
75  std::stringstream sStream;
76  write(sStream);
77  return sStream.str();
78  }
79 
80  virtual bool particleFits(const SimpleParticle &particle) const
81  {
82  return getGenerator().particleFits(particle);
83  }
84 
85  protected:
87  {
88  return *m_pGenerator;
89  }
90 
92  {
93  return *m_pGenerator;
94  }
95  private:
100  };
101 
102  inline std::ostream &operator<<(std::ostream &oStream, const ParticleFitter &fitter)
103  {
104  fitter.write(oStream);
105  return oStream;
106  }
107 
109  {
110  public:
112  : ParticleFitter(blockGenerator)
113  {
114  }
115 
116  virtual std::string getName() const
117  {
118  return "Mv to Surface";
119  }
120 
122  const SimpleParticle &stationary,
123  const SimpleParticle &movable
124  )
125  {
126  SimpleParticle moved = movable;
127  const Vec3 centreDiff = movable.getPos() - stationary.getPos();
128  const double centreDist = centreDiff.norm();
129  if (centreDist > 0.0) {
130  const Vec3 newCentrePos =
131  stationary.getPos()
132  +
133  (centreDiff/(centreDist))*(stationary.getRad() + movable.getRad());
134  moved.moveTo(newCentrePos);
135  }
136  return moved;
137  }
138 
140  const SimpleParticle &particle,
141  const ParticleVector &neighbours,
142  const Plane3D &plane
143  )
144  {
145  incrGetFit();
146  SimpleParticle newParticle = particle;
147  if (neighbours.size() > 0) {
148  if (
149  (particle.getPos() - neighbours[0]->getPos()).norm()
150  <
151  (particle.getRad() + neighbours[0]->getRad())
152  ) {
153  newParticle = moveToSurface(*(neighbours[0]), particle);
154  }
155  }
156  if (newParticle.isValid() && !particleFits(newParticle)) {
157  newParticle = INVALID;
158  incrFailedFit();
159  } else if (newParticle.isValid()) {
161  }
162  return newParticle;
163  }
164  };
165 
167  {
168  public:
170  : ParticleFitter(blockGenerator)
171  {
172  }
173 
174  virtual std::string getName() const
175  {
176  return " 3D Fitter";
177  }
178 
180  const SimpleParticle& Po,
181  const ParticleVector &particleVector
182  )
183  {
184  SimpleParticle fittedParticle = SimpleParticle::INVALID;
185  Vec3 M;
186  double r;
187  int id=Po.getID();
188 
189  if (particleVector.size() > 3) {
190  const bool foundFit =
192  particleVector[0]->getPos(),
193  particleVector[1]->getPos(),
194  particleVector[2]->getPos(),
195  particleVector[3]->getPos(),
196  particleVector[0]->getRad(),
197  particleVector[1]->getRad(),
198  particleVector[2]->getRad(),
199  particleVector[3]->getRad(),
200  M,
201  r
202  );
203  if (foundFit) {
204  fittedParticle = SimpleParticle(M, r, id);
205  }
206  } else {
207  throw std::runtime_error("findAFit: particleVector argument has fewer than 4 elements.");
208  }
209  return fittedParticle;
210  }
211 
213  const SimpleParticle &particle,
214  const ParticleVector &neighbours,
215  const Plane3D &plane
216  )
217  {
218  incrGetFit();
219  SimpleParticle newParticle = INVALID;
220  if(neighbours.size() > 3){ // at least 4 neighbors
221  const SimpleParticle &Pi = *(neighbours[0]);
222  const double ndist=(particle.getPos()-Pi.getPos()).norm();
223  if (ndist > 0.0) {
224  newParticle = particle;
225  if (ndist < Pi.getRad()){ // if Po inside Pi -> move Po to the surface of Pi
226  newParticle.moveTo(
227  Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist))
228  );
229  }
230  const double dist_p = plane.sep(newParticle.getPos());
231  const double dist_3 = (newParticle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad();
232  if (dist_p > dist_3) { // 4th particle closer than plane -> fit 4 particles
233  newParticle = findAFit(newParticle, neighbours);
234  }
235  else {
236  newParticle = INVALID;
237  }
238  }
239  }
240  if (newParticle.isValid() && !particleFits(newParticle)) {
241  newParticle = INVALID;
242  incrFailedFit();
243  } else if (newParticle.isValid()) {
245  }
246  return newParticle;
247  }
248  };
249 
251  {
252  public:
254  : ParticleFitter(blockGenerator)
255  {
256  }
257 
258  virtual std::string getName() const
259  {
260  return " 2D Fitter";
261  }
262 
264  const SimpleParticle& Po,
265  const ParticleVector &particleVector
266  )
267  {
268  SimpleParticle fittedParticle = SimpleParticle::INVALID;
269  Vec3 M;
270  double r;
271  int id=Po.getID();
272 
273  if (particleVector.size() > 2) {
274  const bool foundFit =
276  particleVector[0]->getPos(),
277  particleVector[1]->getPos(),
278  particleVector[2]->getPos(),
279  particleVector[0]->getRad(),
280  particleVector[1]->getRad(),
281  particleVector[2]->getRad(),
282  M,
283  r
284  );
285  if (foundFit) {
286  fittedParticle = SimpleParticle(M, r, id);
287  }
288  } else {
289  throw std::runtime_error("findAFit: particleVector argument has fewer than 3 elements.");
290  }
291  return fittedParticle;
292  }
293 
295  const SimpleParticle &particle,
296  const ParticleVector &neighbours,
297  const Plane3D &plane
298  )
299  {
300  incrGetFit();
301  SimpleParticle newParticle = INVALID;
302  if(neighbours.size() > 2){ // at least 3 neighbors
303  const SimpleParticle &Pi = *(neighbours[0]);
304  const double ndist=(particle.getPos()-Pi.getPos()).norm();
305  if (ndist > 0.0) {
306  newParticle = particle;
307  if (ndist < Pi.getRad()){ // if Po inside Pi -> move Po to the surface of Pi
308  newParticle.moveTo(
309  Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist))
310  );
311  }
312  const double dist_p = plane.sep(newParticle.getPos());
313  const double dist_2 = (newParticle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad();
314  if (dist_p > dist_2) { // 4th particle closer than plane -> fit 4 particles
315  newParticle = findAFit(newParticle, neighbours);
316  }
317  else {
318  newParticle = INVALID;
319  }
320  }
321  }
322  if (newParticle.isValid() && !particleFits(newParticle)) {
323  newParticle = INVALID;
324  incrFailedFit();
325  } else if (newParticle.isValid()) {
327  }
328  return newParticle;
329  }
330  };
331 
333  {
334  public:
336  : ParticleFitter(blockGenerator)
337  {
338  }
339 
340  virtual std::string getName() const
341  {
342  return " 2D Plane";
343  }
344 
346  const SimpleParticle& particle,
347  const ParticleVector &particleVector,
348  const Plane3D& plane
349  )
350  {
351  SimpleParticle fittedParticle = SimpleParticle::INVALID;
352  Vec3 M;
353  double r;
354  const int id = particle.getID();
355 
356  if (particleVector.size() > 1) {
357  const bool foundFit =
359  particleVector[0]->getPos(),
360  particleVector[1]->getPos(),
361  plane.GetO(),
362  plane.GetW(),
363  particleVector[0]->getRad(),
364  particleVector[1]->getRad(),
365  M,
366  r
367  );
368  if (foundFit) {
369  fittedParticle = SimpleParticle(M, r, id);
370  }
371  } else {
372  throw std::runtime_error("findAFit: particleVector vector contains less than 2 particles.");
373  }
374 
375  return fittedParticle;
376  }
377 
379  const SimpleParticle &particle,
380  const ParticleVector &neighbours,
381  const Plane3D &plane
382  )
383  {
384  incrGetFit();
385  SimpleParticle newParticle = INVALID;
386  if (neighbours.size() >= 2) { // 2 neighbors -> try 2 particles + plane
387  const SimpleParticle &Pi = *(neighbours[0]);
388  const double ndist=(particle.getPos()-Pi.getPos()).norm();
389  if (
390  (ndist > 0.0)
391  &&
392  (
393  (neighbours.size() == 2)
394  ||
395  (
396  plane.sep(particle.getPos())
397  <
398  (particle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad()
399  )
400  )
401  ) {
402  newParticle = particle;
403  if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi
404  newParticle.moveTo(
405  Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
406  );
407  }
408  newParticle = findAFit(newParticle, neighbours, plane);
409  }
410  }
411  if (newParticle.isValid() && !particleFits(newParticle)) {
412  newParticle = INVALID;
413  incrFailedFit();
414  } else if (newParticle.isValid()) {
416  }
417  return newParticle;
418  }
419  };
420 
422  {
423  public:
425  : ParticleFitter(blockGenerator)
426  {
427  }
428 
429  virtual std::string getName() const
430  {
431  return " 3D Plane";
432  }
433 
435  const SimpleParticle& particle,
436  const ParticleVector &particleVector,
437  const Plane3D& plane
438  )
439  {
440  SimpleParticle fittedParticle = SimpleParticle::INVALID;
441  Vec3 M;
442  double r;
443  const int id = particle.getID();
444 
445  if (particleVector.size() > 2) {
446  const bool foundFit =
448  particleVector[0]->getPos(),
449  particleVector[1]->getPos(),
450  particleVector[2]->getPos(),
451  plane.GetO(),
452  plane.GetW(),
453  particleVector[0]->getRad(),
454  particleVector[1]->getRad(),
455  particleVector[2]->getRad(),
456  M,
457  r
458  );
459  if (foundFit) {
460  fittedParticle = SimpleParticle(M, r, id);
461  }
462  } else {
463  throw std::runtime_error("findAFit: particleVector vector contains less than 3 particles.");
464  }
465 
466  return fittedParticle;
467  }
468 
470  const SimpleParticle &particle,
471  const ParticleVector &neighbours,
472  const Plane3D &plane
473  )
474  {
475  incrGetFit();
476  SimpleParticle newParticle = INVALID;
477  if (neighbours.size() >= 3) { // 3 neighbors -> try 3 particles + plane
478  const SimpleParticle &Pi = *(neighbours[0]);
479  const double ndist=(particle.getPos()-Pi.getPos()).norm();
480  if (
481  (ndist > 0.0)
482  &&
483  (
484  (neighbours.size() == 3)
485  ||
486  (
487  plane.sep(particle.getPos())
488  <
489  (particle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad()
490  )
491  )
492  ) {
493  newParticle = particle;
494  if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi
495  newParticle.moveTo(
496  Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
497  );
498  }
499  newParticle = findAFit(newParticle, neighbours, plane);
500  }
501  }
502  if (newParticle.isValid() && !particleFits(newParticle)) {
503  newParticle = INVALID;
504  incrFailedFit();
505  } else if (newParticle.isValid()) {
507  }
508  return newParticle;
509  }
510  };
511  }
512 }
513 
514 #endif