MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HorizontalScrew.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 "HorizontalScrew.h"
27 #include "InteractionHandler.h"
28 #include "Math/ExtendedMath.h"
29 #include "Particles/BaseParticle.h"
30 using mathsFunc::square;
31 
38 {
39  start_.setZero();
40  l_ = 1.0;
41  minR_ = 0.0;
42  lowerR_ = 1.0;
43  diffR_ = 0.0;
44  n_ = 1.0;
45  omega_ = 1.0;
46  offset_ = 0.0;
47  thickness_ = 0.0;
48  bladeLength_=0;
49  bladeWidth_=0;
50  bladeMounts_={};
51  logger(DEBUG, "HorizontalScrew() constructor finished.");
52 }
53 
58  : BaseWall(other)
59 {
60  start_ = other.start_;
61  l_ = other.l_;
62  minR_ = other.minR_;
63  lowerR_ = other.lowerR_;
64  diffR_ = other.diffR_;
65  n_ = other.n_;
66  omega_ = other.omega_;
67  thickness_ = other.thickness_;
68  offset_ = other.offset_;
72  logger(DEBUG, "HorizontalScrew(const HorizontalScrew&) copy constructor finished.");
73 }
74 
85 HorizontalScrew::HorizontalScrew(Vec3D start, Mdouble l, Mdouble minR, Mdouble lowerR, Mdouble diffR, Mdouble n, Mdouble omega, Mdouble thickness, const ParticleSpecies* s)
86 {
87  start_ = start;
88  l_ = l;
89  minR_ = minR;
90  lowerR_ = lowerR;
91  diffR_ = diffR;
92  n_ = n;
93  omega_ = omega;
94  thickness_ = thickness;
95  offset_ = 0.0;
96  bladeLength_=0;
97  bladeWidth_=0;
98  bladeMounts_={};
99  setSpecies(s);
100  logger(DEBUG, "HorizontalScrew(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");
101 }
102 
104 {
105  logger(DEBUG, "~HorizontalScrew() finished, destroyed the HorizontalScrew.");
106 }
107 
112 {
113  return new HorizontalScrew(*this);
114 }
115 
127 bool HorizontalScrew::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
128 {
129  Mdouble RSquared = square(p.getPosition().X - start_.X) + square(p.getPosition().Y - start_.Y);
130  Mdouble Z = p.getPosition().Z - start_.Z;
131  //first do a simple check if particle is within the cylindrical hull of the screw
132  Mdouble maxR = std::max(minR_,lowerR_+diffR_*Z/l_)+bladeWidth_;//todo: fix
133  Mdouble interactionRadius = p.getWallInteractionRadius(this) + thickness_;
134  if (RSquared > square(maxR + interactionRadius)
135  || Z > l_ + interactionRadius
136  || Z < - interactionRadius)
137  {
138  return false;
139  }
140 
141  //else:
142  Mdouble R = sqrt(RSquared);
143  Mdouble A = atan2(p.getPosition().Y - start_.Y, p.getPosition().X - start_.X);
144 
145  //after subtracting the start position and transforming the particle position from (XYZ) into (RAZ)
146  //coordinates, we compute the distance to the wall at, located at (r,a,z)=(r,2*pi*(offset+N*q+k/2),q*L), 0<q<1.
147 
148  //To find the contact point we have to minimize (with respect to r and q)
149  //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
150  //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
151  //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
152 
153  //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
154  //Differentiate with respect to r and solve for zero:
155  //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
156  //r=R*cos(A-2*pi*(offset+N*q))
157 
158  //Substitue back
159  //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
160  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
161 
162  //So we have to minimize:
163  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
164  //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])
165  //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
166  //For this we use the Euler algoritm
167 
168  Mdouble q; //Current guess
169  Mdouble dd; //Derivative at current guess
170  Mdouble ddd; //Second derivative at current guess
171  Mdouble q0 = Z / l_; //assume closest point q0 is at same z-location as the particle
172 
173  //Set initial q to the closest position on the screw with the same angle as A
174  //The initial guess will be in the minimum of the sin closest to q0
175  //Minima of the sin are at
176  //A-2*pi*(offset+N*q)=k*pi (k=integer)
177  //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
178 
179  Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
180  q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
181 
182  //this makes the trioliet screw unique: only one turn
183  if (((int)k)%2==0)
184  return false;
185 
186  //Now apply Newton's method to find the rel height q of the contact point
187  do
188  {
189  Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
190  dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (Z - q * l_);
191  ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
192  q -= dd / ddd;
193  } while (fabs(dd / ddd) > 1e-14);
194 
195  //Calculate r
196  Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
197  maxR = std::max(minR_,lowerR_+diffR_*q);//todo: fix
198  for (auto b : bladeMounts_)
199  {
200  Mdouble dq = q-b;
201  if (dq>0 && dq<bladeLength_)
202  {
203  maxR += bladeWidth_*(dq/bladeLength_);
204  }
205  }
206 
207  //Check if the location is actually on the screw:
208  //First possibility is that the radius is too large:
209  if (fabs(r) > maxR) //Left boundary of the coil
210  {
211  //either the contact point is too far from the screw ...
212  if (fabs(r) > maxR+thickness_)
213  return false;
214  //... or we have to compute the contact around the edge
215  r = mathsFunc::sign(r) * maxR;
216  unsigned int steps = 0;
217  //This case reduces to the coil problem
218  do
219  {
220  dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) - 2.0 * l_ * (Z - q * l_);
221  ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) + 2.0 * l_ * l_;
222  q -= dd / ddd;
223  steps++;
224  } while (fabs(dd / ddd) > 1e-14);
225  }
226  //Second possibility is that it occurred before the start of after the end
227  if (q < 0)
228  {
229  q = 0;
230  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
231  if (fabs(r) > maxR)
232  {
233  r = mathsFunc::sign(r) * maxR;
234  }
235  }
236  else if (q > 1)
237  {
238  q = 1;
239  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
240  if (fabs(r) > maxR)
241  {
242  r = mathsFunc::sign(r) * maxR;
243  }
244  }
245 
246  Mdouble distanceSquared = square(R*sin(A - 2 * constants::pi * (offset_ + n_ * q))) + square(Z - q * l_);
247  //If distance is too large there is no contact
248  if (distanceSquared >= square(interactionRadius))
249  {
250  return false;
251  }
252 
253  Vec3D ContactPoint;
254  distance = sqrt(distanceSquared) - thickness_;
255  ContactPoint.X = start_.X + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
256  ContactPoint.Y = start_.Y + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
257  ContactPoint.Z = start_.Z + q * l_;
258  normal_return = (ContactPoint - p.getPosition());
259  normal_return.normalise();
260  return true;
261 }
262 
264 {
265  return l_;
266 }
267 
269 {
270  return start_;
271 }
272 
277 {
278  offset_ += omega_ * dt;
279 }
280 
284 void HorizontalScrew::read(std::istream& is)
285 {
286  BaseWall::read(is);
287  std::string dummy;
288  unsigned n=0;
289  is >> dummy >> start_
290  >> dummy >> l_
291  >> dummy >> minR_
292  >> dummy >> lowerR_
293  >> dummy >> diffR_
294  >> dummy >> n_
295  >> dummy >> omega_
296  >> dummy >> thickness_
297  >> dummy >> offset_
298  >> dummy >> bladeLength_
299  >> dummy >> bladeWidth_
300  >> dummy >> n;
301  for (unsigned i=0; i<n; i++) {
302  Mdouble val;
303  is >> val;
304  bladeMounts_.push_back(val);
305  }
306 }
307 
311 void HorizontalScrew::write(std::ostream& os) const
312 {
313  BaseWall::write(os);
314  os << " start " << start_
315  << " length " << l_
316  << " minRadius " << minR_
317  << " lowerRadius " << lowerR_
318  << " diffRadius " << diffR_
319  << " revolutions " << n_
320  << " omega " << omega_
321  << " thickness " << thickness_
322  << " offset " << offset_
323  << " bladeLength " << bladeLength_
324  << " bladeWidth " << bladeWidth_
325  << " bladeMounts " << bladeMounts_.size();
326  for (auto n : bladeMounts_)
327  os << " " << n;
328 }
329 
333 std::string HorizontalScrew::getName() const
334 {
335  return "HorizontalScrew";
336 }
337 
338 void HorizontalScrew::setBlades(const Mdouble bladeWidth, const Mdouble bladeLength,const std::vector<Mdouble> bladeMounts)
339 {
340  bladeWidth_ = bladeWidth;
341  bladeLength_ = bladeLength;
342  bladeMounts_ = bladeMounts;
343 }
344 
346 {
347  unsigned nr = 10;
348  unsigned nz = 180;
349 
350  unsigned nPoints = vtk.points.size();
351  Vec3D contactPoint;
352  for (unsigned iz=0; iz<nz; iz++) {
353  double maxR = std::max(minR_,lowerR_+(double)iz/nz*diffR_);//todo: fix
354  for (auto b : bladeMounts_)
355  {
356  Mdouble dq = (double)iz/nz-b;
357  if (dq>0 && dq<bladeLength_)
358  {
359  maxR += bladeWidth_*(dq/bladeLength_);
360  }
361  }
362  for (unsigned ir=0; ir<nr; ir++) {
363  double q = (double)iz/nz;
364  double r = (double)ir/nr*maxR;
365  contactPoint.X = start_.X - r * cos(2.0 * constants::pi * (offset_ + n_ * q));
366  contactPoint.Y = start_.Y - r * sin(2.0 * constants::pi * (offset_ + n_ * q));
367  contactPoint.Z = start_.Z + q * l_;
368  vtk.points.push_back(contactPoint);
369  }
370  }
371 
372  unsigned nCells = vtk.triangleStrips.size();
373  //vtk.triangleStrips.reserve(nCells+(nz-1));
374  for (unsigned iz=0; iz<nz-1; iz++) {
375  std::vector<double> cell;
376  cell.reserve(2*nr);
377  for (unsigned ir=0; ir<nr; ir++) {
378  cell.push_back(nPoints+ir+iz*nr);
379  cell.push_back(nPoints+ir+(iz+1)*nr);
380  }
381  vtk.triangleStrips.push_back(cell);
382  }
383 }
~HorizontalScrew()
Default destructor.
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
void read(std::istream &is) override
Reads a HorizontalScrew from an input stream, for example a restart file.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
HorizontalScrew * copy() const final
Copy this screw and return a pointer to the copy.
HorizontalScrew()
Default constructor: make a screw with default parameters.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Vec3D getStart() const
Returns the starting position of the HorizontalScrew.
bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const final
Compute the distance from the HorizontalScrew for a given BaseParticle and return if there is a colli...
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
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
void writeVTK(VTKContainer &vtk) const override
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
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
void setBlades(const Mdouble bladeWidth, const Mdouble bladeLength, const std::vector< Mdouble > bladeMounts)
Mdouble getLength() const
Returns the length of the HorizontalScrew.
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
Mdouble l_
The length of the HorizontalScrew.
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Mdouble thickness_
The thickness of the HorizontalScrew.
This function defines an Archimedes' screw in the z-direction from a (constant) starting point...
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
Mdouble Y
Definition: Vector.h:65
void write(std::ostream &os) const override
Writes this HorizontalScrew to an output stream, for example a restart file.
std::vector< Vec3D > points
Definition: BaseWall.h:38
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:600
void move_time(Mdouble dt)
Rotate the HorizontalScrew for a period dt, so that the offset_ changes with omega_*dt.
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble omega_
Rotation speed in rev/s.
Mdouble n_
The number of revelations.
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
Mdouble minR_
The outer radius of the HorizontalScrew.
Mdouble Z
Definition: Vector.h:65
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171
std::string getName() const final
Returns the name of the object, here the string "HorizontalScrew".