ComponentNeBem3d.hh

Go to the documentation of this file.

#ifndef G_COMPONENT_NEBEM_3D_H
#define G_COMPONENT_NEBEM_3D_H
 
#include <array>
#include <map>
 
#include "Garfield/Component.hh"
#include "Garfield/Solid.hh"
 
namespace Garfield {

 
class ComponentNeBem3d : public Component {
 public:
  ComponentNeBem3d();
  ~ComponentNeBem3d() {}
 
  Medium* GetMedium(const double x, const double y, const double z) override;

  void AddPlaneX(const double x, const double voltage);
  void AddPlaneY(const double y, const double voltage);
  void AddPlaneZ(const double z, const double voltage);
  unsigned int GetNumberOfPlanesX() const;
  unsigned int GetNumberOfPlanesY() const;
  unsigned int GetNumberOfPlanesZ() const;
  bool GetPlaneX(const unsigned int i, double& x, double& v) const;
  bool GetPlaneY(const unsigned int i, double& y, double& v) const;
  bool GetPlaneZ(const unsigned int i, double& z, double& v) const;
 
  unsigned int GetNumberOfPrimitives() const { return m_primitives.size(); }
  bool GetPrimitive(const unsigned int i, double& a, double& b, double& c,
                    std::vector<double>& xv, std::vector<double>& yv,
                    std::vector<double>& zv, int& interface, double& v,
                    double& q, double& lambda) const;
  bool GetPrimitive(const unsigned int i, double& a, double& b, double& c,
                    std::vector<double>& xv, std::vector<double>& yv,
                    std::vector<double>& zv, int& vol1, int& vol2) const;
  bool GetVolume(const unsigned int vol, int& shape, int& material, double& eps,
                 double& potential, double& charge, int& bc);
  int GetVolume(const double x, const double y, const double z);
 
  size_t GetNumberOfElements() const override { return m_elements.size(); }
  bool GetElement(const unsigned int i, std::vector<double>& xv,
                  std::vector<double>& yv, std::vector<double>& zv,
                  int& interface, double& bc, double& lambda) const;

  bool Initialise();

  void SetTargetElementSize(const double length);
  void SetMinMaxNumberOfElements(const unsigned int nmin,
                                 const unsigned int nmax);
 
  void SetNewModel(const unsigned int NewModel);
  void SetNewMesh(const unsigned int NewMesh);
  void SetNewBC(const unsigned int NewBC);
  void SetNewPP(const unsigned int NewPP);
  void SetModelOptions(const unsigned int NewModel, const unsigned int NewMesh,
                       const unsigned int NewBC, const unsigned int NewPP);

  void SetStoreInflMatrix(const unsigned int OptStoreInflMatrix);
  void SetReadInflMatrix(const unsigned int OptReadInflMatrix);
  void SetStoreInvMatrix(const unsigned int OptStoreInvMatrix);
  void SetReadInvMatrix(const unsigned int OptReadInvMatrix);
  void SetStorePrimitives(const unsigned int OptStorePrimitives);
  void SetReadPrimitives(const unsigned int OptReadPrimitives);
  void SetStoreElements(const unsigned int OptStoreElements);
  void SetReadElements(const unsigned int OptReadElements);
  void SetFormattedFile(const unsigned int OptFormattedFile);
  void SetUnformattedFile(const unsigned int OptUnformattedFile);
  void SetStoreReadOptions(const unsigned int OptStoreInflMatrix,
                           const unsigned int OptReadInflMatrix,
                           const unsigned int OptStoreInvMatrix,
                           const unsigned int OptReadInvMatrix,
                           const unsigned int OptStorePrimitives,
                           const unsigned int OptReadPrimitives,
                           const unsigned int OptStoreElements,
                           const unsigned int OptReadElements,
                           const unsigned int OptFormattedFile,
                           const unsigned int OptUnformattedFile);
 
  // Set whether an older model is to be re-used. Expects an inverted matrix
  // stored during an earlier computation that had identical model and mesh.
  // This execution considers only a change in the boundary conditions.
  void SetReuseModel(void);

 
  // Functions that set computation details and constraints
  void SetSystemChargeZero(const unsigned int OptSystemChargeZero);
  void SetValidateSolution(const unsigned int OptValidateSolution);
  void SetForceValidation(const unsigned int OptForceValidation);
  void SetRepeatLHMatrix(const unsigned int OptRepeatLHMatrix);
  void SetComputeOptions(const unsigned int OptSystemChargeZero,
                         const unsigned int OptValidateSolution,
                         const unsigned int OptForceValidation,
                         const unsigned int OptRepeatLHMatrix);
 
