MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AngledPeriodicBoundary.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 #include "AngledPeriodicBoundary.h"
26 #include "ParticleHandler.h"
27 #include "Particles/BaseParticle.h"
28 #include "Math/MatrixSymmetric.h"
30 #include "MpiContainer.h"
31 
36 {
37 #ifdef DEBUG_CONSTRUCTOR
38  std::cout << "AngledPeriodicBoundary::copy() const finished" << std::endl;
39 #endif
40  return new AngledPeriodicBoundary(*this);
41 #ifdef MERCURY_USE_MPI
42  MPIContainer& communicator = MPIContainer::Instance();
43  if (communicator.getNumberOfProcessors() > 1)
44  {
45  logger(WARN,"AngledPeriodicBoundaries are currently not implemented in parallel MercuryDPM");
46  }
47 #endif
48 
49 }
50 
59 void AngledPeriodicBoundary::set(Vec3D normalLeft, Vec3D normalRight, Vec3D origin)
60 {
61  origin_ = origin;
62  leftNormal_ = normalLeft / Vec3D::getLength(normalLeft);
63  rightNormal_ = normalRight / Vec3D::getLength(normalRight);
64  //http://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
68  Matrix3D vx = {0, -v.Z, v.Y,
69  v.Z, 0, -v.X,
70  -v.Y, v.X, 0};
71  Matrix3D eye = {1, 0, 0, 0, 1, 0, 0, 0, 1};
72  rotateLeft = eye - vx + vx * vx * ((1 - c) / s / s);
73  rotateRight = eye + vx + vx * vx * ((1 - c) / s / s);
74  //std::cout << "rotationMatrix" << rotationMatrix << std::endl;
75  commonAxis_ = v / s;
81  // walls; maybe this has to wait till quaternions are implemented.}
82  //angularShift = 0;
83  // std::cout << "common_axis " << common_axis
84  // << ", radialAxis_left " << radialAxis_left
85  // << ", radialAxis_right " << radialAxis_right
86  // << ", angularShift " << angularShift
87  // << std::endl;
88 }
89 
94 {
95  return distance(P.getPosition());
96 }
97 
102 {
103  Vec3D position = P - origin_;
104  Mdouble distance_left = Vec3D::dot(position, leftNormal_);
105  Mdouble distance_right = -Vec3D::dot(position, rightNormal_);
106 
107  if (distance_left < distance_right)
108  {
109  leftWall_ = true;
110  //std::cout << "left wall, " << position << ", distance " << distance_left << "<" << distance_right << std::endl;
111  return distance_left;
112  }
113  else
114  {
115  leftWall_ = false;
116  //std::cout << "right wall, " << position << ", distance " << distance_right << "<" << distance_left << std::endl;
117  return distance_right;
118  }
119 }
120 
125 {
126  Vec3D position = P->getPosition() - origin_;
127  if (leftWall_)
128  {
129  Mdouble normalDistance = Vec3D::dot(position, leftNormal_);
130  Mdouble radialDistance = Vec3D::dot(position, leftRadialAxis_);
131  P->move(normalDistance * differenceNormal_ + radialDistance * differenceRadialAxis_);
132  Mdouble normalVelocity = Vec3D::dot(P->getVelocity(), leftNormal_);
133  Mdouble radialVelocity = Vec3D::dot(P->getVelocity(), leftRadialAxis_);
134  P->accelerate(normalVelocity * differenceNormal_ + radialVelocity * differenceRadialAxis_);
135  Mdouble normalAngularDistance = Vec3D::dot(P->getOrientation().getAxis(), leftNormal_);
136  Mdouble radialAngularDistance = Vec3D::dot(P->getOrientation().getAxis(), leftRadialAxis_);
138  P->rotate(normalAngularDistance * differenceNormal_ + radialAngularDistance * differenceRadialAxis_);
139  Mdouble normalAngularVelocity = Vec3D::dot(P->getAngularVelocity(), leftNormal_);
140  Mdouble radialAngularVelocity = Vec3D::dot(P->getAngularVelocity(), leftRadialAxis_);
141  P->angularAccelerate(normalAngularVelocity * differenceNormal_ + radialAngularVelocity * differenceRadialAxis_);
142  leftWall_ = false;
143  //std::cout << "shift to right wall, " << P->getInteractions().size() << std::endl;
144  for (BaseInteraction* i : P->getInteractions())
145  {
146  if (i->getP()->getIndex() >= P->getIndex())
147  {
148  i->rotateHistory(rotateRight);
149  }
150  }
151  }
152  else
153  {
154  Mdouble normalDistance = Vec3D::dot(position, rightNormal_);
155  Mdouble radialDistance = Vec3D::dot(position, rightRadialAxis_);
156  P->move(-normalDistance * differenceNormal_ - radialDistance * differenceRadialAxis_);
157  Mdouble normalVelocity = Vec3D::dot(P->getVelocity(), rightNormal_);
158  Mdouble radialVelocity = Vec3D::dot(P->getVelocity(), rightRadialAxis_);
159  P->accelerate(-normalVelocity * differenceNormal_ - radialVelocity * differenceRadialAxis_);
160 
162 
163  Mdouble normalAngularVelocity = Vec3D::dot(P->getAngularVelocity(), rightNormal_);
164  Mdouble radialAngularVelocity = Vec3D::dot(P->getAngularVelocity(), rightRadialAxis_);
166  -normalAngularVelocity * differenceNormal_ - radialAngularVelocity * differenceRadialAxis_);
167  leftWall_ = true;
168  //std::cout << "shift to left wall, " << P->getPosition() << std::endl;
169  for (BaseInteraction* i : P->getInteractions())
170  {
171  if (i->getP()->getIndex() >= P->getIndex()) //if it is the deleted interaction
172  {
173  i->rotateHistory(rotateLeft);
174  }
175  }
176  }
177 }
178 
183 {
185  Vec3D position1 = P1 - origin_;
186  Vec3D position2 = P2 - origin_;
187  Mdouble normalDistance;
188  Mdouble radialDistance;
189  if (leftWall_)
190  {
191  normalDistance = Vec3D::dot(position1, leftNormal_);
192  radialDistance = Vec3D::dot(position1, leftRadialAxis_);
193  P1 += (normalDistance * differenceNormal_ + radialDistance * differenceRadialAxis_);
194  normalDistance = Vec3D::dot(position2, leftNormal_);
195  radialDistance = Vec3D::dot(position2, leftRadialAxis_);
196  P2 += (normalDistance * differenceNormal_ + radialDistance * differenceRadialAxis_);
197  leftWall_ = false;
198  }
199  else
200  {
201  normalDistance = Vec3D::dot(position1, rightNormal_);
202  radialDistance = Vec3D::dot(position1, rightRadialAxis_);
203  P1 += (-normalDistance * differenceNormal_ - radialDistance * differenceRadialAxis_);
204  normalDistance = Vec3D::dot(position2, rightNormal_);
205  radialDistance = Vec3D::dot(position2, rightRadialAxis_);
206  P2 += (-normalDistance * differenceNormal_ - radialDistance * differenceRadialAxis_);
207  leftWall_ = true;
208  }
209 }
210 
214 void AngledPeriodicBoundary::read(std::istream& is)
215 {
216  BaseBoundary::read(is);
217  std::string dummy;
218  is >> dummy >> leftNormal_ >> dummy >> rightNormal_ >> dummy >> origin_;
219  set(leftNormal_, rightNormal_, origin_);
220 }
221 
225 void AngledPeriodicBoundary::oldRead(std::istream& is)
226 {
227  std::string dummy;
228  is >> dummy >> leftNormal_ >> dummy >> rightNormal_ >> dummy >> origin_;
229  set(leftNormal_, rightNormal_, origin_);
230 }
231 
235 void AngledPeriodicBoundary::write(std::ostream& os) const
236 {
238  os << " normal_left " << leftNormal_
239  << " normal_right " << rightNormal_
240  << " origin " << origin_;
241 }
242 
247 {
248  return "AngledPeriodicBoundary";
249 }
250 
255 {
256  if (leftWall_)
257  return leftNormal_;
258  else
259  return rightNormal_;
260 }
261 
266 {
267  return acos(Vec3D::dot(leftNormal_, rightNormal_));
268 }
269 
274 {
276  if (distance(*p) < maxDistance)
277  {
278  //std::cout << "Copy particle " << p->getIndex() << " to new ghost particle" << std::endl;
279  //Step 1: Copy the particle to new ghost particle.
280  BaseParticle* pGhost = p->copy();
281 
282  //std::cout << "pGhostInteractions " << pGhost->getInteractions().size();
283  //Step 2: Copy the interactions to the ghost particle.
285  //std::cout << "-> " << pGhost->getInteractions().size() << std::endl;
286 
287  //Step 3: Shift the ghost to the 'reflected' location.
288  shiftPosition(pGhost);
289  //rotateHistory
290  //std::cout << "pGhostPosition " << pGhost->getPosition() << std::endl;
291 
292  // BaseParticle* F0 = P->copy();
293  // shiftPosition(F0);
294 
295  //Step 4: If Particle is double shifted, get correct original particle
296  BaseParticle* from = p;
297  while (from->getPeriodicFromParticle() != nullptr)
298  from = from->getPeriodicFromParticle();
299  pGhost->setPeriodicFromParticle(from);
300 
301  pH.addObject(pGhost);
302  }
303 }
304 
306 {
307  unsigned numberOfParticles = pH.getSize();
308  for (unsigned i = 0; i < numberOfParticles; i++)
309  {
311  }
312 }
313 
318 {
319  if (distance(*P) < 0)
320  {
321  shiftPosition(P);
322  }
323 }
324 
326 {
327  for (auto p = pH.begin(); p != pH.end(); ++p)
328  {
330  }
331 }
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:125
Vec3D origin_
common point of both walls
BaseParticle * getLargestParticle() const
Returns the pointer of the largest particle in the particle handler. When mercury is running in paral...
unsigned int getIndex() const
Returns the index of the object in the handler.
Definition: BaseObject.h:118
AngledPeriodicBoundary * copy() const final
Used to create a copy of the object NB: purely virtual function.
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:130
Mdouble X
the vector components
Definition: Vector.h:65
void addObject(BaseParticle *P) override
Adds a BaseParticle to the ParticleHandler.
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
Definition: BaseHandler.h:655
void shiftPositions(Vec3D &P1, Vec3D &P2)
only needed in StatisticsVector
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
void read(std::istream &is) override
reads wall
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:704
Vec3D leftRadialAxis_
outward unit normal vector for left wall
std::size_t getNumberOfProcessors() const
Get the total number of processors participating in this simulation.
void setPeriodicFromParticle(BaseParticle *p)
Assigns the pointer to the 'original' particle this one's a periodic copy of (used in periodic bounda...
Definition: BaseParticle.h:438
BaseParticle * getPeriodicFromParticle() const
Returns the 'original' particle this one's a periodic copy of.
Definition: BaseParticle.h:338
void oldRead(std::istream &is)
virtual const Vec3D & getAngularVelocity() const
Returns the angular velocity of this interactable.
void read(std::istream &is) override=0
Reads the object's id_ from given istream NB: purely virtual function, overriding the version of Base...
Definition: BaseBoundary.cc:61
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
virtual BaseParticle * copy() const =0
Particle copy method. It calls to copy constructor of this Particle, useful for polymorfism.
void accelerate(const Vec3D &vel)
Increases the particle's velocity_ by the given vector.
Mdouble getOpeningAngle()
angle between walls in radians
Stores information about interactions between two interactable objects; often particles but could be ...
Vec3D getAxis() const
Converts the quaternions into a normal vector by rotating the vector x=(1,0,0); see See Wiki for deta...
Definition: Quaternion.cc:501
void write(std::ostream &os) const override
outputs wall
const std::vector< T * >::const_iterator begin() const
Gets the begin of the const_iterator over all Object in this BaseHandler.
Definition: BaseHandler.h:690
std::string getName() const override
Returns the name of the object.
void createPeriodicParticle(BaseParticle *p, ParticleHandler &pH) override
void createPeriodicParticles(ParticleHandler &pH) override
void checkBoundaryAfterParticleMoved(BaseParticle *P)
void checkBoundaryAfterParticlesMove(ParticleHandler &pH) override
Virtual function that does things to particles, each time step after particles have moved...
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
Vec3D commonAxis_
The normalized cross product of the left and right normal vector. This vector points in the direction...
Mdouble distance(const BaseParticle &P)
Returns the distance of the wall to the particle and sets left_wall = true, if the left wall is the w...
void write(std::ostream &os) const override=0
Adds object's id_ to given ostream NB: purely virtual function, overriding the version of BaseObject...
Definition: BaseBoundary.cc:70
void set(Vec3D normalLeft, Vec3D normalRight, Vec3D origin)
Defines a periodic wall.
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
void copyInteractionsForPeriodicParticles(const BaseInteractable &p)
Copies interactions to this BaseInteractable whenever a periodic copy made.
Container to store all BaseParticle.
Mdouble Y
Definition: Vector.h:65
bool leftWall_
true if closest wall is the left wall
Vec3D rightRadialAxis_
outward unit normal vector for right wall
virtual void rotate(const Vec3D &angularVelocityDt)
Rotates this BaseInteractable.
Mdouble getLength() const
Calculates the length of this Vec3D: .
Definition: Vector.cc:320
const std::vector< BaseInteraction * > & getInteractions() const
Returns a list of interactions which belong to this interactable.
Vec3D rightNormal_
outward unit normal vector for right wall
void angularAccelerate(const Vec3D &angVel)
Increases the particle's angularVelocity_ by the given vector.
Implementation of a 3D matrix.
Definition: Matrix.h:37
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Definition: Vector.h:49
virtual void move(const Vec3D &move)
Moves this BaseInteractable by adding an amount to the position.
void shiftPosition(BaseParticle *P)
shifts the particle (after distance set the left_wall value)
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65