MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SineWall.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 "Walls/SineWall.h"
27 #include "InteractionHandler.h"
28 #include "Particles/BaseParticle.h"
29 
31 {
32  l_ = 1.0;
33  sw_wavn_ = 0.0;
34  sw_phshift_ = 0.0;
35  sw_amp_ = 0.0;
36  logger(DEBUG, "SineWall() constructor finished");
37 }
38 
39 SineWall::SineWall(const SineWall& other) : BaseWall(other)
40 {
41  l_ = other.l_;
42  sw_wavn_ = other.sw_wavn_;
43  sw_phshift_ = other.sw_phshift_;
44  sw_amp_ = other.sw_amp_;
45  logger(DEBUG, "SineWall copy constructor finished");
46 }
47 
48 SineWall::SineWall(Mdouble length, Mdouble sw_wavn, Mdouble sw_phshift, Mdouble sw_amp)
49 {
50  l_ = length;
51  sw_wavn_ = sw_wavn;
52  sw_phshift_ = sw_phshift;
53  sw_amp_ = sw_amp;
54  logger(DEBUG, "SineWall(params) constructor finished");
55 }
56 
58 {
59  logger(DEBUG, "SineWall destructor finished");
60 }
61 
62 void SineWall::set(Mdouble length, Mdouble sw_wavn, Mdouble sw_phshift, Mdouble sw_amp)
63 {
64  l_ = length;
65  sw_wavn_ = sw_wavn;
66  sw_phshift_ = sw_phshift;
67  sw_amp_ = sw_amp;
68 }
69 
71 {
72  return new SineWall(*this);
73 }
74 
75 bool SineWall::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
76 {
77  /* Define some shortcuts */
78  double x0 = p.getPosition().X;
79  double y0 = p.getPosition().Y;
80  double z0 = p.getPosition().Z;
81  double a = p.getRadius();
82 
83  // more shortcuts
84  double A = sw_amp_;
85  double k = sw_wavn_;
86  double ph = sw_phshift_;
87 
88  /* First, check whether the particle is definitely out of contact with the
89  * chute. This will be so if the particle's interaction radius is small and
90  * the particle's z-position is high. */
91  // TODO do this properly
92  if (z0 > A + a)
93  return false;
94 
95  /* Or, if the particle lies under the surface, then it is expelled
96  * perpendicularly. This shouldn't happen very often provided that the step
97  * size is small.*/
98  else if (z0 < A * sin(k * x0 + ph))
99  {
100  fprintf(stderr, "Particle at x0=%f, z0=%f is under SineWall\n", x0, z0);
101  // exit(-1);
102 
103  distance = A * sin(k * x0 + ph) - z0;
104  normal_return = Vec3D(0, 0, 1);
105  return true;
106  }
107 
108  /* Otherwise, find the contact point by using Newton's method to minimise
109  * (half of the squared) distance */
110  else
111  {
112  // Iterate according to Newton-Raphson
113  double q = x0;
114  double correction = 0;
115  // fprintf(stderr, "x0 = %lf, z0 = %lf, q = %lf\n", x0, z0, q);
116  do
117  {
118  double z = A * sin(k * q + ph);
119  double dzdq = A * k * cos(k * q + ph);
120  double d2zdq2 = -A * pow(k, 2) * sin(k * q + ph);
121 
122  // (Half of the square of the) distance. Not used.
123  double R = 0.5 * (pow(q - x0, 2) + pow(z - z0, 2));
124  // derivative of half the squared distance -- to be zeroed
125  double dRdq = q - x0 + (z - z0) * dzdq;
126  // second derivative of half the squared distance
127  double d2Rdq2 = 1 + (z - z0) * d2zdq2 + pow(dzdq, 2);
128  // The required correction to our current estimate. (Note the sign
129  // convention that this is to be subtracted.)
130  correction = dRdq / d2Rdq2;
131  //fprintf(stderr, "x0 = %.12lf, z0 = %.12lf, q = %.12lf, dRdq = %e, d2Rdq2 = %e, correction = %.e\n",
132  // x0, z0, q, dRdq, d2Rdq2, correction);
133  // N-R update
134  q -= correction;
135  // Perhaps this form is useful for reducing rounding errors?
136  // q = (q*d2Rdq2 - dRdq)/d2Rdq2;
137  } while (fabs(correction) > 1e-14);
138  //fprintf(stderr, "correction = %e\n", fabs(correction));
139 
140  // The standard Newton-Raphson iteration won't work, because the problem is
141  // too oscillatory (and N-R will converge on a near-root minimum, not the
142  // root). Neither does the following, a modified N-R which takes second and
143  // third derivatives into account.
144  /*
145  double q = x0;
146  double correction = 0;
147  do
148  {
149  // derivative of half the squared distance
150  double dRdq = q - x0 - A*pow(k,2)*z0*cos(k*q+ph) + 0.5*pow(A,2)*pow(k,3)*sin(2*k*q+2*ph);
151  // second derivative
152  double d2Rdq2 = 1 + pow(A,2)*pow(k,4)*cos(2*k*q+2*ph) + A*pow(k,3)*z0*sin(k*q+ph);
153  // third derivative
154  double d2Rdq2d = A*pow(k,4)*cos(k*q+ph) * (z0 - 4*A*k*sin(k*q+ph));
155 
156  double discriminant = d2Rdq2*d2Rdq2 - 2*dRdq*d2Rdq2d;
157  double newq = discriminant > 0 ?
158  q - (d2Rdq2 - sqrt(discriminant)) / d2Rdq2d
159  : q - dRdq / d2Rdq2;
160  double correction = newq/q - 1;
161  fprintf(stderr, "x0 = %.12lf, z0 = %.12lf, q = %.12lf, newq = %.12lf, correction = %e\n",
162  x0, z0, q, newq, correction);
163  q = newq;
164  } while (fabs(correction) > 1e-7);
165  fprintf(stderr, "converged, fabs(correction) = %e\n", fabs(correction));
166  */
167 
168  Mdouble distanceSquared
169  = pow(q - x0, 2) + pow(sw_amp_ * sin(sw_wavn_ * q + sw_phshift_) - z0, 2);
170 
171  if (distanceSquared >= mathsFunc::square(p.getWallInteractionRadius(this))
172  && z0 > sw_amp_ * sin(sw_wavn_ * q + sw_phshift_))
173  {
174  return false;
175  }
176  else
177  {
178  Vec3D ContactPoint;
179  distance = sqrt(distanceSquared);
180  ContactPoint.X = q;
181  ContactPoint.Y = y0;
182  ContactPoint.Z = sw_amp_ * sin(sw_wavn_ * q + sw_phshift_);
183  // normal_return = ContactPoint - p.getPosition();
184  normal_return = ContactPoint - p.getPosition();
185  normal_return /= normal_return.getLength();
186  return true;
187  }
188  }
189 }
190 
192 SineWall::getInteractionWith(BaseParticle* p, unsigned timeStamp, InteractionHandler* interactionHandler)
193 {
194  Mdouble distance;
195  Vec3D normal;
196  if (getDistanceAndNormal(*p, distance, normal))
197  {
198  BaseInteraction* c = interactionHandler->getInteraction(p, this, timeStamp);
199  c->setNormal(-normal);
200  c->setDistance(distance);
201  c->setOverlap(p->getRadius() - distance);
202  // c->setContactPoint( p->getPosition() - (p->getRadius() - 0.5*c->getOverlap()) * c->getNormal());
203  c->setContactPoint(p->getPosition() + distance * normal);
205  return c;
206  }
207  else
208  return nullptr;
209 }
210 
214 void SineWall::read(std::istream& is)
215 {
216  BaseWall::read(is);
217  std::string dummy;
218  is >> dummy >> l_
219  >> dummy >> sw_wavn_
220  >> dummy >> sw_phshift_
221  >> dummy >> sw_amp_;
222 }
223 
227 void SineWall::write(std::ostream& os) const
228 {
229  BaseWall::write(os);
230  os << " Length " << l_
231  << " sw_wavn " << sw_wavn_
232  << " sw_phshift " << sw_phshift_
233  << " sw_amp " << sw_amp_;
234 }
235 
239 std::string SineWall::getName() const
240 {
241  return "SineWall";
242 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
void setNormal(Vec3D normal)
Sets the normal vector between the two interacting objects.
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
void setOverlap(Mdouble overlap)
Set the overlap between the two interacting object.
double Mdouble
Definition: GeneralDefine.h:34
bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const override
Pure virtual function that computes the distance of a BaseParticle to this wall and returns the norma...
Definition: SineWall.cc:75
Mdouble sw_amp_
Definition: SineWall.h:77
void setContactPoint(Vec3D contactPoint)
Set the location of the contact point between the two interacting objects.
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:102
Mdouble l_
Definition: SineWall.h:76
void setDistance(Mdouble distance)
Sets the interaction distance between the two interacting objects.
Stores information about interactions between two interactable objects; often particles but could be ...
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 getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:383
~SineWall() override
Definition: SineWall.cc:57
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
void read(std::istream &is) override
Definition: SineWall.cc:214
BaseInteraction * getInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
Returns the Interaction between the BaseInteractable's P and I.
SineWall()
Default constructor, sets a chute with default parameters.
Definition: SineWall.cc:30
Mdouble sw_phshift_
Definition: SineWall.h:79
Container to store Interaction objects.
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
Basic class for walls.
Definition: BaseWall.h:47
SineWall * copy() const override
Pure virtual function that can be overwritten in inherited classes in order to copy a BaseWall...
Definition: SineWall.cc:70
BaseInteraction * getInteractionWith(BaseParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler) override
Returns the interaction between this wall and a given particle, nullptr if there is no interaction...
Definition: SineWall.cc:192
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
std::string getName() const override
Definition: SineWall.cc:239
Mdouble sw_wavn_
Definition: SineWall.h:78
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
void set(Mdouble length, Mdouble sw_wavn, Mdouble sw_phshift, Mdouble sw_amp)
Definition: SineWall.cc:62
Mdouble Z
Definition: Vector.h:65
void write(std::ostream &os) const override
Definition: SineWall.cc:227