MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
InfiniteWall.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 
28 #include "InfiniteWall.h"
29 #include "Particles/BaseParticle.h"
30 #include "InteractionHandler.h"
31 #include "WallHandler.h"
32 #include "DPMBase.h"
33 
35 {
36  logger(DEBUG, "InfiniteWall::InfiniteWall ) finished");
37 }
38 
45  : BaseWall(w)
46 {
47  logger(DEBUG, "InfiniteWall::InfiniteWall(const InfiniteWall &p) finished");
48 }
49 
51  : BaseWall()
52 {
53  setSpecies(s);
54  logger(DEBUG, "InfiniteWall::InfiniteWall(const ParticleSpecies* s) finished");
55 }
56 
57 InfiniteWall::InfiniteWall(Vec3D normal, Vec3D point, const ParticleSpecies* species)
58 {
59  setNormal(normal);
60  setPosition(point);
61  setSpecies(species);
62 }
63 
72 InfiniteWall::InfiniteWall(Vec3D PointA, Vec3D PointB, Vec3D PointC, const ParticleSpecies* species)
73 {
74 //function returns infinite wall passing through 3 points
75 //subtract second and third from the first
76  Vec3D SubtB = PointA - PointB;
77  Vec3D SubtC = PointA - PointC;
78 
79  Vec3D WallNormal = Vec3D::cross(SubtB, SubtC);
80  //Check if walls coordinates inline, if true Do not build wall and give error message.
81 
82  if (WallNormal.getLengthSquared() == 0.0)
83  {
84  logger(ERROR,
85  "Error Building InfiniteWall out of 3 coordinates. Coordinates are in line, Wall not constructed.");
86  }
87  else
88  {
89  setNormal(WallNormal);
90  setPosition(PointA);
91  setSpecies(species);
92  }
93 
94 
95 }
96 
98 {
99  logger(DEBUG, "InfiniteWall::~InfiniteWall finished");
100 }
101 
106 {
107  return new InfiniteWall(*this);
108 }
109 
110 /*
111  * \details Defines a standard wall, given a normal vector pointing into the
112  * wall (i.e. out of the flow domain) and a point that the wall goes through,
113  * to give a plane defined by
114  * normal*x=normal*point.
115  * \param[in] normal A Vec3D that represents the normal to the wall.
116  * \param[in] point A Vec3D which is a point on the wall.
117  */
118 void InfiniteWall::set(Vec3D normal, Vec3D point)
119 {
120  setNormal(normal);
121  setPosition(point);
122 }
123 
127 void InfiniteWall::setNormal(const Vec3D normal)
128 {
129  setOrientationViaNormal(normal);
130 }
131 
139 void InfiniteWall::set(Vec3D normal, Mdouble positionInNormalDirection)
140 {
141  logger(WARN, "InfiniteWall::set(Vec3D, Mdouble) is deprecated. Use set(Vec3D, Vec3D) instead.");
142  set(normal, positionInNormalDirection * normal);
143 }
144 
150 {
151  return getOrientation().getDistance(otherPosition, getPosition());
152 }
153 
166 bool InfiniteWall::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
167 {
168  distance = getDistance(p.getPosition());
169  if (distance >= p.getWallInteractionRadius(this))
170  return false;
171  normal_return = getOrientation().getAxis();
172  return true;
173 }
174 
178 void InfiniteWall::read(std::istream& is)
179 {
180  BaseWall::read(is);
181  Vec3D normal;
182  if (helpers::readOptionalVariable(is,"normal",normal)) setNormal(normal);
183 }
184 
188 void InfiniteWall::oldRead(std::istream& is)
189 {
190  std::string dummy;
191  Vec3D velocity;
192  Vec3D position;
193  Vec3D normal;
194  is >> dummy >> normal >> dummy >> position >> dummy >> velocity;
195  setPosition(position);
196  setVelocity(velocity);
197  setOrientation(Quaternion(normal));
198 }
199 
203 std::string InfiniteWall::getName() const
204 {
205  return "InfiniteWall";
206 }
207 
212 {
213  return getOrientation().getAxis();
214 }
215 
222 void InfiniteWall::createVTK(std::vector<Vec3D>& myPoints) const
223 {
224  Vec3D max = getHandler()->getDPMBase()->getMax();
225  Vec3D min = getHandler()->getDPMBase()->getMin();
226  createVTK(myPoints, min, max);
227 }
228 
229 void InfiniteWall::createVTK(std::vector<Vec3D>& myPoints, const Vec3D min, const Vec3D max) const
230 {
231  const Vec3D& n = getOrientation().getAxis();
232  const Vec3D& p = getPosition();
233 
234  if (fabs(n.X) > 0.5)
235  {
236  // If the wall normal has a nonzero x-component,
237  // We first find four intersection points with the four domain edges pointing in x-direction.
238  // Because these points might be outside the domain in x-direction, we use the intersection function have points in the domain
239  myPoints.emplace_back(p.X - ((min.Y - p.Y) * n.Y + (min.Z - p.Z) * n.Z) / n.X, min.Y, min.Z);
240  myPoints.emplace_back(p.X - ((min.Y - p.Y) * n.Y + (max.Z - p.Z) * n.Z) / n.X, min.Y, max.Z);
241  myPoints.emplace_back(p.X - ((max.Y - p.Y) * n.Y + (max.Z - p.Z) * n.Z) / n.X, max.Y, max.Z);
242  myPoints.emplace_back(p.X - ((max.Y - p.Y) * n.Y + (min.Z - p.Z) * n.Z) / n.X, max.Y, min.Z);
243  intersectVTK(myPoints, Vec3D(+1, 0, 0), Vec3D(max.X, 0, 0));
244  intersectVTK(myPoints, Vec3D(-1, 0, 0), Vec3D(min.X, 0, 0));
245  }
246  else if (fabs(n.Y) > 0.5)
247  {
248  // Else, if the wall normal has a nonzero y-component ...
249  myPoints.emplace_back(min.X, p.Y - ((min.X - p.X) * n.X + (min.Z - p.Z) * n.Z) / n.Y, min.Z);
250  myPoints.emplace_back(min.X, p.Y - ((min.X - p.X) * n.X + (max.Z - p.Z) * n.Z) / n.Y, max.Z);
251  myPoints.emplace_back(max.X, p.Y - ((max.X - p.X) * n.X + (max.Z - p.Z) * n.Z) / n.Y, max.Z);
252  myPoints.emplace_back(max.X, p.Y - ((max.X - p.X) * n.X + (min.Z - p.Z) * n.Z) / n.Y, min.Z);
253  intersectVTK(myPoints, Vec3D(0, +1, 0), Vec3D(0, max.Y, 0));
254  intersectVTK(myPoints, Vec3D(0, -1, 0), Vec3D(0, min.Y, 0));
255  }
256  else
257  {
258  // Else, the wall normal has to have a a nonzero z-component ...
259  myPoints.emplace_back(min.X, min.Y, p.Z - ((min.Y - p.Y) * n.Y + (min.X - p.X) * n.X) / n.Z);
260  myPoints.emplace_back(min.X, max.Y, p.Z - ((max.Y - p.Y) * n.Y + (min.X - p.X) * n.X) / n.Z);
261  myPoints.emplace_back(max.X, max.Y, p.Z - ((max.Y - p.Y) * n.Y + (max.X - p.X) * n.X) / n.Z);
262  myPoints.emplace_back(max.X, min.Y, p.Z - ((min.Y - p.Y) * n.Y + (max.X - p.X) * n.X) / n.Z);
263  intersectVTK(myPoints, Vec3D(0, 0, +1), Vec3D(0, 0, max.Z));
264  intersectVTK(myPoints, Vec3D(0, 0, -1), Vec3D(0, 0, min.Z));
265  }
266 }
267 
269 {
270  std::vector<Vec3D> points;
271  createVTK(points);
272  addToVTK(points, vtk);
273 }
274 
275 
276 bool
278  Mdouble& overlap) const
279 {
280  //first check: if the bounding sphere does not touch the wall, there is no contact.
282  {
283  return false;
284  }
285  Vec3D normalBodyFixed = getOrientation().getAxis();
286  p.getOrientation().rotateBack(normalBodyFixed);
287  Vec3D xWallBodyFixed = getPosition() - p.getPosition();
288  p.getOrientation().rotateBack(xWallBodyFixed);
289  Vec3D axes = p.getAxes();
290  Mdouble eps1 = p.getExponentEps1();
291  Mdouble eps2 = p.getExponentEps2();
292 
293  Vec3D furthestPoint = getFurthestPointSuperQuadric(normalBodyFixed, axes, eps1, eps2);
294  overlap = Vec3D::dot(xWallBodyFixed - furthestPoint, -normalBodyFixed);
295  if (overlap > 0)
296  {
297  Vec3D overlapBody = overlap * normalBodyFixed;
298  Vec3D contactPoint = furthestPoint - overlapBody / 2;
299  p.getOrientation().rotate(contactPoint);
300  contactPoint += p.getPosition();
301  distance = (contactPoint - overlapBody / 2 - p.getPosition()).getLength();
302  normal_return = getOrientation().getAxis();
303  return true;
304  }
305  return false;
306 }
307 
308 
310 Vec3D InfiniteWall::getFurthestPointSuperQuadric(const Vec3D& normalBodyFixed, const Vec3D& axes, Mdouble eps1,
311  Mdouble eps2) const
312 {
313  Vec3D furthestPoint;
314  if (std::abs(normalBodyFixed.X) > 1e-10)
315  {
316  Mdouble alpha = std::abs(normalBodyFixed.Y * axes.Y / normalBodyFixed.X / axes.X);
317  Mdouble gamma = std::pow(1 + std::pow(alpha, 2 / eps2), eps2 / eps1 - 1);
318  Mdouble beta = std::pow(gamma * std::abs(normalBodyFixed.Z * axes.Z / normalBodyFixed.X / axes.X),
319  eps1 / (2 - eps1));
320  furthestPoint.X = axes.X * mathsFunc::sign(normalBodyFixed.X) /
321  (std::pow(std::pow(1 + std::pow(alpha, 2 / eps2), eps2 / eps1) + std::pow(beta, 2 / eps1),
322  eps1 / 2));
323  furthestPoint.Y = axes.Y * alpha / axes.X * std::abs(furthestPoint.X) * mathsFunc::sign(normalBodyFixed.Y);
324  furthestPoint.Z = axes.Z * beta / axes.X * std::abs(furthestPoint.X) * mathsFunc::sign(normalBodyFixed.Z);
325  }
326  else if (std::abs(normalBodyFixed.Y) > 1e-10)
327  {
328  Mdouble beta = std::pow(std::abs(normalBodyFixed.Z * axes.Z / normalBodyFixed.Y / axes.Y), eps1 / (2 - eps1));
329  furthestPoint.Y = axes.Y / std::pow((1 + std::pow(beta, 2/eps1)), eps1 / 2) * mathsFunc::sign(normalBodyFixed.Y);
330  furthestPoint.Z = axes.Z / axes.Y * std::abs(furthestPoint.Y) * beta * mathsFunc::sign(normalBodyFixed.Z);
331  }
332  else
333  {
334  furthestPoint.Z = axes.Z * mathsFunc::sign(normalBodyFixed.Z);
335  }
336  return furthestPoint;
337 }
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
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...
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
std::string getName() const override
Writes the InfiniteWall to an output stream, usually a restart file.
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
void createVTK(std::vector< Vec3D > &myPoints) const
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
Vec3D getMin() const
Definition: DPMBase.h:623
Vec3D getNormal() const
Access function for normal.
static void addToVTK(const std::vector< Vec3D > &points, VTKContainer &vtk)
Takes the points provided and adds a triangle strip connecting these points to the vtk container...
Definition: BaseWall.cc:471
int sign(T val)
This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0...
Definition: ExtendedMath.h:95
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
bool getDistanceNormalOverlapSuperquadric(const SuperQuadricParticle &p, Mdouble &distance, Vec3D &normal_return, Mdouble &overlap) const override
Compute the distance from the wall for a given BaseParticle and return if there is a collision...
void oldRead(std::istream &is)
Reads InfiniteWall from an old-style restart file.
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
void writeVTK(VTKContainer &vtk) const override
Vec3D getAxis() const
Converts the quaternions into a normal vector by rotating the vector x=(1,0,0); see See Wiki for deta...
Definition: Quaternion.cc:501
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
static Mdouble getDistance(const Quaternion &a, const Quaternion &b)
Calculates the distance between two Quaternion: .
Definition: Quaternion.cc:188
InfiniteWall()
Default constructor, the normal is infinitely long.
Definition: InfiniteWall.cc:34
~InfiniteWall() override
Default destructor.
Definition: InfiniteWall.cc:97
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
InfiniteWall * copy() const override
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
Basic class for walls.
Definition: BaseWall.h:47
bool readOptionalVariable(std::istream &is, const std::string &name, T &variable)
Reads optional variables in the restart file.
Definition: Helpers.h:247
Mdouble getExponentEps1() const override
Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter ...
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
Mdouble Y
Definition: Vector.h:65
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Vec3D getFurthestPointSuperQuadric(const Vec3D &normalBodyFixed, const Vec3D &axes, Mdouble eps1, Mdouble eps2) const override
Largely untested, use at your own risk for anything other than ellipsoids.
void read(std::istream &is) override
Reads InfiniteWall from a restart file.
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...
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
This is a class defining walls.
Definition: InfiniteWall.h:47
Mdouble getDistance(Vec3D position) const
Returns the distance of the wall to the particle.
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
Definition: Vector.h:49
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
Mdouble getExponentEps2() const override
Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter...
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
void setNormal(Vec3D normal)
Changes the normal of the InfiniteWall.
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171