text: //Copyright (c) 2013-2021, The MercuryDPM Developers Team. All rights reserved.
  //For the list of developers, see <http://www.MercuryDPM.org/Team>.
  //
  //Redistribution and use in source and binary forms, with or without
  //modification, are permitted provided that the following conditions are met:
  //  * Redistributions of source code must retain the above copyright
  //    notice, this list of conditions and the following disclaimer.
  //  * Redistributions in binary form must reproduce the above copyright
  //    notice, this list of conditions and the following disclaimer in the
  //    documentation and/or other materials provided with the distribution.
  //  * Neither the name MercuryDPM nor the
  //    names of its contributors may be used to endorse or promote products
  //    derived from this software without specific prior written permission.
  //
  //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
  //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  #include "MercuryOS.h"
  #include <Walls/AxisymmetricIntersectionOfWalls.h>
  #include <Walls/TriangleWall.h>
  
  /**
   * Case Description:
   * A cylindrical drum is filled with bimodally distributed particles.
   * An initial placement of particles is read from the file Drum/PartCoordinates.txt. This placement is
   * generated by adding 8000 particles of material M2 with diameter of 4 mm, afterwards 30000 particles of
   * material M1 with diameter of 2 mm.
   * The initial placement of the walls is read from the file Drum/walls.txt. One can change to a smooth wall
   * implementation (non-triangulated) by setting useMercuryWalls_ to true.
   * The rotation velocity is 2 rad/s in the counterclockwise direction.
   * During simulation we measure the time-dependent change of number of particles M1 and M2 situated in zone
   * 1 (x>0 and y>0) and zone 2 (x<0 and y<0), where (x,y,z) denotes the particles' center of mass. This
   * information is outputted to the Drum.ene file.
   * The calculations are performed with a time step of 8*10 -7 s for a time period of 5 seconds.
   */
  class Drum : public MercuryOS
  {
      
      // used to set the initial conditions of the particles, walls, species, etc
      void setupInitialConditions() override
      {
          // name of the output files
          setName(soft() ? "DrumSoft" : "Drum");
          
          // turn on gravity
          setGravity({0, 0, -9.8});
          
          // set time step and maximum simulation time
          setTimeStep(soft() ? 1e-4 : 8e-7);
          setTimeMax(5.);
      
          // output frequency
          Mdouble outputPeriod = soft() ? 0.01 : 0.1;
          setSaveCount(static_cast<unsigned>(outputPeriod / getTimeStep()));
          
          // remove files from previous run
          removeOldFiles();
          
          // determine which output files to write
          if (writeOutput()) {
              // write more frequently
              setSaveCount(static_cast<unsigned>(0.01 / getTimeStep()));
              // write ene, data, fstat, restart and vtu files
              setParticlesWriteVTK(true);
              setWallsWriteVTK(FileType::MULTIPLE_FILES);
              fStatFile.writeFirstAndLastTimeStep();
              restartFile.writeFirstAndLastTimeStep();
          } else {
              // only .ene files are written
              setFileType(FileType::NO_FILE);
              eneFile.setFileType(FileType::ONE_FILE);
          }
      
          //use a shortened simulation in test mode
          if (test()) {
              setSaveCount(60);
              setTimeMax(600*getTimeStep());
              logger(INFO,"Test mode, reduced timeMax to %",getTimeMax());
          }
          
          // set domain for visualisation
          setMax(Vec3D(0.1, 0.03, 0.1));
          setMin(Vec3D(-0.1, -0.03, -0.1));
          
          // define the material properties of M1, M2, steel (see MercuryOS.h)
          setMaterialProperties();
          
          // read particle positions and radii from file
          {
              //open file
              std::string fileName = "Drum/PartCoordinates.txt";
              std::ifstream file(fileName);
              logger.assert_always(file.is_open(), "File % could not be opened", fileName);
              // read header line
              std::string header;
              getline(file, header);
              // read particle positions and radii
              SphericalParticle particle;
              double rad;
              Vec3D pos;
              while (file >> rad >> pos) {
                  particle.setSpecies(rad == 0.002 ? m2 : m1);
                  particle.setRadius(rad);
                  particle.setPosition(pos);
                  particleHandler.copyAndAddObject(particle);
              }
              logger(INFO, "Read % particles from %", particleHandler.getSize(), fileName);
          }
          
          // add walls
          if (useMercuryWalls()) {
              // use smooth walls
              InfiniteWall sideWall;
              sideWall.setSpecies(steel);
              sideWall.setAngularVelocity(Vec3D(0, -2, 0));
              sideWall.set(Vec3D(0, 1, 0), Vec3D(0, 3e-2, 0));
              wallHandler.copyAndAddObject(sideWall);
              sideWall.set(Vec3D(0, -1, 0), Vec3D(0, -3e-2, 0));
              wallHandler.copyAndAddObject(sideWall);
              
              AxisymmetricIntersectionOfWalls drumWall;
              drumWall.setSpecies(steel);
              drumWall.setAngularVelocity(Vec3D(0, -2, 0));
              drumWall.setAxis(Vec3D(0, 1, 0));
              drumWall.addObject(Vec3D(1, 0, 0), Vec3D(0.1, 0, 0));
              wallHandler.copyAndAddObject(drumWall);
              logger(INFO, "Created % walls", wallHandler.getSize());
          } else {
              // read triangle walls from file
              // open file
              std::string fileName = "Drum/walls.txt";
              std::ifstream file(fileName);
              logger.assert_always(file.is_open(), "File % could not be opened", fileName);
              // read header line
              std::string header;
              getline(file, header);
              // read triangle vertices
              TriangleWall wall;
              wall.setAngularVelocity(Vec3D(0, -2, 0));
              wall.setSpecies(steel);
              Vec3D a, b, c;
              while (file >> a >> b >> c) {
                  wall.setVertices(a, b, c, Vec3D(0, 0, 0));
                  wallHandler.copyAndAddObject(wall);
              }
              logger(INFO, "Read % walls from %", wallHandler.getSize(), fileName);
          }
      }
      
      // Write requested output to the ene file
      void writeEneTimeStep(std::ostream &os) const override
      {
          if (eneFile.getCounter() == 1) os << "Time zone2mm4 zone1mm4 zone2mm2 zone1mm2\n";
          // count the zones
          unsigned zone2mm2 = 0, zone1mm2 = 0, zone2mm4 = 0, zone1mm4 = 0;
          for (auto particle : particleHandler) {
              Vec3D pos = particle->getPosition();
              if (pos.X > 0 && pos.Z > 0) {
                  if (particle->getRadius() == 0.002) {
                      ++zone1mm4;
                  } else {
                      ++zone1mm2;
                  }
              } else if (pos.X < 0 && pos.Z < 0) {
                  if (particle->getRadius() == 0.002) {
                      ++zone2mm4;
                  } else {
                      ++zone2mm2;
                  }
              }
          }
          os << getTime() << ' '
             << zone2mm4 << ' '
             << zone1mm4 << ' '
             << zone2mm2 << ' '
             << zone1mm2 << std::endl;
      }
      
      // Also write the ene information to the screen
      void printTime() const override
      {
          writeEneTimeStep(std::cout);
      }
      
  };
  
  int main(int argc, char **argv)
  {
      // create an instance of the Drum class
      Drum dpm;
      // command line arguments:
      dpm.setNumberOfOMPThreads(helpers::readFromCommandLine(argc, argv, "-omp", 1));
      // turn on additional output files for viewing/analysing the data
      dpm.test(helpers::readFromCommandLine(argc, argv, "-test"));
      // turn on additional output files for viewing/analysing the data
      dpm.soft(helpers::readFromCommandLine(argc, argv, "-soft"));
      // turn on additional output files for viewing/analysing the data
      dpm.writeOutput(helpers::readFromCommandLine(argc, argv, "-writeOutput"));
      // use the command line argument -useMercuryWalls to turn on the smooth wall implementation
      dpm.useMercuryWalls(helpers::readFromCommandLine(argc, argv, "-useMercuryWalls"));
      // call the solve routine
      dpm.solve();
      return 0;
  }
  