MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ScrewsymmetricIntersectionOfWalls.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 
27 #include "Particles/BaseParticle.h"
28 #include "WallHandler.h"
29 #include "DPMBase.h"
30 using mathsFunc::square;
31 
33 {
34  logger(ERROR, "ScrewsymmetricIntersectionOfWalls() is not yet completely implemented");
35  logger(DEBUG, "ScrewsymmetricIntersectionOfWalls() finished");
36 }
37 
42  : IntersectionOfWalls(other)
43 {
44  logger(ERROR, "ScrewsymmetricIntersectionOfWalls() is not yet completely implemented");
45  logger(DEBUG, "ScrewsymmetricIntersectionOfWalls(const ScrewsymmetricIntersectionOfWalls &p) finished");
46 }
47 
49  std::vector<ScrewsymmetricIntersectionOfWalls::normalAndPosition> walls,
50  const ParticleSpecies* species)
51  : IntersectionOfWalls(walls, species)
52 {
53  logger(ERROR, "ScrewsymmetricIntersectionOfWalls() is not yet completely implemented");
54  setPosition(position);
55  setOrientationViaNormal(orientation);
56 }
57 
59 {
60  logger(DEBUG, "~ScrewsymmetricIntersectionOfWalls() finished.");
61 }
62 
68 {
69  if (this == &other)
70  {
71  return *this;
72  }
73  else
74  {
75  return *(other.copy());
76  }
77 }
78 
83 {
84  return new ScrewsymmetricIntersectionOfWalls(*this);
85 }
86 
98 {
99  // squared radial position of the particle
100  const Vec3D positionLabFrame = p.getPosition() - getPosition();
101  const Mdouble rho2 = positionLabFrame.X * positionLabFrame.X + positionLabFrame.Y * positionLabFrame.Y;
102 
103  // should thickness be added here?
104  const Mdouble wallInteractionRadius = p.getWallInteractionRadius(this) + thickness_;
105 
106  // if the particle is outside the cylinder containing the screw there is no collision
107  if (//rho2 > square(rMax_ + wallInteractionRadius) ||
108  positionLabFrame.Z > length_ + p.getWallInteractionRadius(this) ||
109  positionLabFrame.Z < -p.getWallInteractionRadius(this))
110  {
111  return false;
112  }
113  Vec3D normalVector;
114  Vec3D radialVector;
115  Mdouble deltaN;
116  // positionLabFrame = normalVector + radialVector + axialVector
117  computeNormalRadialDeltaN(positionLabFrame, normalVector, radialVector, deltaN);
118 
119  // radial position of the particle
120  const Mdouble r = std::sqrt(rho2);
121 
122  //determine wall distance, normal and contact in axissymmetric coordinates
123  //and transform from axisymmetric coordinates
124  Vec3D normalRN;
125  if (!IntersectionOfWalls::getDistanceAndNormal(Vec3D(r, 0, deltaN), p.getWallInteractionRadius(this), distance, normalRN))
126  {
127  //if not in contact
128  return false;
129  }
130  else
131  {
132  //if in contact
133  if (r != 0)
134  radialVector /= r;
135  else //in this case the tangential vector is irrelevant
136  logger(WARN, "Warning: Particle % is exactly on the symmetry axis of wall %", p.getIndex(), getIndex());
137  //normal_Return = normalVector.Z * deltaN + normalVector.X * r;
138  normal_return = normalRN.Z * normalVector + normalRN.X * radialVector;
139  normal_return.normalise();
140  return true;
141  }
142 }
143 
153  Vec3D& radialVector, Mdouble& deltaN) const
154 {
155  // radial position of the particle
156  const Mdouble rho = std::sqrt(positionLabFrame.X * positionLabFrame.X
157  + positionLabFrame.Y * positionLabFrame.Y);
158 
159  // The rescaled length of the screw (length_/(2*pi*numberOfTurns_)).
160  const Mdouble h = 0.5 * pitch_ / constants::pi;
161 
162  // The pitch length of the screw (length_/numberOfTurns_).
163  const Mdouble deltaZ = computeDeltaZ(positionLabFrame, h, pitch_);
164 
165 
166  // trigonometric functions relative to the particle angle
167  const Mdouble cosXi = positionLabFrame.X / rho;
168  const Mdouble sinXi = positionLabFrame.Y / rho;
169  if (rightHandedness_)
170  {
171  normalVector.X = h * sinXi;
172  normalVector.Y = -h * cosXi;
173  normalVector.Z = rho;
174  }
175  else
176  {
177  normalVector.X = -h * cosXi;
178  normalVector.Y = h * sinXi;
179  normalVector.Z = rho;
180  }
181 
182  // takes the right direction (+/-) of the vector and normalizes
183  normalVector *= -deltaZ;
184  normalVector /= normalVector.getLength();
185 
186  // radial vector at the particle position
187  radialVector.X = cosXi;
188  radialVector.Y = sinXi;
189  radialVector.Z = 0.0;
190 
191  // The half-thickness of the screw.
192  const Mdouble delta = 0.5 * thickness_;
193 
194  // cosine of the helix angle at the particle position
195  const Mdouble cosEta = 1.0 / std::sqrt(1.0 + (h * h / rho / rho));
196 
197  // component of the particle-blade_center distance normal to the blade surface
198  deltaN = fabs(deltaZ) * cosEta - delta;
199 }
200 
210 {
211  // oriented axial distance between the particle's centre and the blade centre
212  Mdouble deltaZ;
213 
214  // normal to the blade at the particle position
215  if (rightHandedness_) // right-handed thread
216  {
217  // angular coordinate of the particle
218  // IMPORTANT: this angle needs to be defined in the interval [0, +2*pi[ radians!
219  Mdouble xi = atan2(positionLabFrame.Y, positionLabFrame.X);
220  if (xi < 0.0)
221  xi += 2.0 * constants::pi;
222 
223  deltaZ = fmod(positionLabFrame.Z - h * (xi + getOrientation().getAxis().Z) -
224  static_cast<int> (positionLabFrame.Z / pitch), pitch);
225  logger(DEBUG, "xi: %, deltaZ: %", xi, deltaZ);
226  }
227  else // left-handed thread
228  {
229  // angular coordinate of the particle
230  // IMPORTANT: this angle needs to be defined in the interval [0, +2*pi[ radians!
231  Mdouble xi = atan2(-positionLabFrame.Y, positionLabFrame.X);
232  if (xi < 0.0)
233  xi += 2.0 * constants::pi;
234  xi += 0.5 * constants::pi;
235  xi = fmod(xi, 2.0 * constants::pi);
236 
237  deltaZ = fmod(positionLabFrame.Z - h * (xi + 0.5 * constants::pi - getOrientation().getAxis().Z) -
238  static_cast<int> (positionLabFrame.Z / pitch), pitch);
239  logger(DEBUG, "xi: %, deltaZ: %", xi, deltaZ);
240  }
241 
242  if (deltaZ > 0.5 * pitch)
243  {
244  deltaZ -= pitch;
245  }
246  else if (deltaZ < -0.5 * pitch)
247  {
248  deltaZ += pitch;
249  }
250  return deltaZ;
251 }
252 
258 {
260 }
261 
266 void ScrewsymmetricIntersectionOfWalls::write(std::ostream& os) const
267 {
269 }
270 
275 {
276  return "ScrewsymmetricIntersectionOfWalls";
277 }
278 
280 {
282 }
283 
285 {
287  Vec3D rMin = min - getPosition();
288  q.rotateBack(rMin); //set min/max initial values to values of first corner point
289  Vec3D rMax = max - getPosition();
290  q.rotateBack(rMax); //set min/max initial values to values of first corner point
291 
292  Mdouble r = std::sqrt(std::max(rMax.Y * rMax.Y + rMax.Z * rMax.Z, rMin.Y * rMin.Y + rMin.Z * rMin.Z));
293  max = Vec3D(r, 0.001, std::max(rMin.X,rMax.X));
294  min = Vec3D(0, 0, std::min(rMin.X,rMax.X));
295  //std::cout << "r=" << r << std::endl;
296 }
297 
299 {
300  for (auto wall = wallObjects_.begin(); wall != wallObjects_.end(); wall++)
301  {
302  //plot each of the intersecting walls
303  std::vector<Vec3D> myPoints;
304 
305  //first create a slice of non-rotated wall in the xz plane, 0<y<1
306  Vec3D min = getHandler()->getDPMBase()->getMin();
307  Vec3D max = getHandler()->getDPMBase()->getMax();
308  convertLimits(min, max);
309 
310  //create the basic slice for the first wall using the InfiniteWall routine
311  wall->createVTK(myPoints, min, max);
312 
313  //create intersections with the other walls, similar to the IntersectionOfWalls routine
314  for (auto other = wallObjects_.begin(); other != wallObjects_.end(); other++)
315  {
316  if (other != wall)
317  {
318  intersectVTK(myPoints, -other->getNormal(), other->getPosition());
319  }
320  }
321 
322  //only keep the y=0 values
323  std::vector<Vec3D> rzVec;
324  for (auto& p: myPoints)
325  {
326  if (p.Y == 0)
327  {
328  rzVec.push_back(p);
329  }
330  }
331  if (rzVec.empty())
332  return;
333 
334  //create points on the unit circle
335  unsigned nr = 180;
336  struct XY
337  {
338  double X;
339  double Y;
340  };
341  std::vector<XY> xyVec;
342  for (unsigned ir = 0; ir < nr; ir++)
343  {
344  Mdouble angle = 2.0 * constants::pi * ir / nr;
345  xyVec.push_back({cos(angle), sin(angle)});
346  }
347 
348  //now create rings of points on the axisym. shape
350  unsigned long nPoints = vtk.points.size();
351  Vec3D p;
352  //Vec3D o = getOrientation().getAxis();
353  for (auto rz : rzVec)
354  {
355  for (auto xy : xyVec)
356  {
357  p = Vec3D(rz.Z, rz.X * xy.X, rz.X * xy.Y);
358  getOrientation().rotate(p);
359  p += getPosition();
360  vtk.points.push_back(p);
361  }
362  }
363 
364  //finally create the connectivity matri to plot shell-like triangle strips.
365  unsigned long nz = rzVec.size();
366  unsigned long nCells = vtk.triangleStrips.size();
367  for (unsigned iz = 0; iz < nz - 1; iz++)
368  {
369  std::vector<double> cell;
370  cell.reserve(2 * nr + 2);
371  for (unsigned ir = 0; ir < nr; ir++)
372  {
373  cell.push_back(nPoints + ir + iz * nr);
374  cell.push_back(nPoints + ir + (iz + 1) * nr);
375  }
376  cell.push_back(nPoints + iz * nr);
377  cell.push_back(nPoints + (iz + 1) * nr);
378  vtk.triangleStrips.push_back(cell);
379  }
380  }
381 
382 }
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
unsigned int getIndex() const
Returns the index of the object in the handler.
Definition: BaseObject.h:118
const Vec3D getAxis() const
Definition: BaseWall.cc:503
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
A IntersectionOfWalls is convex polygon defined as an intersection of InfiniteWall's.
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 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.
bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const final
Computes the distance from the wall for a given BaseParticle and returns true if there is a collision...
Vec3D getMin() const
Definition: DPMBase.h:623
Mdouble pitch_
axial length of one rotation
ScrewsymmetricIntersectionOfWalls * copy() const final
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
void read(std::istream &is) final
reads wall
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 write(std::ostream &os) const final
outputs wall
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
std::vector< InfiniteWall > wallObjects_
The wall "segments"/directions that together make up the finite wall.
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
void convertLimits(Vec3D &min, Vec3D &max) const
Use ScrewsymmetricIntersectionOfWalls to define screwsymmetric walls, such as cylinders, cones, etc.
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
ScrewsymmetricIntersectionOfWalls & operator=(const ScrewsymmetricIntersectionOfWalls &other)
Copy assignment operator.
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Mdouble computeDeltaZ(const Vec3D &positionLabFrame, Mdouble h, Mdouble pitch) const
void computeNormalRadialDeltaN(const Vec3D &positionLabFrame, Vec3D &normalVector, Vec3D &radialVector, Mdouble &deltaN) const
Vec3D getMax() const
Definition: DPMBase.h:629
Mdouble Y
Definition: Vector.h:65
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
std::vector< Vec3D > points
Definition: BaseWall.h:38
Mdouble thickness_
The thickness of the screw blade.
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
void writeVTK(VTKContainer &vtk) const override
bool rightHandedness_
The right handedness of the screw, i.e. the direction of the screw-blade.
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
void write(std::ostream &os) const override
Writes an IntersectionOfWalls to an output stream, for example a restart file.
std::string getName() const final
Returns the name of the object.