MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ParticleAnalytics2Mercury.cpp
Go to the documentation of this file.
1 //Copyright (c) 2015, 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 <iostream>
27 #include <fstream>
28 #include <Logger.h>
31 #include <Walls/InfiniteWall.h>
32 #include "Mercury3D.h"
33 using constants::pi;
34 
40 {
41 public:
42 
46  explicit FileReader(std::string name) {
47  pFile_.open(name+".p3p");
48  if (pFile_) {
50  wFile_.open(name+".p3w");
51  cFile_.open(name+".p3c");
52  if (!wFile_) logger(ERROR, "Could not open %", name + ".p3w");
53  if (!cFile_) logger(ERROR, "Could not open %", name + ".p3c");
54  } else {
56  pFile_.open(name+".p4p");
57  wFile_.open(name+".p4w");
58  cFile_.open(name+".p4c");
59  if (!pFile_) logger(ERROR, "Could not open % or %", name + ".p3p",name + ".p4p");
60  if (!wFile_) logger(ERROR, "Could not open %", name + ".p4w");
61  if (!cFile_) logger(ERROR, "Could not open %", name + ".p4c");
62  }
63  logger(INFO,"Input files opened");
64 
65  dpm.setName(name);
66  }
67 
68  bool read() {
72 
73  std::string line;
74 
75  //read first line p3p
76  std::getline(pFile_, line);
77  if (!line.compare(""))
78  return false;
79 
80  //read second line p3p
81  Mdouble time;
82  unsigned N;
83  pFile_ >> time >> N;
84  dpm.setTime(time);
86  std::getline(pFile_, line);
87 
88  //read third line p3p
89  std::getline(pFile_, line);
90 
91  //read next lines p3p
92  logger(INFO,"Reading % particles at time %",N,time);
93  {
95  for (unsigned i = 0; i < N; ++i) {
96  unsigned id, species;
97  Mdouble volume, mass;
98  Vec3D position, velocity;
99  pFile_ >> id >> species >> volume >> mass >> position >> velocity;
100  std::getline(pFile_, line);
101  if (version_==Version::P3) species--;
102  //logger(INFO,"% % % % % %",id, species, volume, mass, position, velocity);
103  // add new species as necessary
104  while (dpm.speciesHandler.getNumberOfObjects() <= species) {
106  species.setDensity(mass / volume);
107  //use an infinite interaction radius
108  species.setAdhesionForceMax(1e20);
109  species.setAdhesionStiffness(1);
110  logger(INFO, "Adding species of density %", species.getDensity());
112  }
114  p.setRadius(cbrt(0.75 / pi * volume));
115  p.setPosition(position);
116  p.setVelocity(velocity);
117  dpm.particleHandler.copyAndAddObject(p)->setId(id);
118  }
119  }
120  //logger(INFO,"Read % particles",dpm.particleHandler.getNumberOfObjects());
121 
122  //read first line p3c
123  std::getline(cFile_, line);
124  if (!line.compare(""))
125  return false;
126 
127  //read second line p3c
128  cFile_ >> time >> N;
129  if (fabs(dpm.getTime()/time-1)>0.01)
130  logger(ERROR,"Timesteps in p3c and p3p do not agree: % %",dpm.getTime(),time);
132  std::getline(cFile_, line);
133 
134  //read third line p3c
135  std::getline(cFile_, line);
136 
137  //read next lines p3c
138  logger(INFO,"Reading % contacts",N,time);
139  for (unsigned i=0; i<N; ++i) {
140  unsigned id1, id2;
141  Vec3D force, contact;
142  cFile_ >> id1 >> id2;
145  logger.assert(p1!=nullptr,"Particle % does not exist",id1);
146  logger.assert(p2!=nullptr,"Particle % does not exist",id1);
147  Vec3D P1ToP2 = p2->getPosition()-p1->getPosition();
149  logger.assert(c!= nullptr,"Particle-particle interaction % % does not exist",p1,p2);
150  c->setDistance(P1ToP2.getLength());
151  c->setNormal(P1ToP2/c->getDistance());
152  c->setOverlap(c->getDistance()-p1->getRadius()-p2->getRadius());
153  if (version_==Version::P3) {
154  cFile_ >> force;
155  contact = p1->getPosition()-P1ToP2*((p1->getRadius()-0.5*c->getOverlap())/c->getDistance());
156  } else {
157  cFile_ >> contact >> force;
158  }
159  std::getline(cFile_, line);
160  c->setContactPoint(contact);
161  c->setForce(force);
162  if (i%(N/10)==0) {std::cout << "\r " << std::round((double)i/N*100) << '%'; std::cout.flush();}
163  }
164  std::cout << '\n';
165 
166  //read first line p3w
167  std::getline(wFile_, line);
168  if (!line.compare(""))
169  return false;
170 
171  //read second line p3w
172  wFile_ >> time >> N;
173  if (fabs(dpm.getTime()/time-1)>0.01)
174  logger(ERROR,"Timesteps in p3w and p3p do not agree");
176  std::getline(wFile_, line);
177 
178  //read third line p3w
179  std::getline(wFile_, line);
180 
181  //create wall
182  InfiniteWall wall;
184  auto w = dpm.wallHandler.copyAndAddObject(wall);
185 
186  //read next lines p3w
187  logger(INFO,"Reading % wall contacts",N,time);
188  for (unsigned i=0; i<N; ++i) {
189  unsigned id;
190  Vec3D force, contact, particleToContact;
191  wFile_ >> id;
193  logger.assert(p!=nullptr,"Particle % does not exist",id);
194  BaseInteraction* c = w->getInteractionWith(p,time,&dpm.interactionHandler);
195  logger.assert(c!= nullptr,"Particle-wall interaction % % does not exist",p,w);
196  if (version_==Version::P3) {
197  wFile_ >> force >> particleToContact;
198  contact = p->getPosition()-particleToContact;
199  } else {
200  wFile_ >> contact >> force;
201  particleToContact = p->getPosition()-contact;
202  }
203  std::getline(wFile_, line);
204  c->setContactPoint(contact);
205  c->setDistance(particleToContact.getLength());
206  c->setNormal(particleToContact/c->getDistance());
207  c->setOverlap(c->getDistance()-p->getRadius());
208  c->setForce(force);
209  }
210 
211  logger(INFO,"Writing output files");
212  for (const auto p : dpm.particleHandler) {
213  dpm.setMin(Vec3D::min(dpm.getMin(),p->getPosition()));
214  dpm.setMax(Vec3D::max(dpm.getMax(),p->getPosition()));
215  }
216  for (const auto c : dpm.interactionHandler) {
217  dpm.setMin(Vec3D::min(dpm.getMin(),c->getContactPoint()));
218  dpm.setMax(Vec3D::max(dpm.getMax(),c->getContactPoint()));
219  }
221  return true;
222  }
223 
228  pFile_.close();
229  wFile_.close();
230  cFile_.close();
231  logger(INFO,"Input files closed");
232  }
233 
234 private:
235 
237  std::ifstream pFile_, cFile_, wFile_;
239  enum class Version {P3, P4} version_;
240 
242 };
243 
244 int main(int argc, char** argv)
245 {
246  helpers::writeToFile("XRT","xrt");
247 
248  //Check to see if we actually received two arguments
249  if (argc < 2) {
250  //We didn't. Print a usage and exit the program.
251  logger(FATAL,"This program converts Particle Analytics (.p3* or p4*) to MercuryDPM files.\n"
252  //"These file can then be used in MercuryCG to analyse your data.\n"
253  "Usage: Call the executable with the base name as argument.\n"
254  "E.g. to convert name.p* call\n"
255  " ./ParticleAnalytics2Mercury name\n", argv[0]);
256  }
257  FileReader fileReader(argv[1]);
258  while (fileReader.read());
259 
260 }
void setTime(Mdouble time)
Sets a new value for the current simulation time.
Definition: DPMBase.cc:825
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: Helpers.cc:446
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
void setNormal(Vec3D normal)
Sets the normal vector between the two interacting objects.
A basic particle.
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
void setOverlap(Mdouble overlap)
Set the overlap between the two interacting object.
double Mdouble
Definition: GeneralDefine.h:34
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
void setMax(const Vec3D &max)
Sets the maximum coordinates of the problem domain.
Definition: DPMBase.cc:1043
BaseInteraction * getInteractionWith(BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler) override
Checks if particle is in interaction with given particle P, and if so, returns vector of pointer to t...
void setStorageCapacity(const unsigned int N)
Sets the storage capacity of this BaseHandler.
Definition: BaseHandler.h:669
T * getObjectById(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:565
Vec3D getMin() const
Definition: DPMBase.h:623
Version
The version number of the particle analytics files.
std::ifstream pFile_
Pointers for the input files.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void setContactPoint(Vec3D contactPoint)
Set the location of the contact point between the two interacting objects.
void setForce(Vec3D force)
set total force (this is used by the normal force, tangential forces are added use addForce) ...
void setSpecies(const ParticleSpecies *species)
void clear() override
Empties the whole ParticleHandler by removing all BaseParticle.
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
int main(int argc, char **argv)
void forceWriteOutputFiles()
Writes output files immediately, even if the current time step was not meant to be written...
Definition: DPMBase.cc:3825
~FileReader()
Destructor, simple closes the file.
static Vec3D min(const Vec3D &a, const Vec3D &b)
Calculates the pointwise minimum of two Vec3D.
Definition: Vector.cc:102
const Mdouble pi
Definition: ExtendedMath.h:45
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:36
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
void setMin(const Vec3D &min)
Sets the minimum coordinates of the problem domain.
Definition: DPMBase.cc:1079
void setDensity(Mdouble density)
Allows density_ to be changed.
This gives functionality to read information from binary formats like STL etc. This class is complete...
LL< Log::FATAL > FATAL
Fatal log level.
Definition: Logger.cc:52
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
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
virtual unsigned int getNumberOfObjects() const
Gets the number of real Object in this BaseHandler. (i.e. no mpi or periodic particles) ...
Definition: BaseHandler.h:648
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1319
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1359
Vec3D getMax() const
Definition: DPMBase.h:629
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1339
FileReader(std::string name)
Default constuction, requires to users to prove the name of the file that will be opened...
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:600
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:412
Mdouble getDistance() const
Returns an Mdouble which is the norm (length) of distance vector.
This is a class defining walls.
Definition: InfiniteWall.h:47
Mdouble getDensity() const
Allows density_ to be accessed.
Contains material and contact force properties.
Definition: Interaction.h:42
Definition: Vector.h:49
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:797
enum FileReader::Version version_
virtual void clear()
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
Definition: BaseHandler.h:528
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171
static Vec3D max(const Vec3D &a, const Vec3D &b)
Calculates the pointwise maximum of two Vec3D.
Definition: Vector.cc:89