  // Fast volume related information (physical potential and fields)
  void SetFastVolOptions(const unsigned int OptFastVol,
                         const unsigned int OptCreateFastPF,
                         const unsigned int OptReadFastPF);
  void SetFastVolVersion(const unsigned int VersionFV);
  void SetFastVolBlocks(const unsigned int NbBlocksFV);
 
  // Needs to include IdWtField information for each of these WtFld functions
  // Weighting potential and field related Fast volume information
  void SetWtFldFastVolOptions(const unsigned int IdWtField,
                              const unsigned int OptWtFldFastVol,
                              const unsigned int OptCreateWtFldFastPF,
                              const unsigned int OptReadWtFldFastPF);
  void SetWtFldFastVolVersion(const unsigned int IdWtField,
                              const unsigned int VersionWtFldFV);
  void SetWtFldFastVolBlocks(const unsigned int IdWtField,
                             const unsigned int NbBlocksWtFldFV);
 
  // Known charge options
  void SetKnownChargeOptions(const unsigned int OptKnownCharge);
 
  // Charging up options
  void SetChargingUpOptions(const unsigned int OptChargingUp);

  void UseLUInversion() { m_inversion = Inversion::LU; }
  void UseSVDInversion() { m_inversion = Inversion::SVD; }

  void SetPeriodicCopies(const unsigned int nx, const unsigned int ny,
                         const unsigned int nz);
  void GetPeriodicCopies(unsigned int& nx, unsigned int& ny,
                         unsigned int& nz) const {
    nx = m_nCopiesX;
    ny = m_nCopiesY;
    nz = m_nCopiesZ;
  }

  void SetPeriodicityX(const double s);
  void SetPeriodicityY(const double s);
  void SetPeriodicityZ(const double s);
  void SetMirrorPeriodicityX(const double s);
  void SetMirrorPeriodicityY(const double s);
  void SetMirrorPeriodicityZ(const double s);
  bool GetPeriodicityX(double& s) const;
  bool GetPeriodicityY(double& s) const;
  bool GetPeriodicityZ(double& s) const;

  void SetNumberOfThreads(const unsigned int n) { m_nThreads = n > 0 ? n : 1; }

  void SetPrimAfter(const int n) { m_primAfter = n; }

  void SetWtFldPrimAfter(const int n) { m_wtFldPrimAfter = n; }

  void SetOptRmPrim(const unsigned int n) { m_optRmPrim = n; }
 
  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;
 
  void WeightingField(const double x, const double y, const double z,
                      double& wx, double& wy, double& wz,
                      const std::string& label) override;
  double WeightingPotential(const double x, const double y, const double z,
                            const std::string& label) override;
 
 protected:
  void Reset() override;
  void UpdatePeriodicity() override;
 
 private:
  struct Primitive {
    double a, b, c;
    std::vector<double> xv;
    std::vector<double> yv;
    std::vector<double> zv;
    int interface = 0;
    double v = 0.;
    double q = 0.;
    double lambda = 0.;
    double elementSize;
    int vol1, vol2;
  };

  std::vector<Primitive> m_primitives;
 
  struct Element {
    std::array<double, 3> origin;
    double lx;
    double lz;
    double dA;
    std::array<std::array<double, 3>, 3> dcos;
    std::vector<double> xv;
    std::vector<double> yv;
    std::vector<double> zv;
    int interface = 0;
    double lambda = 0.;
    std::array<double, 3> collocationPoint;
    double bc = 0.;
    double assigned = 0.;
    double solution = 0.;
  };

  std::vector<Element> m_elements;

  std::array<bool, 6> m_ynplan{{false, false, false, false, false, false}};
  std::array<double, 6> m_coplan{{0., 0., 0., 0., 0., 0.}};
  std::array<double, 6> m_vtplan{{0., 0., 0., 0., 0., 0.}};
 
  // Model specifications
  unsigned int m_newModel = 1;
  unsigned int m_newMesh = 1;
  unsigned int m_newBC = 1;
  unsigned int m_newPP = 1;
 
  // Store and read options
  unsigned int m_optStoreInflMatrix = 0;
  unsigned int m_optReadInflMatrix = 0;
  unsigned int m_optStoreInvMatrix = 1;
  unsigned int m_optReadInvMatrix = 0;
  unsigned int m_optStorePrimitives = 0;
  unsigned int m_optReadPrimitives = 0;
  unsigned int m_optStoreElements = 0;
  unsigned int m_optReadElements = 0;
  unsigned int m_optStoreFormatted = 1;
  unsigned int m_optStoreUnformatted = 0;
 
  // Plot options
  // unsigned int m_optGnuplotPrimitives = 0;
  // unsigned int m_optGnuplotElements = 0;
  // unsigned int m_optPrimitiveFiles = 0;
  // unsigned int m_optElementFiles = 0;
 
