MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BasicIntersectionOfWalls.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 <limits>
27 
29 #include "Particles/BaseParticle.h"
31 #include "InteractionHandler.h"
32 #include "WallHandler.h"
33 #include "DPMBase.h"
34 
36 {
37  logger(DEBUG, "BasicIntersectionOfWalls::BasicIntersectionOfWalls ) finished");
38 }
39 
46  : BaseWall(b)
47 {
48  for (auto& w : b.walls_)
49  {
50  walls_.push_back(w->copy());
51  }
52  logger(DEBUG, "BasicIntersectionOfWalls::BasicIntersectionOfWalls(const BasicIntersectionOfWalls &p) finished");
53 }
54 
56 {
57  for (auto& w : walls_)
58  {
59  delete w;
60  }
61  logger(DEBUG, "BasicIntersectionOfWalls::~BasicIntersectionOfWalls finished");
62 }
63 
68 {
69  return new BasicIntersectionOfWalls(*this);
70 }
71 
78 {
79  return walls_.size();
80 }
81 
82 
83 /*
84  * \param[in] normal A Vec3D that represents the normal to the wall.
85  * \param[in] point A Vec3D which is a point on the wall.
86  * \details Sets the wall such that for all points x on the wall it holds that
87  * normal*x=normal*point.
88  */
91 {
92  walls_.push_back(wall.copy());
93  walls_.back()->setId(walls_.size());
94 }
95 
104 {
105  if (walls_.empty())
106  {
107  logger(DEBUG, "Empty BasicIntersectionOfWalls");
108  return false;
109  }
110 
111  distance = -constants::inf; //distance of the closest wall
112  Mdouble distance2 = -constants::inf; //distance of the second-closest wall
113  Mdouble distance3 = -constants::inf; //distance of the third-closest wall
114  Mdouble distanceCurrent; //distance of teh current wall
115  Vec3D normal2; //normal of the second-closest wall
116  Vec3D normal3; //normal of the third-closest wall
117  Vec3D normalCurrent; //normal of the current wall
118  unsigned int id = static_cast<unsigned int>(-1); //id of the closest wall
119  unsigned int id2 = static_cast<unsigned int>(-1); //id of the second-closest wall
120  unsigned int id3 = static_cast<unsigned int>(-1); //id of the third-closest wall
121  Mdouble wallInteractionRadius = p.getWallInteractionRadius(this);
122  Vec3D position = p.getPosition() - getPosition();
123  getOrientation().rotateBack(position);
124  SphericalParticle shifted;
125  shifted.setSpecies(p.getSpecies());
126  shifted.setPosition(position);
127  shifted.setRadius(p.getRadius());
128 
129  //The object has to touch each wall each wall (distanceCurrent) and keep the minimum distance (distance) and wall index (id)
130  for (unsigned int i = 0; i < walls_.size(); i++)
131  {
132  //immediately stop if one of teh walls is not in interaction distance
133  if (!walls_[i]->getDistanceAndNormal(shifted, distanceCurrent, normalCurrent))
134  return false;
135  // Find out which of the walls is interacting with the particle.
136  // Keep the minimum distance (distance) and wall index (id)
137  // and store up to two walls (id2, id3) and their distances (distance2, distance3),
138  // if the possible contact point is near the intersection between id and id2 (and id3)
139  if (distanceCurrent > distance)
140  {
141  if (distance > -wallInteractionRadius)
142  {
143  if (distance2 > -wallInteractionRadius)
144  {
145  distance3 = distance;
146  normal3 = normal;
147  id3 = id;
148  }
149  else
150  {
151  distance2 = distance;
152  normal2 = normal;
153  id2 = id;
154  }
155  }
156  distance = distanceCurrent;
157  normal = normalCurrent;
158  id = i;
159  }
160  else if (distanceCurrent > -wallInteractionRadius)
161  {
162  if (distance2 > -wallInteractionRadius)
163  {
164  distance3 = distanceCurrent;
165  normal3 = normalCurrent;
166  id3 = i;
167  }
168  else
169  {
170  distance2 = distanceCurrent;
171  normal2 = normalCurrent;
172  id2 = i;
173  }
174  }
175  }
176  //logger(INFO, "particle %: contact with % % %, n % % %", p.getId(), distance, distance2, distance3, normal, normal2, normal3);
177 
178 
179  //If we are here, the closest wall is id;
180  //if distance2>-P.Radius, it's a possible face or vertex contact
181  //if distance3>-P.Radius, it's a possible vertex contact
182  if (distance2 > -wallInteractionRadius)
183  {
184  //contact point on wall id
185  SphericalParticle contact;
186  contact.setSpecies(p.getSpecies());
187  contact.setPosition(position + distance * normal);
188  contact.setRadius(0);
189 
190  //If the distance of D to id2 is positive, the contact is with the intersection
191  bool contactPointOutsideID2 = !walls_[id2]->getDistanceAndNormal(contact, distanceCurrent, normalCurrent);
192 
193  if (distance3 > -wallInteractionRadius &&
194  !walls_[id3]->getDistanceAndNormal(contact, distanceCurrent, normalCurrent))
195  {
196  if (contactPointOutsideID2)
197  {
198  //possible contact is with intersection of id,id2,id3
199  //we know id2<id3
200  unsigned int index =
201  (id < id2) ? ((id3 - 2) * (id3 - 1) * id3 / 6 + (id2 - 1) * id2 / 2 + id) :
202  (id < id3) ? ((id3 - 2) * (id3 - 1) * id3 / 6 + (id - 1) * id / 2 + id2) :
203  ((id - 2) * (id - 1) * id / 6 + (id3 - 1) * id3 / 2 + id2);
204  //find vertex C
205  Matrix3D N(normal.X, normal.Y, normal.Z, normal2.X, normal2.Y, normal2.Z, normal3.X, normal3.Y,
206  normal3.Z);
207  Vec3D P(Vec3D::dot(position, normal) + distance, Vec3D::dot(position, normal2) + distance2,
208  Vec3D::dot(position, normal3) + distance3);
209  Vec3D vertex = N.ldivide(P);
210  normal = position - vertex;
211  distance = sqrt(normal.getLengthSquared());
212  if (distance >= wallInteractionRadius)
213  return false; //no contact
214  normal /= -distance;
215  //logger(INFO, "vertex contact with % % %: pos %, n %, d %", id, id2, id3, p.getPosition(), normal, distance);
216  getOrientation().rotate(normal);
217  return true;
218  }
219  else
220  {
221  contactPointOutsideID2 = true;
222  distance2 = distance3;
223  id2 = id3;
224  }
225  }
226 
227  if (contactPointOutsideID2)
228  {
229  //find contact point C
230  normal3 = Vec3D::cross(normal, normal2);
231  Matrix3D N(normal.X, normal.Y, normal.Z, normal2.X, normal2.Y, normal2.Z, normal3.X, normal3.Y,
232  normal3.Z);
233  Vec3D P(Vec3D::dot(position, normal) + distance, Vec3D::dot(position, normal2) + distance2,
234  Vec3D::dot(position, normal3));
235  Vec3D C = N.ldivide(P);
236  normal = position - C;
237  distance = sqrt(normal.getLengthSquared());
238  if (distance >= wallInteractionRadius)
239  return false; //no contact
240  normal /= -distance;
241  //logger(INFO, "edge contact with % %: normal %", id, id2, normal);
242  getOrientation().rotate(normal);
243  return true;
244  }
245  }
246  //logger(INFO, "face contact with %", id);
247  getOrientation().rotate(normal);
248  return true;
249 }
250 
254 void BasicIntersectionOfWalls::read(std::istream& is)
255 {
256  BaseWall::read(is);
257  std::string dummy;
258  unsigned size;
259  // read e.g. "numIntersectionOfWalls 2"
260  is >> dummy >> size;
261  for (unsigned i = 0; i < size; i++)
262  {
263  std::string type;
264  // read e.g. "IntersectionOfWalls"
265  is >> type;
266  BaseWall* wall = getHandler()->createObject(type);
267  wall->setHandler(getHandler());
268  wall->read(is);
269  walls_.push_back(wall);
270  walls_.back()->setId(walls_.size());
271  }
272 }
273 
277 void BasicIntersectionOfWalls::write(std::ostream& os) const
278 {
279  BaseWall::write(os);
280  os << " numIntersectionOfWalls " << walls_.size();
281  for (auto w : walls_)
282  {
283  os << " ";
284  w->write(os);
285  }
286 }
287 
292 {
293  return "BasicIntersectionOfWalls";
294 }
295 
296 void BasicIntersectionOfWalls::getVTK(std::vector<Vec3D>& points, std::vector<std::vector<double>>& triangleStrips)
297 {
298 // for (auto w : walls_) {
299 // w->writeVTK (points, triangleStrips);
300 // }
301  //writes points and strips for all walls; points are added to the global point vector, but the strips are held back
302  std::vector<std::vector<double>> myTriangleStrips;
303  unsigned long n = points.size();
304  for (auto w : walls_)
305  {
306  w->getVTK(points, myTriangleStrips);
307  }
308  //add position of the BasicIntersectionOfWalls to the point
309  for (std::vector<Vec3D>::iterator p = points.begin() + n; p != points.end(); p++)
310  {
311  getOrientation().rotate(*p);
312  }
313  //create a vector which points are in the wall (actually, only the new points are necessary)
314  std::vector<bool> pointInWall;
315  pointInWall.reserve(points.size());
316  SphericalParticle particle;
317  particle.setSpecies(getSpecies());
318  particle.setRadius(1e-10); //points within that distance are declared part of the wall
319  Mdouble distance;
320  Vec3D normal;
321  for (auto p : points)
322  {
323  particle.setPosition(p);
324  pointInWall.push_back(getDistanceAndNormal(particle, distance, normal));
325  }
326  //now loop through myTriangleStrips to find the strip parts that are fully inside the wall
327  std::vector<double> strip;
328  for (auto t : myTriangleStrips)
329  {
330  for (unsigned i : t)
331  {
332  if (pointInWall[i] == true)
333  {
334  strip.push_back(i);
335  }
336  else
337  {
338  if (strip.size() > 2)
339  triangleStrips.push_back(strip);
340  strip.clear();
341  }
342  }
343  if (strip.size() > 2)
344  triangleStrips.push_back(strip);
345  strip.clear();
346  }
348 }
349 
351 {
352  logger.assert_always(walls_.size() > i, "Index % exceeds number of walls %", i, walls_.size());
353  return walls_[i];
354 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
A basic particle.
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
BasicIntersectionOfWalls()
Default constructor, the normal is infinitely long.
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
This is a class defining walls.
unsigned long getNumberOfObjects()
Returns the number of objects.
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
void setSpecies(const ParticleSpecies *species)
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
const Mdouble inf
Definition: GeneralDefine.h:44
std::vector< BaseWall * > walls_
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
static BaseWall * createObject(const std::string &type)
Create a new wall, with the type given as a string (required for restarting).
Definition: WallHandler.cc:122
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
void write(std::ostream &os) const override
Writes the BasicIntersectionOfWalls to an output stream, usually a restart file.
BaseWall * getObject(unsigned i)
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
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:80
void read(std::istream &is) override
Reads BasicIntersectionOfWalls from a restart file.
Mdouble Y
Definition: Vector.h:65
Vec3D ldivide(const Vec3D &b)
A.ldivide(b) computes the solution x to A*x=b.
Definition: Matrix.cc:373
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
void add(BaseWall &wall)
Defines a standard wall, given an outward normal vector s.t. normal*x=normal*point for all x of the w...
BasicIntersectionOfWalls * copy() const override
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
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...
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
~BasicIntersectionOfWalls() override
Default destructor.
void getVTK(std::vector< Vec3D > &points, std::vector< std::vector< double >> &triangleStrips)
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
std::string getName() const override
Returns the name of the object, in this case the string "BasicIntersectionOfWalls".
Mdouble Z
Definition: Vector.h:65
virtual BaseWall * copy() const =0
Pure virtual function that can be overwritten in inherited classes in order to copy a BaseWall...