ComponentNeBem2d.hh

Go to the documentation of this file.

#ifndef G_COMPONENT_NEBEM_2D_H
#define G_COMPONENT_NEBEM_2D_H
 
#include <array>
 
#include "Garfield/Component.hh"
 
namespace Garfield {

 
class ComponentNeBem2d : public Component {
 public:
  ComponentNeBem2d();
  ~ComponentNeBem2d() {}

  void SetMedium(Medium* medium) { m_medium = medium; }

  bool AddSegment(const double x0, const double y0, const double x1,
                  const double y1, const double v, const int ndiv = -1);
  bool AddWire(const double x, const double y, const double d, const double v,
               const int ntrap = 5);
  bool AddRegion(const std::vector<double>& xp, const std::vector<double>& yp,
                 Medium* medium, const unsigned int bctype = 4,
                 const double v = 0., const int ndiv = -1);
  void AddChargeDistribution(const double x, const double y, const double a,
                             const double b, const double rho);
  void SetRangeZ(const double zmin, const double zmax);

  bool Initialise();

  void SetNumberOfDivisions(const unsigned int ndiv);
 
  void SetNumberOfCollocationPoints(const unsigned int ncoll);
  void EnableAutoResizing(const bool on = true) { m_autoSize = on; }
  void EnableRandomCollocation(const bool on = true) {
    m_randomCollocation = on;
  }
  void SetMaxNumberOfIterations(const unsigned int niter);

  unsigned int GetNumberOfRegions() const { return m_regions.size(); }
  bool GetRegion(const unsigned int i, std::vector<double>& xv,
                 std::vector<double>& yv, Medium*& medium, unsigned int& bctype,
                 double& v);
  unsigned int GetNumberOfSegments() const { return m_segments.size(); }
  bool GetSegment(const unsigned int i, double& x0, double& y0, double& x1,
                  double& x2, double& v) const;
  unsigned int GetNumberOfWires() const { return m_wires.size(); }
  bool GetWire(const unsigned int i, double& x, double& y, double& d, double& v,
               double& q) const;
  size_t GetNumberOfElements() const override { return m_elements.size(); }
  bool GetElement(const unsigned int i, double& x0, double& y0, double& x1,
                  double& y1, double& q) const;
 
  Medium* GetMedium(const double x, const double y, const double z) override;
 
  void ElectricField(const double x, const double y, const double z, double& ex,
                     double& ey, double& ez, Medium*& m, int& status) override;
  void ElectricField(const double x, const double y, const double z, double& ex,
                     double& ey, double& ez, double& v, Medium*& m,
                     int& status) override;
  using Component::ElectricField;
  bool GetVoltageRange(double& vmin, double& vmax) override;
 
  bool GetBoundingBox(double& xmin, double& ymin, double& zmin, double& xmax,
                      double& ymax, double& zmax) override;
  bool GetElementaryCell(double& xmin, double& ymin, double& zmin, double& xmax,
                         double& ymax, double& zmax) override;
 
  bool CrossedWire(const double x0, const double y0, const double z0,
                   const double x1, const double y1, const double z1,
                   double& xc, double& yc, double& zc, const bool centre,
                   double& rc) override;
  bool InTrapRadius(const double q0, const double x0, const double y0,
                    const double z0, double& xw, double& yx,
                    double& rw) override;
 
 private:
  static const double InvEpsilon0;
  static const double InvTwoPiEpsilon0;

  unsigned int m_nDivisions = 5;
  unsigned int m_nCollocationPoints = 1;
  bool m_autoSize = false;
  bool m_randomCollocation = false;
  unsigned int m_nMaxIterations = 3;

  Medium* m_medium = nullptr;
  bool m_useRangeZ = false;
  double m_zmin = -1.;
  double m_zmax = 1.;
  enum BC {
    Voltage = 1,  
    Charge,       
    Floating,     
    Dielectric    
  };
  struct Region {
    std::vector<double> xv;    
    std::vector<double> yv;    
    Medium* medium;            
    std::pair<BC, double> bc;  
    unsigned int depth;        
    int ndiv;                  
  };

  std::vector<Region> m_regions;
 
  struct Segment {
    std::array<double, 2> x0;  
    std::array<double, 2> x1;  
    int region1;               
    int region2;               
    std::pair<BC, double> bc;  
    int ndiv;                  
  };

  std::vector<Segment> m_segments;
 
  struct Wire {
    double x, y;  
    double r;     
    double v;     
    double q;     
    int ntrap;    
  };

  std::vector<Wire> m_wires;
 
  struct Element {
    double x, y;  
    double a;     
    double cphi;  
    double sphi;  
    double q;     
    std::pair<BC, double> bc;  
    double lambda;             
  };

  std::vector<Element> m_elements;
 
  struct SpaceCharge {
    double x, y;  
    double a, b;  
    double q;     
    double v0;    
  };
  std::vector<SpaceCharge> m_spaceCharge;

  void EliminateOverlaps(std::vector<Segment>& segments);
  bool Discretise(const Segment& segment, std::vector<Element>& elements,
                  const double lambda, const unsigned int ndiv);
 
  bool ComputeInfluenceMatrix(std::vector<std::vector<double> >& infmat) const;
  bool InvertMatrix(std::vector<std::vector<double> >& influenceMatrix,
                    std::vector<std::vector<double> >& inverseMatrix) const;
  bool LUDecomposition(std::vector<std::vector<double> >& mat,
                       std::vector<int>& index) const;
  void LUSubstitution(const std::vector<std::vector<double> >& mat,
                      const std::vector<int>& index,
                      std::vector<double>& col) const;
 
  bool Solve(const std::vector<std::vector<double> >& inverseMatrix,
             const std::vector<double>& bc);
  bool CheckConvergence(const double tol, std::vector<bool>& ok);
  void SplitElement(Element& oldElement, std::vector<Element>& elements);

  double LinePotential(const double a, const double x, const double y) const;
  double WirePotential(const double r0, const double x, const double y) const;
  void LineField(const double a, const double x, const double y, double& ex,
                 double& ey) const;
  void WireField(const double r0, const double x, const double y, double& ex,
                 double& ey) const;

  double BoxPotential(const double a, const double b, const double x,
                      const double y, const double v0) const;
  void BoxField(const double a, const double b, const double x, const double y,
                double& ex, double& ey) const;
 
  void Reset() override;
  void UpdatePeriodicity() override;
  void ToLocal(const double xIn, const double yIn, const double cphi,
               const double sphi, double& xOut, double& yOut) const;
  void ToGlobal(const double xIn, const double yIn, const double cphi,
                const double sphi, double& xOut, double& yOut) const;

  int Field(const double x, const double y, const double z, double& ex,
            double& ey, double& ez, double& v, Medium*& m, const bool opt);
};
}  // namespace Garfield
 
#endif