MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Coil.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 "Coil.h"
27 #include "InteractionHandler.h"
28 #include "Particles/BaseParticle.h"
29 
36 {
37  start_.setZero();
38  l_ = 1.0;
39  r_ = 1.0;
40  n_ = 1.0;
41  omega_ = 1.0;
42  offset_ = 0.0;
43  thickness_ = 0.0;
44  logger(DEBUG, "Coil() constructor finished");
45 }
46 
50 Coil::Coil(const Coil& other) : BaseWall(other)
51 {
52  start_ = other.start_;
53  l_ = other.l_;
54  r_ = other.r_;
55  n_ = other.n_;
56  omega_ = other.omega_;
57  offset_ = other.offset_;
58  thickness_ = other.thickness_;
59  logger(DEBUG, "Coil(const Coil&) copy constructor finished");
60 }
61 
72 Coil::Coil(Vec3D start, Mdouble l, Mdouble r, Mdouble n, Mdouble omega, Mdouble thickness)
73 {
74  start_ = start;
75  l_ = l;
76  r_ = r;
77  n_ = n;
78  omega_ = omega;
79  thickness_ = thickness;
80  offset_ = 0.0;
81  logger(DEBUG, "Coil(params) constructor with parameters finished.");
82 }
83 
85 {
86  logger(DEBUG, "~Coil() destructor finished.");
87 }
88 
101 void
102 Coil::set(Vec3D start, Mdouble length, Mdouble radius, Mdouble numberOfRevelations, Mdouble omega, Mdouble thickness)
103 {
104  start_ = start;
105  l_ = length;
106  r_ = radius;
107  n_ = numberOfRevelations;
108  omega_ = omega;
109  thickness_ = thickness;
110  offset_ = 0.0;
111 }
112 
116 Coil* Coil::copy() const
117 {
118  return new Coil(*this);
119 }
120 
132 bool Coil::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
133 {
134  const Mdouble interactionRadius = p.getWallInteractionRadius(this) + thickness_;
135  Mdouble Rsqr = pow(p.getPosition().X - start_.X, 2) + pow(p.getPosition().Y - start_.Y, 2);
136  if (Rsqr > pow(r_ + interactionRadius, 2) ||
137  Rsqr < pow(r_ - interactionRadius, 2) ||
138  p.getPosition().Z > l_ + start_.Z + interactionRadius ||
139  p.getPosition().Z < start_.Z - interactionRadius)
140  {
141  //std::cout<<"Particle is out of first bound checking"<<std::endl;
142  //std::cout<<"Rule 1: "<<pow(r-P.getRadius()-thickness_,2)<<"<"<<Rsqr<<"<"<<pow(r+P.getRadius()+thickness_,2)<<std::endl;
143  //std::cout<<"Rule 2: "<<Start.Z-P.getRadius()-thickness_<<"<"<<P.getPosition().Z<<"<"<<L+Start.Z+P.getRadius()+thickness_<<std::endl;
144  return false;
145  }
146  Mdouble R = sqrt(Rsqr);
147  Mdouble alpha = atan2(p.getPosition().Y - start_.Y, p.getPosition().X - start_.X);
148  Mdouble dz = p.getPosition().Z - start_.Z;
149 
150  //To find the contact point we have to minimize (with respect to q)
151  //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
152  //Using polar coordinated (i.e. x-x0=R*cos(alpha), y-y0=R*sin(alpha) and dz=z-z0)
153  //Distance^2=R^2+r^2-2*R*r*cos(alpha-2*Pi*(offset+N*q))+(dz-q*L)^2
154 
155  //So we have to minimize:
156  //Distance^2=R^2+r^2-2*R*r*cos(alpha-2*Pi*(offset+N*q))+(dz-q*L)^2
157  //For this we use the Euler algoritm
158 
159  Mdouble q; //Current guess
160  Mdouble dd; //Derivative at current guess
161  Mdouble ddd; //Second derivative at current guess
162  Mdouble q0 = dz / l_; //Minimum of the parabolic part
163 
164  //The initial guess will be in the maximum of the cos closest to the minimum of the parabolic part
165  //Minima of the cos are at
166  //alpha-2*Pi*(offset+N*q)=2*k*Pi (k=integer)
167  //q=alpha/(2*Pi*N)-k/N-offset/N (k=integer)
168 
169  Mdouble k = round(alpha / 2.0 / constants::pi - (offset_ + n_ * q0));
170  q = alpha / (2 * constants::pi * n_) - k / n_ - offset_ / n_;
171 
172  //Now apply Newton's method
173  do
174  {
175  dd = -4.0 * R * r_ * constants::pi * n_ * sin(alpha - 2.0 * constants::pi * (n_ * q + offset_)) -
176  2.0 * l_ * (dz - q * l_);
177  ddd = 8.0 * R * r_ * constants::sqr_pi * n_ * n_ * cos(alpha - 2.0 * constants::pi * (n_ * q + offset_)) +
178  2.0 * l_ * l_;
179  q -= dd / ddd;
180  } while (fabs(dd / ddd) > 1e-14);
181 
182  //Check if the location is actually on the coil, otherwise a point collision with the end of the coil calculated
183  if (q < 0) //Left boundary of the coil
184  {
185  q = 0;
186  }
187  else if (q > 1) //right boundary of the coil
188  {
189  q = 1;
190  }
191 
192  Mdouble distanceSquared =
193  R * R + r_ * r_ - 2 * R * r_ * cos(alpha - 2 * constants::pi * (offset_ + n_ * q)) + pow(dz - q * l_, 2);
194  //If distance is too large there is no contact
195  if (distanceSquared >= mathsFunc::square(p.getWallInteractionRadius(this) + thickness_))
196  {
197  //std::cout<<"Particle is out of second bound checking, distance^2="<<distance<<" max="<<(P.getRadius()+thickness_)*(P.getRadius()+thickness_)<<std::endl;
198  return false;
199  }
200 
201  Vec3D ContactPoint;
202  distance = sqrt(distanceSquared) - thickness_;
203  ContactPoint.X = start_.X + r_ * cos(2.0 * constants::pi * (offset_ + n_ * q));
204  ContactPoint.Y = start_.Y + r_ * sin(2.0 * constants::pi * (offset_ + n_ * q));
205  ContactPoint.Z = start_.Z + q * l_;
206  normal_return = ContactPoint - p.getPosition();
207  normal_return /= normal_return.getLength();
208  return true;
209 }
210 
215 {
216  offset_ += omega_ * dt;
217 }
218 
222 void Coil::read(std::istream& is)
223 {
224  BaseWall::read(is);
225  std::string dummy;
226  is >> dummy >> start_
227  >> dummy >> l_
228  >> dummy >> r_
229  >> dummy >> n_
230  >> dummy >> omega_
231  >> dummy >> thickness_
232  >> dummy >> offset_;
233 }
234 
238 void Coil::oldRead(std::istream& is)
239 {
240  std::string dummy;
241  is >> dummy >> start_
242  >> dummy >> l_
243  >> dummy >> r_ >> dummy >> n_
244  >> dummy >> omega_
245  >> dummy >> offset_;
246 }
247 
251 void Coil::write(std::ostream& os) const
252 {
253  BaseWall::write(os);
254  os << " Start " << start_
255  << " Length " << l_
256  << " Radius " << r_
257  << " Revolutions " << n_
258  << " Omega " << omega_
259  << " Thickness " << thickness_
260  << " Offset " << offset_;
261 }
262 
266 std::string Coil::getName() const
267 {
268  return "Coil";
269 }
~Coil() override
Default destructor.
Definition: Coil.cc:84
Vec3D start_
The centre of the lower end of the Coil.
Definition: Coil.h:111
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
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 set(Vec3D Start, Mdouble length, Mdouble radius, Mdouble numberOfRevelations, Mdouble omega, Mdouble thickness)
Set all parameters of this Coil.
Definition: Coil.cc:102
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble offset_
Definition: Coil.h:136
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
void read(std::istream &is) override
Reads a Coil from an input stream, for example a restart file.
Definition: Coil.cc:222
const Mdouble pi
Definition: ExtendedMath.h:45
void write(std::ostream &os) const override
Writes a Coil to an output stream, for example a restart file.
Definition: Coil.cc:251
bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const override
Compute the distance from the Coil for a given BaseParticle and return if there is a collision...
Definition: Coil.cc:132
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
Mdouble l_
The length of the Coil.
Definition: Coil.h:116
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:600
Mdouble n_
Definition: Coil.h:126
void move_time(Mdouble dt)
Rotate the Coil for a period dt, so that the offset_ changes with omega_*dt.
Definition: Coil.cc:214
Mdouble r_
Definition: Coil.h:121
Coil * copy() const override
Copy this Coil and return a pointer to the copy, useful for polymorphism.
Definition: Coil.cc:116
Coil()
Default constructor, sets a coil with default parameters.
Definition: Coil.cc:35
Mdouble thickness_
Definition: Coil.h:141
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
MERCURY_DEPRECATED void oldRead(std::istream &is)
Reads an old-style Coil from an input stream, for example an old restart file.
Definition: Coil.cc:238
Mdouble Z
Definition: Vector.h:65
This class defines a coil in the z-direction from a (constant) starting point, a (constant) length L...
Definition: Coil.h:40
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
Mdouble omega_
Definition: Coil.h:131
std::string getName() const override
Returns the name of the object, in this case the string "Coil".
Definition: Coil.cc:266