MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CircularPeriodicBoundary.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 
27 #include "ParticleHandler.h"
28 #include "Particles/BaseParticle.h"
29 
31  : BaseBoundary()
32 {
33  innerRadius_ = 1.0;
34 #ifdef DEBUG_CONSTRUCTOR
35  std::cout << "CircularPeriodicBoundary::CircularPeriodicBoundary() finished" << std::endl;
36 #endif
37 #ifdef MERCURY_USE_MPI
38  MPIContainer& communicator = MPIContainer::Instance();
39  if (communicator.getNumberOfProcessors() > 1)
40  {
41  logger(WARN,"CircularPeriodicBoundaries are currently not implemented in parallel MercuryDPM");
42  }
43 #endif
44 }
45 
47  : BaseBoundary()
48 {
49  this->innerRadius_ = innerRadius;
50 #ifdef DEBUG_CONSTRUCTOR
51  std::cout << "CircularPeriodicBoundary::CircularPeriodicBoundary(double innerRadius) finished" << std::endl;
52 #endif
53 }
54 
56 {
57 #ifdef DEBUG_DESTRUCTOR
58  std::cerr << "CircularPeriodicBoundary::~CircularPeriodicBoundary() finished" << std::endl;
59 #endif
60 }
61 
63 {
64 #ifdef DEBUG_CONSTRUCTOR
65  std::cerr << "CircularPeriodicBoundary::copy() const finished" << std::endl;
66 #endif
67  return new CircularPeriodicBoundary(*this);
68 }
69 
71 {
72  double R = sqrt(pow(P->getPosition().X, 2) + pow(P->getPosition().Y, 2));
73  double alphaPos = atan2(P->getPosition().Y, P->getPosition().X);
74  double V = sqrt(pow(P->getVelocity().X, 2) + pow(P->getVelocity().Y, 2));
75  double alphaVel = atan2(P->getVelocity().Y, P->getVelocity().X);
76  alphaPos += angle;
77  alphaVel += angle;
78 
79  P->setPosition(Vec3D(cos(alphaPos) * R, sin(alphaPos * R), P->getPosition().Z));
80  P->setVelocity(Vec3D(cos(alphaVel) * V, sin(alphaVel * V), P->getVelocity().Z));
82 }
83 
85 {
86  double R = sqrt(pow(p->getPosition().X, 2) + pow(p->getPosition().Y, 2));
87  double alpha = atan2(p->getPosition().Y, p->getPosition().X);
88  unsigned int i = static_cast<unsigned int>(std::floor(std::log(R / innerRadius_) / std::log(2.0))) + 1;
89  double pieSize = 2.0 / pow(2.0, i) * constants::pi;
90  //std::cout<<"R="<<R<<" alpha="<<alpha<<" i="<<i<<" pieSize="<<pieSize<<std::endl;
91 
92  //Check if the particle is close to it's inner Radius or is is close to zero alpha (small y)
95  if (i > 0 && (R - maxDistance <
96  pow(2.0, i - 1) * innerRadius_ ||
97  p->getPosition().Y < maxDistance))
98  {
99  //std::cout<<"Going to shift because "<<R-P->getRadius()<<"<"<<pow(2,i-1)*innerRadius<<" or "<<P->getPosition().Y<<"<"<<P->getRadius()<<std::endl;
100  //std::cout<<*P<<" has been shifted"<<std::endl;
101 
102  BaseParticle* F0 = p->copy();
103  rotateParticle(F0, pieSize);
104 
105  //If Particle is Mdouble shifted, get correct original particle
106  BaseParticle* From = p;
107  while (From->getPeriodicFromParticle() != nullptr)
108  From = From->getPeriodicFromParticle();
109  F0->setPeriodicFromParticle(From);
110 
111  //std::cout<<*F0<<" is the result"<<std::endl;
112  pH.addObject(F0);
113  }
114  //Check here only for i>0 becuase for i=1 they both give the same particle
115  if (i > 1 && R * R * (1 - pow(cos(alpha - pieSize), 2)) < maxDistance)
116  {
117  //std::cout<<*P<<" has been shifted back"<<std::endl;
118 
119  BaseParticle* F0 = p->copy();
120  rotateParticle(F0, -pieSize);
121 
122  //If Particle is Mdouble shifted, get correct original particle
123  BaseParticle* From = p;
124  while (From->getPeriodicFromParticle() != nullptr)
125  From = From->getPeriodicFromParticle();
126  F0->setPeriodicFromParticle(From);
127 
128  //std::cout<<*F0<<" is the result"<<std::endl;
129  pH.addObject(F0);
130  }
131 }
132 
134 {
135  unsigned numberOfParticles = pH.getSize();
136  for (unsigned i = 0; i < numberOfParticles; i++)
137  {
139  }
140 }
141 
143 {
144  double R = sqrt(pow(P->getPosition().X, 2) + pow(P->getPosition().Y, 2));
145  double alpha = atan2(P->getPosition().Y, P->getPosition().X);
146  int i = static_cast<int>(std::floor(std::log(R / innerRadius_) / std::log(2.0))) + 1;
147  double pieSize = 2.0 / pow(2.0, i) * constants::pi;
148 
149  double oldR = sqrt(
150  pow(P->getPosition().X - P->getDisplacement().X, 2) + pow(P->getPosition().Y - P->getDisplacement().Y, 2));
152  int oldI = static_cast<int>(std::floor(std::log(oldR / innerRadius_) / std::log(2.0))) + 1;
153 
154  if (i > 0 && i > oldI) //Particle moves outward so it may have to be deleted
155  {
156  //std::cout<<"Particle="<<P->getIndex()<<" moving outward with alpha="<<alpha<<" and pieSize="<<pieSize<<" ";
157  if (alpha < 0 || alpha > pieSize)
158  {
159  if (alpha > 2.0 * pieSize)
160  {
161  //std::cout<<"and it is rotated into the pie"<<std::endl;
162  rotateParticle(P, -pieSize);
163  }
164  else
165  {
166  //Delete particle
167  //std::cout<<"and it should be deleted"<<std::endl;
168  pH.removeObject(P->getIndex());
169  return true;
170  }
171  }
172  //else
173  //std::cout<<"and nothing happens"<<std::endl;
174  }
175  else if (i >= 0 && i < oldI) //Particle moves inward so it has to be coppied
176  {
177  //std::cout<<"Particle="<<P->getIndex()<<" moving inward and is thus coppied with alpha="<<alpha<<" and pieSize="<<pieSize<<std::endl;
178  //std::cout<<"i="<<i<<" oldI="<<oldI<<" R="<<R<<" oldR="<<oldR<<std::endl;
179  //std::cout<<"Position="<<P->getPosition()<<" Displacement="<<P->get_Displacement()<<std::endl;
180  BaseParticle* F0 = P->copy();
181  F0->setDisplacement(Vec3D(0.0, 0.0, 0.0));
182  if (alpha < 0)
183  {
184  rotateParticle(P, pieSize);
185  rotateParticle(F0, 0.5 * pieSize);
186  }
187  else if (alpha < 0.5 * pieSize)
188  {
189  rotateParticle(F0, 0.5 * pieSize);
190  }
191  else if (alpha < pieSize)
192  {
193  rotateParticle(F0, -0.5 * pieSize);
194  }
195  else
196  {
197  rotateParticle(P, -pieSize);
198  rotateParticle(F0, -0.5 * pieSize);
199  }
200  pH.addObject(F0);
201  }
202  else if (i > 0 && alpha < 0)
203  {
204  //std::cout<<"Particle="<<P->getIndex()<<" i="<<i<<" R="<<R<<" alpha="<<alpha<<" positive rotated pieSize="<<pieSize<<std::endl;
205  rotateParticle(P, pieSize);
206  }
207  else if (i > 0 && alpha > pieSize)
208  {
209  //std::cout<<"Particle="<<P->getIndex()<<" i="<<i<<" R="<<R<<" alpha="<<alpha<<" negative rotated pieSize="<<pieSize<<std::endl;
210  rotateParticle(P, -pieSize);
211  }
212  return false;
213 }
214 
216 {
217  for (auto p = pH.begin(); p != pH.end(); ++p)
218  {
220  {
221  p--;
222  }
223  }
224 }
225 
226 void CircularPeriodicBoundary::read(std::istream& is)
227 {
228  BaseBoundary::read(is);
229  std::string dummy;
230  is >> dummy >> innerRadius_;
231 }
232 
233 void CircularPeriodicBoundary::oldRead(std::istream& is)
234 {
235  std::string dummy;
236  is >> dummy >> innerRadius_;
237 }
238 
239 void CircularPeriodicBoundary::write(std::ostream& os) const
240 {
242  os << " innerRadius " << innerRadius_;
243 }
244 
246 {
247  return "CircularPeriodicBoundary";
248 }
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:125
void read(std::istream &is) override
reads the CircularPeriodicBoundary
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
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 removeObject(unsigned int index) override
Removes a BaseParticle from the ParticleHandler.
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
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
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
CircularPeriodicBoundary()
default constructor
double Mdouble
Definition: GeneralDefine.h:34
void setDisplacement(const Vec3D &disp)
Sets the particle's displacement (= difference between current position and that of the previous time...
std::string getName() const override
Returns the name of the object.
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
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
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 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
virtual BaseParticle * copy() const =0
Particle copy method. It calls to copy constructor of this Particle, useful for polymorfism.
Mdouble log(Mdouble Power)
used to create a circular periodic boundary
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
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
void rotateParticle(BaseParticle *P, double angle)
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
~CircularPeriodicBoundary() override
destructor
void checkBoundaryAfterParticlesMove(ParticleHandler &pH) override
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
Container to store all BaseParticle.
Mdouble Y
Definition: Vector.h:65
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
void createPeriodicParticle(BaseParticle *p, ParticleHandler &pH) override
CircularPeriodicBoundary * copy() const override
Used to create a copy of the object NB: purely virtual function.
const Vec3D & getDisplacement() const
Returns the particle's displacement relative to the previous time step.
Definition: BaseParticle.h:392
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
void createPeriodicParticles(ParticleHandler &pH) override
Definition: Vector.h:49
void write(std::ostream &os) const override
outputs the CircularPeriodicBoundary
Mdouble Z
Definition: Vector.h:65
bool checkBoundaryAfterParticleMoved(BaseParticle *P, ParticleHandler &pH)