inflowFromPeriodic.h
Go to the documentation of this file.
1 //Copyright (c) 2013-2023, 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 "scr/Chute.h"
27 
28 class inflowFromPeriodic : public Chute
29 {
30 public:
32  //load particles and chute setup into periodic chute
33  setName("../../H20A22L0.5M0.5");
34  load_restart_data();
35  setRestarted(false);
36 
37  // adds new species
39 
40  // creates bottom outside periodic chute of species 1
41  set_Nmax(get_N()*12.);
42  int N=get_N();
43  for (int j=0; j<5; j++) {
44  for (int i=0; i<N; i++) {
45  if (Particles[i].is_fixed()) {
46  Particles.push_back(Particles[i]);
47  Particles.back().indSpecies=1;
48  Particles.back().Position.X+=20*j;
49  }
50  }
51  }
52 
53  Chute SlowBottom;
54  SlowBottom.setName("../../H20A22L0.75M0.5");
55  SlowBottom.load_restart_data();
56  for (int j=5; j<10; j++) {
57  for (int i=0; i<SlowBottom.get_N(); i++) {
58  if (SlowBottom.getObjects()[i].is_fixed()) {
59  Particles.push_back(SlowBottom.getObjects()[i]);
60  Particles.back().indSpecies=1;
61  Particles.back().Position.X+=20*j;
62  }
63  }
64  }
65  setXMax(10*getXMax());
66 
67  Walls.resize(0);
68  }
69 
70  inflowFromPeriodic(string restart_file) {
71  //load particles and chute setup into periodic chute
72  setName(restart_file.c_str());
73  load_restart_data();
74  setName("restarted");
75  set_HGRID_num_buckets_to_power();
76  }
77 
78  //Do not add or remove particles
79  void actionsBeforeTimeStep() override {cleanChute();};
80 
81  void cleanChute() {
82  //clean outflow every 100 time steps
83  static int count = 0, maxcount = 100;
84  if (count>maxcount)
85  {
86  count = 0;
87  // delete all outflowing particles
88  for (unsigned int i=0;i<Particles.size();)
89  {
90  if (Particles[i].Position.Z<getZMin()-10)//||Particles[i].Position.Z+Particles[i].Radius<zmin)
91  {
92  #ifdef DEBUG_OUTPUT_FULL
93  cout << "erased:" << Particles[i] << endl;
94  #endif
95  removeParticle(i);
96  }
97  else i++;
98  }
99  } else count++;
100  }
101 
102  //Do not add bottom
103  void setupInitialConditions() override {}
104 
106  {
107  Particles[i].Velocity += Particles[i].Force*Particles[i].get_invmass()*0.5*getTimeStep();
108  Particles[i].Position += Particles[i].Velocity*getTimeStep();
109  if (speciesHandler.getObject(0)->mu)
110  {
111  Particles[i].AngularVelocity += Particles[i].Torque*Particles[i].get_invinertia()*0.5*getTimeStep();
112  Particles[i].Angle += Particles[i].AngularVelocity*getTimeStep();
113  }
114 
115  // This shifts particles that moved through periodic walls
116  if (Particles[i].indSpecies==0 && WallsPeriodic[0].distance(Particles[i])<0) {
117  if (Particles[i].Position.X>18 && Particles[i].Position.X<22) {
118  if (get_Nmax()<=get_N()) set_Nmax(get_Nmax()+10000);
119  Particles.push_back(Particles[i]);
120  Particles.back().indSpecies=1;
121  }
122  WallsPeriodic[0].shift_position(Particles[i].Position);
123  }
124  if (WallsPeriodic[1].distance(Particles[i])<0)
125  WallsPeriodic[1].shift_position(Particles[i].Position);
126  }
127 
128  //~ void computeExternalForces(int CI) {
129  //~ MD::computeExternalForces(CI);
130  //limit force
131  //~ if (Particles[CI].indSpecies==1 && Particles[CI].Position.X<40) {
132  //~ Particles[CI].Force *= max(0.,(Particles[CI].Position.X-22.)/18.);
133  //~ if (std::max(0.,(Particles[CI].Position.X-22.)/18.))
134  //~ cout << max(0.,(Particles[CI].Position.X-22.)/18.) << endl;
135  //~ }
136  //~ }
137 
138  void printTime() {
139  int counter=0;
140  for (int i=0; i<get_N(); i++)
141  if (Particles[i].indSpecies==0) counter++;
142 
143  cout << "t=" << setprecision(3) << left << setw(6) << getTime()
144  << ", Nmax=" << setprecision(3) << left << setw(6) << get_Nmax()
145  << ", counter=" << setprecision(3) << left << setw(6) << counter << endl;
146  }
147 
148  void add_forces(int PI, int PJreal) {
149  // Add and subtract this force to the force acting on particles i, j
150  if (Particles[PI ].indSpecies!=Particles[PJreal].indSpecies) {
151  if (Particles[PI ].indSpecies==0) {
152  if (Particles[PI ].Position.X>18) {
153  Particles[PJreal].Force-=force;
154  }
155  } else {
156  if (Particles[PJreal].Position.X>18) {
157  Particles[PI ].Force+=force;
158  }
159  }
160  } else {
161  Particles[PI ].Force+=force;
162  Particles[PJreal].Force-=force;
163  }
164  }
165 
166  int Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)
167  {
168  int C=0; //Number of particles created
169  if (nWallPeriodic==0 && Particles[i].indSpecies==0 && WallsPeriodic[0].distance(Particles[i])<Particles[i].Radius+get_radius_of_largest_particle())
170  {
171  CParticle F0=Particles[i];
172  WallsPeriodic[0].shift_position(F0.Position);
173 
174  //If Particle is double shifted, get correct original particle
175  int From=i;
176  while(Particles[From].get_periodicFromParticle()!=-1)
177  From=Particles[From].get_periodicFromParticle();
178  F0.set_periodicFromParticle(From);
179 
180  Particles.push_back(F0);
181  HGridInsertParticleToHgrid(Particles[get_N()-1]);
182  C++;
183 
184  //Check for double shifted particles
185  C+=Check_and_Duplicate_Periodic_Particle(get_N()-1, 0+1);
186  }
187  for (int k=max(1,nWallPeriodic); k<get_NWallPeriodic(); k++) { //Loop over all still posible walls
189  if (WallsPeriodic[k].distance(Particles[i])<Particles[i].Radius+get_radius_of_largest_particle())
190  {
191  CParticle F0=Particles[i];
192  WallsPeriodic[k].shift_position(F0.Position);
193 
194  //If Particle is double shifted, get correct original particle
195  int From=i;
196  while(Particles[From].get_periodicFromParticle()!=-1)
197  From=Particles[From].get_periodicFromParticle();
198  F0.set_periodicFromParticle(From);
199 
200  Particles.push_back(F0);
201  HGridInsertParticleToHgrid(Particles[get_N()-1]);
202  C++;
203 
204  //Check for double shifted particles
205  C+=Check_and_Duplicate_Periodic_Particle(get_N()-1, k+1);
206  }
207  }
208  return(C);
209  }
210 
212  {
213 
214  stringstream file_name;
215  ofstream script_file;
216  file_name << problem_name.str() <<".disp";
217  script_file.open((file_name.str()).c_str());
218 
220  script_file << "#!/bin/bash" << endl;
221  script_file << "x=$(echo $0 | cut -c2-)" << endl;
222  script_file << "file=$PWD$x" << endl;
223  script_file << "dirname=`dirname \"$file\"`" << endl;
224  script_file << "cd $dirname" << endl;
225 
226  double scale;
227  int format;
228 
229  int width=1700-140, height=1000-140;
230  int ratio=getSystemDimensions()<3?(getXMax()-getXMin())/(getYMax()-getYMin()):(getXMax()-getXMin())/(getZMax()-getZMin());
231  if (ratio>width/height) height=width/ratio;
232  else width=height*ratio;
233 
234  if (getSystemDimensions()<3)
235  { // dim = 1 or 2
236  format = 8;
237  if (getXBallsScale()<0)
238  {
239  scale = 1.0 / max( getYMax()-getYMin(), getXMax()-getXMin() );
240  }
241  else
242  {
243  scale=getXBallsScale();
244  }
245  }
246  else
247  { //dim==3
248  format = 14;
249  if (getXBallsScale()<0)
250  {
251  scale = width/480/(getXMax()-getXMin());
252  }
253 
254  else
255  {
256  scale=getXBallsScale();
257  }
258 
259  }
260 
261  script_file << "../xballs -format " << format
262  << " -f " << data_filename.str() << ((get_options_data()==2)?"'.0000":"")
263  << " -s " << scale
264  << " -w " << width+140
265  << " -h " << height+140
266  << " -cmode " << getXBallsColourMode()
267  << " -cmax -scala 4 -sort "
269  << " $*";
270  if (getXBallsVectorScale()>-1)
271  {
272  script_file << " -vscale " << getXBallsVectorScale();
273  }
274  script_file.close();
275 
276  //This line changes teh file permision and give the owener (i.e. you) read, write and excute permission to the file.
277  chmod((file_name.str().c_str()),S_IRWXU);
278 
279  }
280 
281 
282  //~ void outputXBallsData()
283  //~ {
284  //~ format=14;
285  //~ // This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
286  //~ int counter=0;
287  //~ for (unsigned int i = 0;i<Particles.size();i++) if (Particles[i].indSpecies==0) counter++;
288  //~ if (format!=14) // dim = 1 or 2
289  //~ data_file << Particles.size()-counter << " " <<t << " "
290  //~ << xmin << " " << ymin << " "
291  //~ << xmax << " " << ymax << " " << endl;
292  //~ else //dim==3
293  //~ data_file << Particles.size()-counter << " " <<t << " "
294  //~ << xmin << " " << ymin << " " << zmin << " "
295  //~ << xmax << " " << ymax << " " << zmax << " " << endl;
296  //~ // This outputs the particle data
297  //~ for (unsigned int i = 0;i<Particles.size();i++) {
298  //~ if (Particles[i].indSpecies==1) outputXBallsDataParticle(i);
299  //~ }
300  //~ }
301 
302 };
std::enable_if<!std::is_pointer< U >::value, U * >::type copyAndAddObject(const U &object)
Creates a copy of a Object and adds it to the BaseHandler.
Definition: BaseHandler.h:379
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
Creates chutes with different bottoms. Inherits from Mercury3D (-> MercuryBase -> DPMBase).
Definition: Chute.h:65
int getXBallsColourMode() const
Get the xballs colour mode (CMode).
Definition: DPMBase.cc:1310
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin.
Definition: DPMBase.h:619
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax.
Definition: DPMBase.h:626
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1427
double getXBallsVectorScale() const
Returns the scale of vectors used in xballs.
Definition: DPMBase.cc:1330
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin.
Definition: DPMBase.h:632
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:422
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1250
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:808
std::string getXBallsAdditionalArguments() const
Returns the additional arguments for xballs.
Definition: DPMBase.cc:1355
void setRestarted(bool newRestartedFlag)
Allows to set the flag stating if the simulation is to be restarted or not.
Definition: DPMBase.cc:1501
void setXMax(Mdouble newXMax)
Sets the value of XMax, the upper bound of the problem domain in the x-direction.
Definition: DPMBase.cc:1165
double getXBallsScale() const
Returns the scale of the view in xballs.
Definition: DPMBase.cc:1372
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax.
Definition: DPMBase.h:638
unsigned int getSystemDimensions() const
Returns the system dimensionality.
Definition: DPMBase.cc:1430
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax.
Definition: DPMBase.h:650
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin.
Definition: DPMBase.h:644
Definition: inflowFromPeriodic.h:29
int Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)
Definition: inflowFromPeriodic.h:166
void integrateBeforeForceComputation(int i)
Definition: inflowFromPeriodic.h:105
void add_forces(int PI, int PJreal)
Definition: inflowFromPeriodic.h:148
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: inflowFromPeriodic.h:103
void actionsBeforeTimeStep() override
A virtual function which allows to define operations to be executed before the new time step.
Definition: inflowFromPeriodic.h:79
inflowFromPeriodic()
Definition: inflowFromPeriodic.h:31
void writeXBallsScript()
Definition: inflowFromPeriodic.h:211
void cleanChute()
Definition: inflowFromPeriodic.h:81
inflowFromPeriodic(string restart_file)
Definition: inflowFromPeriodic.h:70
void printTime()
Definition: inflowFromPeriodic.h:138
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51