example-002-dynamics.cpp

#include <elec.hpp>
#include <vector>
#include <iterator>
#include <cmath>
#include <cstdlib>

// This splits the electrons into subsets. A single simulation step
// consist of SIM_SPLIT repetitions of the following : compute E, move
// the ith subset. The higher SIM_SPLIT, the smoother particles
// motion.... but the higher time consumption.

#define SIM_SPLIT     10
#define SIM_DT         1

#define MAX_PLOT_FIELD_NORM  100000
#define PLOT_FIELD_SCALE     5e-6
#define MESH_SIZE            40

#define NB_V_CONTOURS 20

#define RESISTIVITY 1 // Good conductor

#define ROTATION_VIEW .1

int main(int argc, char* argv[]) {

  if(argc != 4) {
    std::cout << "Usage : " << argv[0] << " <particle|field|both> <a|b|..> <nb-frames> | elec-plot.py" << std::endl
              << "        " << argv[0] << " <particle|field|both> <a|b|..> <nb-frames> | elec-plotV.py" << std::endl;
      
    return 0;
  }

  bool part  = true;
  bool field = true;
  std::string mode(argv[1]);
  std::string cond(argv[2]);
  unsigned int nb_steps = (unsigned int)(atoi(argv[3]));
  double vmin,vmax;
  double a = 20;
  
  if(mode == "particle") {
    part  = true;
    field = false;
  }
  else if(mode == "field") {
    part  = false;
    field = true;
  }

  // The world handles the simulation, i.e. its particles, the
  // conductors, as well as a generator (not used here).

  elec::World world;
  elec::conductor::Disk      da1({-.5 ,    0},          .20, RESISTIVITY);
  elec::conductor::Disk      da2({ .5 ,    0},          .30, RESISTIVITY);
  elec::conductor::Rectangle ra1({-.5 , -.14},  {.5 ,  .15}, RESISTIVITY);
  elec::conductor::Disk      db1({-.2 , -.2 },           .4, RESISTIVITY);
  elec::conductor::Rectangle rb1({-.25,    0}, {-.15,  .6 }, RESISTIVITY);
  elec::conductor::Rectangle rb2({0   , -.25}, {.6  , -.15}, RESISTIVITY);
  elec::conductor::Disk      d  ({0   ,    0},           .5, RESISTIVITY);
  
  elec::Point rand_min, rand_max;
  bool uniform_electron_init = false;
  double electron_ratio = 1;

  // cond values allows to generate different conductor
  // configurations.
  
  if(cond == "a") {

    // First, add conductors in the world. Be careful, in case of
    // overlaps, the first added conductor defines the space
    // properties at overlapping areas.

    world += da1;
    world += da2;
    world += ra1;
    rand_min = {-.05,-.05};
    rand_max = { .05, .05};
    electron_ratio = 3;
    vmin = -10000;
    vmax =  0;
  }
  else if(cond == "b") {
    world += db1;
    world += rb1;
    world += rb2;
    uniform_electron_init = true;
    electron_ratio = 3;
    vmin = -10000;
    vmax =  0;
  }
  else if(cond == "c") {
    world += da1;
    world += da2;
    world += ra1;
    rand_min = {-.05,-.05};
    rand_max = { .05, .05};
    electron_ratio = .33;
    vmin = -1000;
    vmax =  2000;
  }
  else {
    world += d;
    rand_min = {-.05,-.05};
    rand_max = { .05, .05};
    electron_ratio = 3;
    vmin = -10000;
    vmax =  0;
  }


  // Now, we have to set up protons. They will not move in the
  // simulation.
  world.add_protons();

  // Let us set the same number of electrons in case of
  // uniform_electron_init. This correspond to an electrically neutral
  // piece of metal.
  unsigned int nb_electrons = (unsigned int)(world.nb_protons()*electron_ratio);
  if(!uniform_electron_init) 
    for(unsigned int i = 0; i < nb_electrons; ++i)
      world.add_electron(elec::uniform(rand_min,rand_max));
  else
    world.add_electrons(3*world.nb_protons()); 

  // Notify the end of particle definition each time yu add particles. 
  world.end_of_particles();
  

  elec::Point min = {-1,-1};
  elec::Point max = { 1, 1};

  std::vector<elec::Point> mesh;
  std::vector<elec::Point> Efield;
  std::vector<double>      Vfield;

  double                   max_norm;
  std::pair<double,double> vbound;

  elec::mesh(min,max,MESH_SIZE,MESH_SIZE,std::back_inserter(mesh));

  auto plot = elec::plot::figure("Conductor", min, max);
  
  Efield.clear();
  Vfield.clear();
  if(field) {
    max_norm = elec::E(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
    vbound = elec::V(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
  }

  for(unsigned int step = 0; step < nb_steps; ++step, a+=ROTATION_VIEW) {
    
    if(field)
      std::cerr << "max |E| = " << max_norm << std::endl
                << "V in [" << vbound.first << ", " << vbound.second << ']' << std::endl;

    plot.begin_frame("frame","jpg"); 
    plot.set_view(50,a);
    if(field) {
      plot.__plot_potential(mesh.begin(), mesh.end(), Vfield.begin(), Vfield.end(), vmin, vmax, NB_V_CONTOURS);
      plot.__plot_field(mesh.begin(), mesh.end(), Efield.begin(), Efield.end(),PLOT_FIELD_SCALE,MAX_PLOT_FIELD_NORM);
    }

    if(part) 
      plot.__plot_particles(world.particles().begin(), world.particles().end());
    plot.end_frame();
    
    world.sim(SIM_DT,SIM_SPLIT);
    if(field) {
      Efield.clear();
      Vfield.clear();
      max_norm = elec::E(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
      vbound = elec::V(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
    }
  }

  return 0;
}