MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Screw.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 "Screw.h"
27 #include "array"
28 #include "InteractionHandler.h"
29 #include "WallHandler.h"
30 #include "DPMBase.h"
31 #include "Math/ExtendedMath.h"
32 #include "Particles/BaseParticle.h"
33 
34 using mathsFunc::square;
35 
42 {
43  start_.setZero();
44  l_ = 1.0;
45  maxR_ = 1.0;
46  n_ = 1.0;
47  omega_ = 1.0;
48  offset_ = 0.0;
49  thickness_ = 0.0;
51  setOrientationViaNormal({0, 0, 1});//default screw is in z-direction
52  logger(DEBUG, "Screw() constructor finished.");
53 }
54 
58 Screw::Screw(const Screw& other)
59  : BaseWall(other)
60 {
61  start_ = other.start_;
62  l_ = other.l_;
63  maxR_ = other.maxR_;
64  n_ = other.n_;
65  omega_ = other.omega_;
66  thickness_ = other.thickness_;
67  offset_ = other.offset_;
68  screwType_ = other.screwType_;
69  logger(DEBUG, "Screw(const Screw&) copy constructor finished.");
70 }
71 
82 Screw::Screw(Vec3D start, Mdouble l, Mdouble r, Mdouble n, Mdouble omega, Mdouble thickness, ScrewType screwType)
83 {
84  start_ = start;
85  l_ = l;
86  maxR_ = r;
87  n_ = n;
88  omega_ = omega;
89  thickness_ = thickness;
90  offset_ = 0.0;
91  screwType_ = screwType;
92  logger(DEBUG, "Screw(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");
93 }
94 
96 {
97  logger(DEBUG, "~Screw() finished, destroyed the Screw.");
98 }
99 
104 {
105  return new Screw(*this);
106 }
107 
108 bool Screw::getDistanceAndNormal(const BaseParticle& P, Mdouble& distance, Vec3D& normal_return) const
109 {
110  //transform coordinates into position-orientation frame
111  Vec3D position = P.getPosition() - getPosition();
112  getOrientation().rotateBack(position);
115  if (getDistanceAndNormalLabCoordinates(position, P.getRadius() + s->getInteractionDistance(), distance,
116  normal_return))
117  {
118  getOrientation().rotate(normal_return);
119  return true;
120  }
121  else
122  {
123  return false;
124  }
125 }
126 
138 bool Screw::getDistanceAndNormalLabCoordinates(Vec3D position, Mdouble wallInteractionRadius, Mdouble& distance,
139  Vec3D& normal_return) const
140 {
141  // compute the square of the radial distance (in yz-plane) between particle and screw.
142  const Mdouble RSquared = square(position.Y - start_.Y) + square(position.Z - start_.Z);
143  const Mdouble X = position.X - start_.X;
144 
145  //first do a simple check if particle is within the cylindrical hull of the screw
146  if (RSquared > square(maxR_ + wallInteractionRadius + thickness_)
147  || X > l_ + wallInteractionRadius + thickness_
148  || X < -wallInteractionRadius - thickness_)
149  {
150  return false;
151  }
152 
153  //else compute radial distance, and angle in yz-plane
154  const Mdouble R = sqrt(RSquared);
155  const Mdouble A = atan2(position.Z - start_.Z, position.Y - start_.Y);
156 
157  //after subtracting the start position and transforming the particle position from (XYZ) into (XRA)
158  //coordinates, we compute the distance to the wall at, located at (r,a,z)=(q*L,r,2*pi*(offset+N*q+k/2)), 0<q<1.
159 
160  //To find the contact point we have to minimize (with respect to r and q)
161  //distance^2=(x-x0-r*cos(2*pi*(offset+N*q)))^2+(y-y0-r*sin(2*pi*(offset+N*q)))^2+(z-z0-q*L)^2
162  //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
163  //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
164 
165  //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
166  //Differentiate with respect to r and solve for zero:
167  //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
168  //r=R*cos(A-2*pi*(offset+N*q))
169 
170  //Substitue back
171  //distance^2=R^2+R^2*cos^2(A-2*pi*(offset+N*q))-2*R^2*cos^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
172  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
173 
174  //So we have to minimize:
175  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
176  //f'(q)=(-2*pi*N)*R^2*sin(2*A-4*pi*(offset+N*q)) + 2*L*(Z-q*L) (D[Sin[x]^2,x]=Sin[2x])
177  //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
178  //For this we use the Euler algoritm
179 
180  Mdouble q; //Current guess
181  Mdouble dd; //Derivative at current guess
182  Mdouble ddd; //Second derivative at current guess
183 
184  //Set initial q to the closest position on the screw with the same angle as A
185  //The initial guess will be in the minimum of the sin closest to q0
186  //Minima of the sin are at
187  //A-2*pi*(offset+N*q)=k*pi (k=integer)
188  //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
189  {
190  const Mdouble q0 = X / l_; //assume closest point q0 is at same z-location as the particle
191  const Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
192  q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
193  }
194 
195  //Now apply Newton's method
196  unsigned count = 0;
197  const unsigned maxCount = 30;
198  do
199  {
200  Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
201  dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (X - q * l_);
202  ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
203  q -= dd / ddd;
204  if(++count>maxCount) {
205  logger(WARN,"Screw % contact detection did not converge (err=%); continuing with approximate contact point",getId(),fabs(dd/ddd));
206  break;
207  }
208  } while (fabs(dd / ddd) > 5e-14);
209 
210  //Calculate r
211  Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
212 
213  //If the particle touches the single screw center
214  if (screwType_==ScrewType::singleHelix && r>0) {
215  distance = R-thickness_;
216  if (distance>wallInteractionRadius) {
217  return false;
218  } else {
219  normal_return = Vec3D(0,-position.Y/R,-position.Z/R);
220  return true;
221  }
222  }
223 
224  //Check if the location is actually on the screw:
225  //First posibility is that the radius is too large:
226  Mdouble distanceSquared = 0;
227  if (fabs(r) > maxR_) //Left boundary of the coil
228  {
229  r = mathsFunc::sign(r) * maxR_;
230  distanceSquared = mathsFunc::square(R-maxR_);
231  count = 0;
232  //This case reduces to the coil problem
233  do
234  {
235  dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) -
236  2.0 * l_ * (X - q * l_);
237  ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) +
238  2.0 * l_ * l_;
239  q -= dd / ddd;
240  if(++count>maxCount) {
241  logger(WARN,"Screw % contact detection at left boundary did not converge (err=%); continuing with approximate contact point",getId(),fabs(dd/ddd));
242  break;
243  }
244  } while (fabs(dd / ddd) > 1e-13);
245  }
246  //Second possibility is that it occured before the start of after the end
247  if (q < 0)
248  {
249  q = 0;
250  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
251  if (fabs(r) > maxR_)
252  {
253  r = mathsFunc::sign(r) * maxR_;
254  }
255  }
256  else if (q > 1)
257  {
258  q = 1;
259  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
260  if (fabs(r) > maxR_)
261  {
262  r = mathsFunc::sign(r) * maxR_;
263  }
264  }
265 
266  //compute squared distance between particle and contact point
267  distanceSquared += square(X - q * l_) + square(R*sin(A - 2 * constants::pi * (offset_ + n_ * q)));
268 
269  //If distance is too large there is no contact
270  if (distanceSquared >= square(wallInteractionRadius + thickness_))
271  {
272  return false;
273  }
274 
275  //Else
276  Vec3D ContactPoint;
277  distance = sqrt(distanceSquared) - thickness_;
278  ContactPoint.Y = start_.Y + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
279  ContactPoint.Z = start_.Z + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
280  ContactPoint.X = start_.X + q * l_;
281  normal_return = (ContactPoint - position);
282  normal_return.normalise();
283  return true;
284 }
285 
290 {
291  //offset_ += omega_ * dt;
292 }
293 
298 void Screw::rotate(const Vec3D& angularVelocityDt)
299 {
300  BaseInteractable::rotate(angularVelocityDt);
302 }
303 
304 
308 void Screw::read(std::istream& is)
309 {
310  BaseWall::read(is);
311  std::string dummy;
312  unsigned screwType;
313  is >> dummy >> start_
314  >> dummy >> l_
315  >> dummy >> maxR_
316  >> dummy >> n_
317  >> dummy >> omega_
318  >> dummy >> thickness_
319  >> dummy >> offset_
320  >> dummy >> screwType;
321  screwType_ = static_cast<ScrewType>(screwType);
322 }
323 
330 void Screw::oldRead(std::istream& is)
331 {
332  std::string dummy;
333  is >> dummy >> start_
334  >> dummy >> l_
335  >> dummy >> maxR_
336  >> dummy >> n_
337  >> dummy >> omega_
338  >> dummy >> offset_;
339 }
340 
344 void Screw::write(std::ostream& os) const
345 {
346  BaseWall::write(os);
347  os << " start " << start_
348  << " Length " << l_
349  << " Radius " << maxR_
350  << " Revolutions " << n_
351  << " Omega " << omega_
352  << " Thickness " << thickness_
353  << " Offset " << offset_
354  << " ScrewType " << static_cast<unsigned>(screwType_);
355 }
356 
360 std::string Screw::getName() const
361 {
362  return "Screw";
363 }
364 
366 {
367  //number of points in radial direction (for one side of the screw)
368  unsigned nr = 5;
369  //number of points in axial direction
370  unsigned nz = static_cast<unsigned int>(99 * fabs(n_));
371 
372  unsigned long nPoints = vtk.points.size();
373  Vec3D contactPoint;
374  // either one or two helices
375  for (Mdouble offset = offset_; offset<=offset_+(screwType_==ScrewType::doubleHelix?0.5:0); offset+=0.5)
376  {
377  for (unsigned iz = 0; iz < nz; iz++)
378  {
379  for (unsigned ir = 0; ir < nr; ir++)
380  {
381  double q = (double) iz / nz;
382  double r = (double) ir / (nr - 1) * maxR_;
383  contactPoint.Y = start_.Y - r * cos(2.0 * constants::pi * (offset + n_ * q));
384  contactPoint.Z = start_.Z - r * sin(2.0 * constants::pi * (offset + n_ * q));
385  contactPoint.X = start_.X + q * l_ - thickness_;
386  getOrientation().rotate(contactPoint);
387  contactPoint += getPosition();
388  vtk.points.push_back(contactPoint);
389  }
390  for (unsigned ir = 0; ir < nr; ir++)
391  {
392  double q = (double) iz / nz;
393  double r = (double) (nr - 1 - ir) / (nr - 1) * maxR_;
394  contactPoint.Y = start_.Y - r * cos(2.0 * constants::pi * (offset + n_ * q));
395  contactPoint.Z = start_.Z - r * sin(2.0 * constants::pi * (offset + n_ * q));
396  contactPoint.X = start_.X + q * l_ + thickness_;
397  getOrientation().rotate(contactPoint);
398  contactPoint += getPosition();
399  vtk.points.push_back(contactPoint);
400  }
401  }
402  }
403 
404  unsigned long nCells = vtk.triangleStrips.size();
405  //vtk.triangleStrips.reserve(nCells + (nz - 1)*(screwType_==ScrewType::doubleHelix?2:1));
406  for (unsigned iz = 0; iz < (screwType_==ScrewType::doubleHelix?2:1)*nz-1; iz++)
407  {
408  //skip step that would connect the two screw parts
409  if (iz==nz-1) ++iz;
410  std::vector<double> cell;
411  cell.reserve(2 * nr);
412  for (unsigned ir = 0; ir < 2*nr; ir++)
413  {
414  cell.push_back(nPoints + ir + 2*iz * nr);
415  cell.push_back(nPoints + ir + 2*(iz + 1) * nr);
416  }
417  vtk.triangleStrips.push_back(cell);
418  }
419 }
420 
421 void Screw::writeVTK(std::string filename) const
422 {
423  VTKContainer vtk;
424  writeVTK(vtk);
425 
426  std::stringstream file;
427  file << "# vtk DataFile Version 2.0\n"
428  << getName() << "\n"
429  "ASCII\n"
430  "DATASET UNSTRUCTURED_GRID\n"
431  "POINTS " << vtk.points.size() << " double\n";
432  for (const auto& vertex : vtk.points)
433  file << vertex << '\n';
434  file << "\nCELLS " << vtk.triangleStrips.size() << ' ' << 4 * vtk.triangleStrips.size() << "\n";
435  for (const auto& face : vtk.triangleStrips)
436  file << "3 " << face[0] << ' ' << face[1] << ' ' << face[2] << '\n';
437  file << "\nCELL_TYPES " << vtk.triangleStrips.size() << "\n";
438  for (const auto& face : vtk.triangleStrips)
439  file << "5\n";
440  helpers::writeToFile(filename, file.str());
441 }
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: Helpers.cc:446
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
BaseSpecies is the class from which all other species are derived.
Definition: BaseSpecies.h:49
void read(std::istream &is) override
Reads a Screw from an input stream, for example a restart file.
Definition: Screw.cc:308
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
This function defines an Archimedes' screw in the z-direction from a (constant) starting point...
Definition: Screw.h:43
void writeVTK(VTKContainer &vtk) const override
Definition: Screw.cc:365
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
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
bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const final
Compute the distance from the Screw for a given BaseParticle and return if there is a collision...
Definition: Screw.cc:108
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:102
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
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
void move_time(Mdouble dt)
Rotate the Screw for a period dt, so that the offset_ changes with omega_*dt.
Definition: Screw.cc:289
Mdouble getInteractionDistance() const
returns the largest separation distance at which adhesive short-range forces can occur.
Definition: BaseSpecies.h:146
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
std::string getName() const final
Returns the name of the object, here the string "Screw".
Definition: Screw.cc:360
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
const Mdouble pi
Definition: ExtendedMath.h:45
~Screw() override
Default destructor.
Definition: Screw.cc:95
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Screw * copy() const final
Copy this screw and return a pointer to the copy.
Definition: Screw.cc:103
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
bool getDistanceAndNormalLabCoordinates(Vec3D position, Mdouble wallInteractionRadius, Mdouble &distance, Vec3D &normal_return) const
Definition: Screw.cc:138
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:80
void write(std::ostream &os) const override
Writes this Screw to an output stream, for example a restart file.
Definition: Screw.cc:344
Mdouble Y
Definition: Vector.h:65
std::vector< Vec3D > points
Definition: BaseWall.h:38
void rotate(const Vec3D &angularVelocityDt) override
Definition: Screw.cc:298
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:600
virtual void rotate(const Vec3D &angularVelocityDt)
Rotates this BaseInteractable.
Mdouble n_
The number of revelations.
Definition: Screw.h:128
ScrewType
used to define if the screw consists of a single or double helical thread
Definition: Screw.h:35
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
Screw()
Default constructor: make a screw with default parameters.
Definition: Screw.cc:41
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
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
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1211
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
void oldRead(std::istream &is)
Reads a Screw in the old style from an input stream, for example a restart file old style...
Definition: Screw.cc:330