  // Compute options
  unsigned int m_optSystemChargeZero = 1;
  unsigned int m_optValidateSolution = 1;
  unsigned int m_optForceValidation = 0;
  unsigned int m_optRepeatLHMatrix = 0;
 
  // Fast volume information (physical potential and fields)
  unsigned int m_optFastVol = 0;
  unsigned int m_optCreateFastPF = 0;
  unsigned int m_optReadFastPF = 0;
  unsigned int m_versionFV = 0;
  unsigned int m_nbBlocksFV = 0;
 
  // Weighting potential and field related Fast volume information
  unsigned int m_idWtField = 0;
  unsigned int m_optWtFldFastVol[11];
  unsigned int m_optCreateWtFldFastPF[11];
  unsigned int m_optReadWtFldFastPF[11];
  unsigned int m_versionWtFldFV[11];
  unsigned int m_nbBlocksWtFldFV[11];
 
  // Known charge options
  unsigned int m_optKnownCharge = 0;
 
  // Charging up options
  unsigned int m_optChargingUp = 0;
 
  // Number of threads to be used by neBEM.
  unsigned int m_nThreads = 1;
 
  // Number of repetitions, after which only primitive properties are used.
  // a negative value implies elements are used always.
  int m_primAfter = -1;
 
  // Number of repetitions, after which only primitive properties are used
  // for weighting field calculations.
  // A negative value implies only elements are used.
  int m_wtFldPrimAfter = -1;
 
  // Option for removing primitives from a device geometry.
  // Zero implies none to be removed.
  unsigned int m_optRmPrim = 0;
 
  static constexpr double MinDist = 1.e-6;
  double m_targetElementSize = 50.0e-4;
  unsigned int m_minNbElementsOnLength = 1;
  unsigned int m_maxNbElementsOnLength = 100;
  std::array<double, 3> m_periodicLength{{0., 0., 0.}};
  unsigned int m_nCopiesX = 5;
  unsigned int m_nCopiesY = 5;
  unsigned int m_nCopiesZ = 5;
 
  enum class Inversion { LU = 0, SVD };
  Inversion m_inversion = Inversion::LU;

  std::map<std::string, int> m_wfields;
 
  void InitValues();
  void ShiftPanels(std::vector<Panel>& panels) const;
  bool EliminateOverlaps(const Panel& panel1, const Panel& panel2,
                         std::vector<Panel>& panelsOut,
                         std::vector<int>& itypo);
 
  bool TraceEnclosed(const std::vector<double>& xl1,
                     const std::vector<double>& yl1,
                     const std::vector<double>& xl2,
                     const std::vector<double>& yl2, const Panel& originalPanel,
                     std::vector<Panel>& newPanels) const;
 
  void TraceNonOverlap(
      const std::vector<double>& xp1, const std::vector<double>& yp1,
      const std::vector<double>& xl1, const std::vector<double>& yl1,
      const std::vector<double>& xl2, const std::vector<double>& yl2,
      const std::vector<int>& flags1, const std::vector<int>& flags2,
      const std::vector<int>& links1, const std::vector<int>& links2,
      std::vector<bool>& mark1, int ip1, const Panel& originalPanel,
      std::vector<Panel>& newPanels) const;
 
  void TraceOverlap(
      const std::vector<double>& xp1, const std::vector<double>& yp1,
      const std::vector<double>& xp2, const std::vector<double>& yp2,
      const std::vector<double>& xl1, const std::vector<double>& yl1,
      const std::vector<double>& xl2, const std::vector<double>& yl2,
      const std::vector<int>& flags1, const std::vector<int>& links1,
      const std::vector<int>& links2, std::vector<bool>& mark1, int ip1,
      int ip2, const Panel& originalPanel, std::vector<Panel>& newPanels) const;

  bool MakePrimitives(const Panel& panelIn,
                      std::vector<Panel>& panelsOut) const;

  bool SplitTrapezium(const Panel panelIn, std::vector<Panel>& stack,
                      std::vector<Panel>& panelsOut, const double epsang) const;
 
  unsigned int NbOfSegments(const double length, const double target) const;
  bool DiscretizeWire(const Primitive& primitive, const double targetSize,
                      std::vector<Element>& elements) const;
  bool DiscretizeTriangle(const Primitive& primitive, const double targetSize,
                          std::vector<Element>& elements) const;
  bool DiscretizeRectangle(const Primitive& prim, const double targetSize,
                           std::vector<Element>& elements) const;
  int InterfaceType(const Solid::BoundaryCondition bc) const;
};
 
extern ComponentNeBem3d* gComponentNeBem3d;
 
}  // namespace Garfield
 
#endif