MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
IntersectionOfWalls.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "IntersectionOfWalls.h"
27 #include "InteractionHandler.h"
28 #include "WallHandler.h"
29 #include "DPMBase.h"
30 //#include "Species/BaseSpecies.h"
31 #include "Particles/BaseParticle.h"
32 
34 {
35  logger(DEBUG, "IntersectionOfWalls() constructed.");
36 }
37 
42  : BaseWall(other)
43 {
44  wallObjects_ = other.wallObjects_;
45  for (auto& wall : wallObjects_)
46  {
47  if (getHandler() != nullptr)
48  wall.setHandler(getHandler());
49  }
50  A_ = other.A_;
51  AB_ = other.AB_;
52  C_ = other.C_;
53  logger(DEBUG, "IntersectionOfWalls(IntersectionOfWalls&) constructed.");
54 }
55 
56 IntersectionOfWalls::IntersectionOfWalls(const std::vector<normalAndPosition>& walls,
57  const ParticleSpecies* species)
58 {
59  setSpecies(species);
60  for (auto wall : walls)
61  {
62  addObject(wall.normal, wall.position);
63  }
64 }
65 
66 
68 {
69  logger(DEBUG, "~IntersectionOfWalls() has been called.");
70 }
71 
73 {
74  BaseWall::setSpecies(species);
75  for (auto wall : wallObjects_)
76  {
77  wall.setSpecies(species);
78  }
79 }
80 
85 {
86  logger(DEBUG, "IntersectionOfWalls::operator= called.");
87  if (this == &other)
88  {
89  return *this;
90  }
91  return *(other.copy());
92 }
93 
98 {
99  return new IntersectionOfWalls(*this);
100 }
101 
103 {
104  wallObjects_.clear();
105  A_.clear();
106  AB_.clear();
107  C_.clear();
108 }
109 
111 {
112  BaseWall::setHandler(wallHandler);
113  for (InfiniteWall& w : wallObjects_)
114  {
115  w.setHandler(wallHandler);
116  }
117 }
118 
125 {
126  return wallObjects_.size();
127 }
128 
138 {
139  normal.normalise();
140 
141  //n is the index of the new wall
142  std::size_t n = wallObjects_.size();
143  InfiniteWall w;
144  if (getSpecies()) w.setSpecies(getSpecies());
145  w.set(normal, point);
146  wallObjects_.push_back(w);
148 }
149 
159 {
160  Vec3D SubtB = PointA - PointB;
161  Vec3D SubtC = PointA - PointC;
162 
163  Vec3D WallNormal = Vec3D::cross(SubtB, SubtC);
164  //Check if walls coordinates inline, if true Do not build wall and give error message.
165  if (WallNormal.getLengthSquared() == 0.0)
166  {
167  logger(ERROR,
168  "Error Building IntersectionOfWalls::add3PointObject out of 3 coordinates. Coordinates are in line, Wall not constructed.");
169  }
170  else
171  {
172  addObject(WallNormal, PointA);
173  }
174 
175 }
176 
177 
190 void IntersectionOfWalls::addTetraSTL(Vec3D PointA, Vec3D PointB, Vec3D PointC, Vec3D WallNormal, Mdouble Thickness,
191  int wallidentifier)
192 {
193  //Check if wall coordinates inline, if true Do not build wall and give error message. But keep continuing
194  if (WallNormal.getLengthSquared() == 0.0)
195  {
196  std::cout << "Error Building Plate number " << wallidentifier
197  << "out of 3 coordinates. Coordinates are in line, Wall not constructed." << std::endl;
198  }
199  else
200  {
201  //do a check whether the normal follows the RHR (Right Hand Rule) or not
202  //todo: Bert Use this check a lot, possibly make a function of this
203  Vec3D SubtB = PointA - PointB;
204  Vec3D SubtC = PointA - PointC;
205 
206  Vec3D WallNormalRHR = Vec3D::cross(SubtB, SubtC);
207 
208  //normalise for easy check
209  WallNormalRHR.normalise();
210  WallNormal.normalise();
211 
212  //if RHRchecl is 1, wall normal and RHR normal are in same direction, if -1 theyre not, then point B and C need to be swapped
213  Mdouble RHRcheck = Vec3D::dot(WallNormalRHR, WallNormal);
214  //todo: Bert Officially need to check for other answers than 1, however, it will either be -1 or 1 in ideal case
215 
216  if (RHRcheck == 1)
217  {
218  //calculate centroid
219  Vec3D mid = (PointA + PointB + PointC) / 3.0;
220  //shift centroid in normal direction
221  //mental note: if Same direction as STL normal it should be subtracted to go IN the wall
222  Vec3D midT = mid - WallNormalRHR * Thickness;
223  //generate the base wall first through the input coordinates
224  add3PointObject(PointA, PointC, PointB);
225 
226  //add sidewalls
227  add3PointObject(PointA, midT, PointC);
228  add3PointObject(PointC, midT, PointB);
229  add3PointObject(PointB, midT, PointA);
230  }
231  else
232  {
233  //calculate centroid
234  Vec3D mid = (PointA + PointB + PointC) / 3.0;
235  //shift centroid in normal direction
236  //mental note: if opposite direction as STL normal it should be added to go IN the wall
237  Vec3D midT = mid + WallNormalRHR * Thickness;
238  //generate the base wall first through the input coordinates
239  add3PointObject(PointA, PointB, PointC);
240 
241  //add sidewalls
242  add3PointObject(PointA, midT, PointB);
243  add3PointObject(PointB, midT, PointC);
244  add3PointObject(PointC, midT, PointA);
245  }
246 
247  }
248 }
249 
250 
259 void IntersectionOfWalls::addTetra(const Vec3D& PointA, const Vec3D& PointB, const Vec3D& PointC, Mdouble& Thickness)
260 {
261 
262  //generate other coordinate by finding the centre and shift it the distance thickness in the normal direction
263  //calculate normal
264  Vec3D SubtB = PointA - PointB;
265  Vec3D SubtC = PointA - PointC;
266 
267  Vec3D WallNormal = Vec3D::cross(SubtB, SubtC);
268 
269  //Check if walls coordinates inline, if true Do not build wall and give error message.
270  if (WallNormal.getLengthSquared() == 0.0)
271  {
272  logger(ERROR,
273  "Error Building IntersectionOfWalls::addTetra out of 3 coordinates. "
274  "Coordinates are in line, Wall not constructed.");
275  }
276  else
277  {
278  //calculate centroid
279  Vec3D mid = (PointA + PointB + PointC) / 3.0;
280  //shift centroid in normal direction
281  Vec3D midT = mid + WallNormal * Thickness;
282  //generate the base wall first through the input coordinates
283  add3PointObject(PointC, PointB, PointA);
284 
285  //add sidewalls
286  add3PointObject(PointA, midT, PointB);
287  add3PointObject(PointB, midT, PointC);
288  add3PointObject(PointC, midT, PointA);
289  }
290 }
291 
292 
293 //funtion builds a triangle wall through the 3 points with a thickness
294 //note normal is defined OPPOSITE as to normal Mercury convention due to STL convention
295 void
296 IntersectionOfWalls::addPlate(const Vec3D& PointA, const Vec3D& PointB, const Vec3D& PointC, const Vec3D& WallNormal,
297  const Mdouble& Thickness, int wallidentifier)
298 {
299 
300  //Check if walls coordinates inline, if true Do not build wall and give error message. But keep continuing
301  if (WallNormal.getLengthSquared() == 0.0)
302  {
303  std::cout << "Error Building Plate number " << wallidentifier
304  << "out of 3 coordinates. Coordinates are in line, Wall not constructed." << std::endl;
305  }
306  else
307  {
308  Vec3D PointAT = PointA - WallNormal * Thickness;
309  Vec3D PointBT = PointB - WallNormal * Thickness;
310  Vec3D PointCT = PointC - WallNormal * Thickness;
311 
312  //generate the base wall first through the input coordinates
313  add3PointObject(PointC, PointB, PointA);
314 
315  //add sidewalls
316  add3PointObject(PointA, PointB, PointAT);
317  add3PointObject(PointB, PointC, PointBT);
318  add3PointObject(PointC, PointA, PointCT);
319 
320 
321  //add opposite wall
322  add3PointObject(PointAT, PointBT, PointCT);
323 
324 
325  }
326 
327 
328 }
329 
339 {
340  //n is the index of the new wall
341  std::size_t n = wallObjects_.size();
342  InfiniteWall w;
343  if (getSpecies()) w.setSpecies(getSpecies());
344  w.setOrientation(orientation);
345  w.setPosition(position);
346  wallObjects_.push_back(w);
348 }
349 
351 {
352  // AB[n*(n-1)/2+m] is the direction of the intersecting line between walls m and n, m<n
353  // A[n*(n-1)/2+m] is a point on the intersecting line between walls m and n, m<n
354  // See http://www.netcomuk.co.uk/~jenolive/vect18d.html for finding the line where two planes meet
355  AB_.resize(n * (n + 1) / 2); //0 + 1 + 2 + ... + indexNew, total number of walls you need
356  A_.resize(n * (n + 1) / 2);
357  for (std::size_t m = 0; m < n; m++)
358  {
359  std::size_t id = (n - 1) * n / 2 + m;
360  //first we cross the wall normals and normalise to obtain AB
361  AB_[id] = Vec3D::cross(wallObjects_[m].getNormal(), wallObjects_[n].getNormal());
362  AB_[id] /= sqrt(AB_[id].getLengthSquared());
363  //then we find a point A (using AB*x=0 as a third plane)
364  Mdouble invdet = 1.0 / (+wallObjects_[n].getNormal().X *
365  (wallObjects_[m].getNormal().Y * AB_[id].Z - AB_[id].Y * wallObjects_[m].getNormal().Z)
366  - wallObjects_[n].getNormal().Y * (wallObjects_[m].getNormal().X * AB_[id].Z -
367  wallObjects_[m].getNormal().Z * AB_[id].X)
368  + wallObjects_[n].getNormal().Z * (wallObjects_[m].getNormal().X * AB_[id].Y -
369  wallObjects_[m].getNormal().Y * AB_[id].X));
370 
371  A_[id] = Vec3D(+(wallObjects_[m].getNormal().Y * AB_[id].Z - AB_[id].Y * wallObjects_[m].getNormal().Z) *
372  Vec3D::dot(wallObjects_[n].getPosition(), wallObjects_[n].getNormal())
373  - (wallObjects_[n].getNormal().Y * AB_[id].Z - wallObjects_[n].getNormal().Z * AB_[id].Y) *
374  Vec3D::dot(wallObjects_[m].getPosition(), wallObjects_[m].getNormal())
375  + (wallObjects_[n].getNormal().Y * wallObjects_[m].getNormal().Z -
376  wallObjects_[n].getNormal().Z * wallObjects_[m].getNormal().Y) * 0.0,
377  -(wallObjects_[m].getNormal().X * AB_[id].Z - wallObjects_[m].getNormal().Z * AB_[id].X) *
378  Vec3D::dot(wallObjects_[n].getPosition(), wallObjects_[n].getNormal())
379  + (wallObjects_[n].getNormal().X * AB_[id].Z - wallObjects_[n].getNormal().Z * AB_[id].X) *
380  Vec3D::dot(wallObjects_[m].getPosition(), wallObjects_[m].getNormal())
381  - (wallObjects_[n].getNormal().X * wallObjects_[m].getNormal().Z -
382  wallObjects_[m].getNormal().X * wallObjects_[n].getNormal().Z) * 0.0,
383  +(wallObjects_[m].getNormal().X * AB_[id].Y - AB_[id].X * wallObjects_[m].getNormal().Y) *
384  Vec3D::dot(wallObjects_[n].getPosition(), wallObjects_[n].getNormal())
385  - (wallObjects_[n].getNormal().X * AB_[id].Y - AB_[id].X * wallObjects_[n].getNormal().Y) *
386  Vec3D::dot(wallObjects_[m].getPosition(), wallObjects_[m].getNormal())
387  + (wallObjects_[n].getNormal().X * wallObjects_[m].getNormal().Y -
388  wallObjects_[m].getNormal().X * wallObjects_[n].getNormal().Y) * 0.0) * invdet;
389  }
390 
391  // C[(n-2)*(n-1)*n/6+(m-1)*m/2+l] is a point intersecting walls l, m and n, l<m<n
392  C_.resize((n - 1) * n * (n + 1) / 6);
393  for (std::size_t m = 0; m < n; m++)
394  {
395  for (std::size_t l = 0; l < m; l++)
396  {
397  std::size_t id = (n - 2) * (n - 1) * n / 6 + (m - 1) * m / 2 + l;
398  Mdouble invdet = 1.0 / (+wallObjects_[n].getNormal().X *
399  (wallObjects_[m].getNormal().Y * wallObjects_[l].getNormal().Z -
400  wallObjects_[l].getNormal().Y * wallObjects_[m].getNormal().Z)
401  - wallObjects_[n].getNormal().Y *
402  (wallObjects_[m].getNormal().X * wallObjects_[l].getNormal().Z -
403  wallObjects_[m].getNormal().Z * wallObjects_[l].getNormal().X)
404  + wallObjects_[n].getNormal().Z *
405  (wallObjects_[m].getNormal().X * wallObjects_[l].getNormal().Y -
406  wallObjects_[m].getNormal().Y * wallObjects_[l].getNormal().X));
407  C_[id] = Vec3D(+(wallObjects_[m].getNormal().Y * wallObjects_[l].getNormal().Z -
408  wallObjects_[l].getNormal().Y * wallObjects_[m].getNormal().Z) *
409  Vec3D::dot(wallObjects_[n].getPosition(), wallObjects_[n].getNormal())
410  - (wallObjects_[n].getNormal().Y * wallObjects_[l].getNormal().Z -
411  wallObjects_[n].getNormal().Z * wallObjects_[l].getNormal().Y) *
412  Vec3D::dot(wallObjects_[m].getPosition(), wallObjects_[m].getNormal())
413  + (wallObjects_[n].getNormal().Y * wallObjects_[m].getNormal().Z -
414  wallObjects_[n].getNormal().Z * wallObjects_[m].getNormal().Y) *
415  Vec3D::dot(wallObjects_[l].getPosition(), wallObjects_[l].getNormal()),
416  -(wallObjects_[m].getNormal().X * wallObjects_[l].getNormal().Z -
417  wallObjects_[m].getNormal().Z * wallObjects_[l].getNormal().X) *
418  Vec3D::dot(wallObjects_[n].getPosition(), wallObjects_[n].getNormal())
419  + (wallObjects_[n].getNormal().X * wallObjects_[l].getNormal().Z -
420  wallObjects_[n].getNormal().Z * wallObjects_[l].getNormal().X) *
421  Vec3D::dot(wallObjects_[m].getPosition(), wallObjects_[m].getNormal())
422  - (wallObjects_[n].getNormal().X * wallObjects_[m].getNormal().Z -
423  wallObjects_[m].getNormal().X * wallObjects_[n].getNormal().Z) *
424  Vec3D::dot(wallObjects_[l].getPosition(), wallObjects_[l].getNormal()),
425  +(wallObjects_[m].getNormal().X * wallObjects_[l].getNormal().Y -
426  wallObjects_[l].getNormal().X * wallObjects_[m].getNormal().Y) *
427  Vec3D::dot(wallObjects_[n].getPosition(), wallObjects_[n].getNormal())
428  - (wallObjects_[n].getNormal().X * wallObjects_[l].getNormal().Y -
429  wallObjects_[l].getNormal().X * wallObjects_[n].getNormal().Y) *
430  Vec3D::dot(wallObjects_[m].getPosition(), wallObjects_[m].getNormal())
431  + (wallObjects_[n].getNormal().X * wallObjects_[m].getNormal().Y -
432  wallObjects_[m].getNormal().X * wallObjects_[n].getNormal().Y) *
433  Vec3D::dot(wallObjects_[l].getPosition(), wallObjects_[l].getNormal())) * invdet;
434  }
435  }
436 
437  logger(VERBOSE, "%", *this);
438  for (const InfiniteWall& w : wallObjects_)
439  logger(VERBOSE, "wallObject %, %", w.getNormal(), w.getPosition());
440  for (Vec3D v : A_)
441  logger(VERBOSE, "A %", v);
442  for (Vec3D v : AB_)
443  logger(VERBOSE, "AB %", v);
444  for (Vec3D v : C_)
445  logger(VERBOSE, "C %", v);
446 }
447 
453 {
454  logger(WARN, "This function is deprecated, use IntersectionOfWalls::addObject(Vec3D, Vec3D) instead.");
455  addObject(normal, position * normal);
456 }
457 
466 void IntersectionOfWalls::createOpenPrism(std::vector<Vec3D> points, Vec3D prismAxis)
467 {
468  clear();
469  //note: use i+1 < points.size() instead of i < points.size()-1, otherwise it creates havoc if point has zero entries.
470  for (unsigned int i = 0; i + 1 < points.size(); i++)
471  addObject(Vec3D::cross(points[i] - points[i + 1], prismAxis), points[i]);
472 }
473 
481 void IntersectionOfWalls::createPrism(std::vector<Vec3D> points, Vec3D prismAxis)
482 {
483  createOpenPrism(points, prismAxis);
484  addObject(Vec3D::cross(points.back() - points.front(), prismAxis), points.front());
485 }
486 
494 void IntersectionOfWalls::createOpenPrism(std::vector<Vec3D> points)
495 {
496  Vec3D prismAxis = Vec3D::cross(
497  Vec3D::getUnitVector(points[1] - points[0]),
498  Vec3D::getUnitVector(points[2] - points[0]));
499  createOpenPrism(points, prismAxis);
500 }
501 
509 void IntersectionOfWalls::createPrism(std::vector<Vec3D> points)
510 {
511  Vec3D prismAxis = Vec3D::cross(
512  Vec3D::getUnitVector(points[1] - points[0]),
513  Vec3D::getUnitVector(points[2] - points[0]));
514  createPrism(points, prismAxis);
515 }
516 
531 bool IntersectionOfWalls::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
532 {
533  //transform coordinates into position-orientation frame
534  Vec3D position = p.getPosition() - getPosition();
535  getOrientation().rotateBack(position);
538  if (getDistanceAndNormal(position, p.getRadius() + s->getInteractionDistance(), distance, normal_return))
539  {
540  getOrientation().rotate(normal_return);
541  return true;
542  }
543  else
544  {
545  return false;
546  }
547 }
548 
566 bool IntersectionOfWalls::getDistanceAndNormal(const Vec3D& position, Mdouble wallInteractionRadius, Mdouble& distance,
567  Vec3D& normal_return) const
568 {
569  if (wallObjects_.empty())
570  {
571  logger(DEBUG, "Empty IntersectionOfWalls");
572  return false;
573  }
574 
575  distance = -1e20;
576  Mdouble distance2 = -1e20;
577  Mdouble distance3 = -1e20;
578  Mdouble distanceCurrent;
579  unsigned int id = 0;
580  unsigned int id2 = 0;
581  unsigned int id3 = 0;
582 
583  //For each wall we calculate the distance (distanceCurrent). To compute The walls with the minimal overlap
584  //The object has to touch each wall each wall (distanceCurrent) and keep the minimum distance (distance) and wall index (id
585  for (unsigned int i = 0; i < wallObjects_.size(); i++)
586  {
587  // Calculate distance to each wall (distanceCurrent);
588  distanceCurrent = wallObjects_[i].getDistance(position);
589  // The object has to touch each wall (distanceCurrent >= wallInteractionRadius), otherwise return false (i.e. no contact)
590  // This means that for each InfiniteWall in wallObjects_, the particle is either "inside"
591  // the wall or touching it. If not, there is no interaction.
592  if (distanceCurrent >= wallInteractionRadius)
593  {
594  return false;
595  }
596  // Find out which of the InfiniteWalls is interacting with the particle.
597  // Keep the minimum distance (distance) and wall index (id)
598  // and store up to two walls (id2, id3) and their distances (distance2, distance3),
599  // if the possible contact point is near the intersection between id and id2 (and id3)
600  if (distanceCurrent > distance)
601  {
602  if (distance > -wallInteractionRadius) //if distance was set previously
603  {
604  if (distance2 > -wallInteractionRadius) //if distance2 was set previously
605  {
606  distance3 = distance2;
607  id3 = id2;
608  }
609  distance2 = distance;
610  id2 = id;
611  }
612  distance = distanceCurrent;
613  id = i;
614  }
615  else if (distanceCurrent < -wallInteractionRadius)
616  {
617  continue;
618  }
619  else if (distanceCurrent > distance2)
620  {
621  if (distance2 > -wallInteractionRadius) //if distance2 was set previously
622  {
623  distance3 = distance2;
624  id3 = id2;
625  }
626  distance2 = distanceCurrent;
627  id2 = i;
628  }
629  else if (distanceCurrent > distance3)
630  {
631  distance3 = distanceCurrent;
632  id3 = i;
633  }
634  }
635 
636  //If we are here, the closest wall is id;
637  //if distance2>-P.Radius (and distance3>-P.Radius), the possible contact point
638  // is near the intersection between id and id2 (and id3)
639  if (distance2 > -wallInteractionRadius)
640  {
641  //D is the point on wall id closest to P
642  Vec3D D = position + wallObjects_[id].getNormal() * distance;
643  //If the distance of D to id2 is positive, the contact is with the intersection
644  bool intersection_with_id2 = (wallObjects_[id2].getDistance(D) > 0.0);
645 
646  if (distance3 > -wallInteractionRadius && (wallObjects_[id3].getDistance(D) > 0.0))
647  {
648  //E is the point on wall id2 closest to P
649  Vec3D E = position + wallObjects_[id2].getNormal() * distance2;
650  //If the distance of D to id2 is positive, the contact is with the intersection
651  bool intersection_with_id3 = (wallObjects_[id3].getDistance(E) > 0.0);
652 
653  if (intersection_with_id2)
654  {
655  if (intersection_with_id3)
656  {
657  //possible contact is with intersection of id,id2,id3
658  //we know id2<id3
659  if (id2 > id3)
660  {
661  auto id = id2;
662  id2 = id3;
663  id3 = id;
664  }
665  unsigned int index =
666  (id < id2) ? ((id3 - 2) * (id3 - 1) * id3 / 6 + (id2 - 1) * id2 / 2 + id) :
667  (id < id3) ? ((id3 - 2) * (id3 - 1) * id3 / 6 + (id - 1) * id / 2 + id2) :
668  ((id - 2) * (id - 1) * id / 6 + (id3 - 1) * id3 / 2 + id2);
669  normal_return = position - C_[index];
670  distance = sqrt(normal_return.getLengthSquared());
671  if (distance <= wallInteractionRadius) //note what if nan?
672  {
673  normal_return /= -distance;
674  return true; //contact with id,id2,id3
675  }
676  else
677  {
678  if (distance == distance)
679  return false;
680  }
681  }
682  }
683  else
684  {
685  intersection_with_id2 = true;
686  distance2 = distance3;
687  id2 = id3;
688  }
689  }
690 
691  if (intersection_with_id2)
692  { //possible contact is with intersection of id,id2
693  unsigned int index = (id > id2) ? ((id - 1) * id / 2 + id2) : ((id2 - 1) * id2 / 2 + id);
694  Vec3D AC = position - A_[index];
695  normal_return = AC - AB_[index] * Vec3D::dot(AC, AB_[index]);
696  distance = sqrt(normal_return.getLengthSquared());
697  if (distance <= wallInteractionRadius) //note what if nan?
698  {
699  normal_return /= -distance;
700  return true; //contact with id,id2,id3
701  }
702  else
703  {
704  if (distance == distance)
705  return false;
706  }
707  }
708  }
709  //contact is with id
710  normal_return = wallObjects_[id].getNormal();
711  return true;
712 }
713 
715 // * \param[in] move A reference to a Vec3D that denotes the direction and length
716 // * it should be moved with.
717 // * \details A function that moves the InterSectionOfWalls in a certain direction
718 // * by both moving the walls and all intersections. Note that the directions of the
719 // * intersections are not moved since they don't change when moving the IntersectionOfWalls
720 // * as a whole.
721 // * \todo We should use the position_ and orientation_ of the IntersectionOfWalls;
722 // * that way, IntersectionOfWalls can be moved with the standard BaseInteractable::move function,
723 // * getting rid of an anomaly in the code and removing the virtual from the move function. \author weinhartt
724 // */
725 //void IntersectionOfWalls::move(const Vec3D& move)
726 //{
727 // BaseInteractable::move(move);
728 // for (Vec3D& a : A_)
729 // {
730 // a += move;
731 // }
732 // for (Vec3D& c : C_)
733 // {
734 // c += move;
735 // }
736 // for (InfiniteWall& o : wallObjects_)
737 // {
738 // o.move(move);
739 // }
740 //}
741 
745 void IntersectionOfWalls::read(std::istream& is)
746 {
747  BaseWall::read(is);
748  std::string dummy;
749  int n;
750  is >> dummy >> n;
751 
752  Vec3D normal;
753  Vec3D position;
754  for (int i = 0; i < n; i++)
755  {
756  is >> dummy;
757  if (dummy != "normal")
758  {
759  Quaternion orientation;
760  is >> position >> dummy >> orientation;
761  addObject(orientation, position);
762  }
763  else
764  {
765  is >> normal >> dummy >> position;
766  addObject(normal, position);
767  }
768  }
769 }
770 
775 void IntersectionOfWalls::write(std::ostream& os) const
776 {
777  BaseWall::write(os);
778  os << " numIntersectionOfWalls " << wallObjects_.size();
779  for (const auto& wallObject : wallObjects_)
780  {
781  os << " position " << wallObject.getPosition() << " orientation " << wallObject.getOrientation();
782  }
783 }
784 
788 std::string IntersectionOfWalls::getName() const
789 {
790  return "IntersectionOfWalls";
791 }
792 
794 {
795  Vec3D max = getHandler()->getDPMBase()->getMax()-getPosition();
796  Vec3D min = getHandler()->getDPMBase()->getMin()-getPosition();
797  for (auto wall = wallObjects_.begin(); wall != wallObjects_.end(); wall++)
798  {
799  std::vector<Vec3D> points;
800  wall->createVTK(points, min, max);
801  for (auto other = wallObjects_.begin(); other != wallObjects_.end(); other++)
802  {
803  if (other != wall)
804  {
805  intersectVTK(points, -other->getNormal(), other->getPosition());
806  }
807  }
808  //rotate into real frame
809  for (auto& p : points)
810  {
811  getOrientation().rotate(p);
812  p += getPosition();
813  }
814  wall->addToVTK(points, vtk);
815  }
816 }
817 
818 //void IntersectionOfWalls::rotate(const Vec3D& rotate)
819 //{
820 // size_t n = 0;
821 // for (auto& wall : wallObjects_)
822 // {
823 // Vec3D p = wall.getPosition() - getPosition();
824 // wall.getOrientation().rotateBack(p);
825 // wall.setOrientation(wall.getOrientation().updateAngularDisplacement(rotate));
826 // wall.getOrientation().rotate(p);
827 // wall.setPosition(p + getPosition());
828 // setPointsAndLines(n);
829 // n++;
830 // }
831 //}
832 
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
void add3PointObject(Vec3D PointA, Vec3D PointB, Vec3D PointC)
void addTetraSTL(Vec3D PointA, Vec3D PointB, Vec3D PointC, Vec3D WallNormal, Mdouble Thickness, int wallidentifier)
constructs a tetrahedron for an STL file input
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
BaseSpecies is the class from which all other species are derived.
Definition: BaseSpecies.h:49
static Vec3D getUnitVector(const Vec3D &a)
Returns a unit Vec3D based on a.
Definition: Vector.cc:345
A IntersectionOfWalls is convex polygon defined as an intersection of InfiniteWall's.
void writeVTK(VTKContainer &vtk) const override
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
void intersectVTK(std::vector< Vec3D > &points, Vec3D normal, Vec3D position) const
Definition: BaseWall.cc:243
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
void read(std::istream &is) override
Move the IntersectionOfWalls to a new position, which is a Vec3D from the old position.
Vec3D getMin() const
Definition: DPMBase.h:623
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
virtual void setHandler(WallHandler *handler)
A function which sets the WallHandler for this BaseWall.
Definition: BaseWall.cc:127
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:102
bool getDistanceAndNormal(const BaseParticle &p, Mdouble &distance, Vec3D &normal_return) const override
Compute the distance from the wall for a given BaseParticle and return if there is a collision...
std::string getName() const override
Returns the name of the object, here the string "IntersectionOfWalls".
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
void addObject(Vec3D normal, Vec3D point)
Adds a wall to the set of infinite walls, given a normal vector pointing into the wall (i...
Mdouble getInteractionDistance() const
returns the largest separation distance at which adhesive short-range forces can occur.
Definition: BaseSpecies.h:146
~IntersectionOfWalls() override
Destructor.
std::vector< InfiniteWall > wallObjects_
The wall "segments"/directions that together make up the finite wall.
void createOpenPrism(std::vector< Vec3D > points, Vec3D prismAxis)
Creates an open prism which is a polygon between the points, except the first and last point...
void setSpecies(const ParticleSpecies *species)
sets species of subwalls as well
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
IntersectionOfWalls()
Default constructor.
std::vector< Vec3D > C_
A vector that stores the intersection point of three different InfiniteWall.
std::vector< Vec3D > A_
A vector that stores a point for each intersecting line between two different InfiniteWall.
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
Basic class for walls.
Definition: BaseWall.h:47
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1319
IntersectionOfWalls * copy() const override
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:80
Vec3D getMax() const
Definition: DPMBase.h:629
void clear()
Removes all parts of the walls.
Container to store all BaseWall.
Definition: WallHandler.h:42
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
void addPlate(const Vec3D &PointA, const Vec3D &PointB, const Vec3D &PointC, const Vec3D &WallNormal, const Mdouble &Thickness, int wallidentifier)
std::vector< Vec3D > AB_
A vector that stores the direction of the intersecting lines between two different InfiniteWall...
void set(Vec3D normal, Vec3D point)
Defines a standard wall, given an outward normal vector s.t. normal*x=normal*point for all x of the w...
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
void createPrism(std::vector< Vec3D > points, Vec3D prismAxis)
Creates an open prism which is a polygon between the points and extends infinitely in the PrismAxis d...
This is a class defining walls.
Definition: InfiniteWall.h:47
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
unsigned int getNumberOfObjects()
Returns the number of objects.
void addTetra(const Vec3D &PointA, const Vec3D &PointB, const Vec3D &PointC, Mdouble &Thickness)
constructs a tetrahedron from 3 input coordinates
Definition: Vector.h:49
IntersectionOfWalls & operator=(const IntersectionOfWalls &other)
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
void setPointsAndLines(unsigned int n)
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
void setHandler(WallHandler *wallHandler) override
A function which sets the WallHandler for this BaseWall.
void write(std::ostream &os) const override
Writes an IntersectionOfWalls to an output stream, for example a restart file.
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171