ComponentFieldMap.hh

Go to the documentation of this file.

#if defined(__GPUCOMPILE__) || !defined(G_COMPONENT_FIELD_MAP_H)
 
#if !defined(__GPUCOMPILE__) && !defined(G_COMPONENT_FIELD_MAP_H)
#define G_COMPONENT_FIELD_MAP_H
#endif
 
#ifdef __GPUCOMPILE__
#else
 
#include <array>
#include <map>
#include <memory>
#include <string>
#include <vector>
 
#include "Garfield/Component.hh"
#include "Garfield/TetrahedralTree.hh"
#include "TMatrixD.h"
#include "TVectorD.h"
 
#endif
 
#ifndef __GPUCOMPILE__
 
namespace Garfield {

 
class ComponentFieldMap : public Component {
 public:
  ComponentFieldMap() = delete;
  ComponentFieldMap(const std::string& name);
  virtual ~ComponentFieldMap();

  bool Check();

  void PrintRange();

  void PrintMaterials();
  void DriftMedium(const size_t imat);
  void NotDriftMedium(const size_t imat);
  size_t GetNumberOfMaterials() const { return m_materials.size(); }
  double GetPermittivity(const size_t imat) const;
  double GetConductivity(const size_t imat) const;
  void SetMedium(const size_t imat, Medium* medium);
  Medium* GetMedium(const size_t imat) const;
  using Component::GetMedium;
  void SetGas(Medium* medium);
 
  size_t GetNumberOfElements() const override { return m_elements.size(); }
  bool GetElementNodes(const size_t i,
                       std::vector<size_t>& nodes) const override;
  bool GetElementRegion(const size_t i, size_t& mat,
                        bool& drift) const override;
  bool GetElement(const size_t i, double& vol, double& dmin,
                  double& dmax) const;
  size_t GetNumberOfNodes() const override { return m_nodes.size(); }
  bool GetNode(const size_t i, double& x, double& y, double& z) const override;
  double GetPotential(const size_t i) const;
 
  // Options
  void EnableCheckMapIndices(const bool on = true) {
    m_checkMultipleElement = on;
  }

  void EnableDeleteBackgroundElements(const bool on = true) {
    m_deleteBackground = on;
  }

  void EnableTetrahedralTreeForElementSearch(const bool on = true) {
    m_useTetrahedralTree = on;
  }

  void EnableConvergenceWarnings(const bool on = true) {
    m_printConvergenceWarnings = on;
  }
 
  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;
#endif
  __DEVICE__
  void WeightingField(const double x, const double y, const double z,
                      double& wx, double& wy, double& wz,
#ifdef __GPUCOMPILE__
                      size_t label);
#else
                    const std::string& label) override;
double WeightingPotential(const double x, const double y, const double z,
                          const std::string& label) override;
 
double DelayedWeightingPotential(double x, double y, double z, const double t,
                                 const std::string& label) override;
void DelayedWeightingPotentials(const double x, const double y, const double z,
                                const std::string& label,
                                std::vector<double>& dwp) override;
 
bool IsInBoundingBox(const double x, const double y, const double z) const {
  return x >= m_minBoundingBox[0] && x <= m_maxBoundingBox[0] &&
         y >= m_minBoundingBox[1] && y <= m_maxBoundingBox[1] &&
         z >= m_minBoundingBox[2] && y <= m_maxBoundingBox[2];
}
 
bool Is3d() override { return m_is3d; }
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;
 
std::map<std::string, std::vector<double> > GetWeightingPotentials() {
  return m_wpot;
}
 
bool GetVoltageRange(double& vmin, double& vmax) override {
  vmin = m_mapvmin;
  vmax = m_mapvmax;
  return true;
}

void CopyWeightingPotential(const std::string& label,
                            const std::string& labelSource, const double x,
                            const double y, const double z, const double alpha,
                            const double beta, const double gamma);

double CreateGPUTransferObject(ComponentGPU*& comp_gpu) override;
#endif
 protected:
  bool m_is3d = true;
 
  enum class ElementType {
    Unknown = 0,
    Serendipity = 5,
    CurvedTetrahedron = 13
  };
  ElementType m_elementType = ElementType::CurvedTetrahedron;
 
  // Elements
  struct Element {
    // Nodes
    int emap[10];
    // Material
    unsigned int matmap;
  };
#ifdef __GPUCOMPILE__
  Element* m_elements = nullptr;
  int* m_elementIndices = nullptr;
  int numElements = 0;
#else
std::vector<Element> m_elements;
std::vector<int> m_elementIndices;
#endif
 
// Degeneracy flags.
#ifdef __GPUCOMPILE__
  bool* m_degenerate;
#else
std::vector<bool> m_degenerate;
#endif
 
// Bounding boxes of the elements.
#ifdef __GPUCOMPILE__
  cuda_t** m_bbMin = nullptr;
  cuda_t** m_bbMax = nullptr;
#else
std::vector<std::array<double, 3> > m_bbMin;
std::vector<std::array<double, 3> > m_bbMax;
#endif
 
#ifdef __GPUCOMPILE__
  cuda_t*** m_w12 = nullptr;
#else
std::vector<std::array<std::array<double, 3>, 4> > m_w12;
#endif
 
  // Nodes
  struct Node {
    // Coordinates
    double x, y, z;
  };
#ifdef __GPUCOMPILE__
  Node* m_nodes = nullptr;
  int numNodes = 0;
#else
std::vector<Node> m_nodes;
#endif
 
// TODO GPU: m_dwpot not yet implemented
#ifndef __GPUCOMPILE__
  // Potentials.
  std::vector<double> m_pot;
  // Weighting potentials.
  std::map<std::string, std::vector<double> > m_wpot;
  // Delayed weighting potentials.
  std::map<std::string, std::vector<std::vector<double> > > m_dwpot;
#else
double* m_pot = nullptr;
int m_numpot = 0;
 
// The number of weighting potentials
int m_num_wpots = 0;
// The number of entries in each weighting potential
int* m_num_entries_wpot = nullptr;
double** m_wpot = nullptr;
#endif
 
  // Materials
  struct Material {
    // Permittivity
    double eps;
    // Resistivity
    double ohm;
    bool driftmedium;
#ifdef __GPUCOMPILE__
    MediumGPU* medium;
#else
  // Associated medium
  Medium* medium;
#endif
  };
 
#ifdef __GPUCOMPILE__
  Material* m_materials;
  int numMaterials;
#else
std::vector<Material> m_materials;
#endif
 
#ifndef __GPUCOMPILE__
  // Weighting potential copy
  struct WeightingFieldCopy {
    // Source
    std::string source;
    TMatrixD rot = TMatrixD(3, 3);
    TVectorD trans = TVectorD(3);
  };
 
  // Weighting potential copies.
  std::map<std::string, WeightingFieldCopy> m_wfieldCopies;
 
  // Bounding box
  bool m_hasBoundingBox = false;
  std::array<double, 3> m_minBoundingBox = {{0., 0., 0.}};
  std::array<double, 3> m_maxBoundingBox = {{0., 0., 0.}};

  bool m_cacheElemBoundingBoxes = false;
#endif
 
#ifdef __GPUCOMPILE__
  double m_mapmin[3];
  double m_mapmax[3];
  double m_mapamin[3];
  double m_mapamax[3];
#else
// Ranges and periodicities
std::array<double, 3> m_mapmin = {{0., 0., 0.}};
std::array<double, 3> m_mapmax = {{0., 0., 0.}};
std::array<double, 3> m_mapamin = {{0., 0., 0.}};
std::array<double, 3> m_mapamax = {{0., 0., 0.}};
#endif
 
#ifndef __GPUCOMPILE__
  std::array<double, 3> m_mapna = {{0., 0., 0.}};
  std::array<double, 3> m_cells = {{0., 0., 0.}};
 
  double m_mapvmin = 0.;
  double m_mapvmax = 0.;
 
  std::array<bool, 3> m_setang;
 
  // Option to delete meshing in conductors
  bool m_deleteBackground = true;
 
  // Warnings flag
  bool m_warning = false;
  unsigned int m_nWarnings = 0;
 
  // Print warnings about failed convergence when refining
  // isoparametric coordinates.
  bool m_printConvergenceWarnings = true;
 
  // Get the scaling factor for a given length unit.
  static double ScalingFactor(std::string unit);
 
  // Reset the component.
  void Reset() override;
 
  void Prepare();
 
  // Calculate x, y, z, V and angular ranges.
  virtual void SetRange();
 
  // Periodicities
  void UpdatePeriodicity() override {
    if (!m_is3d) UpdatePeriodicity2d();
    UpdatePeriodicityCommon();
  }
  void UpdatePeriodicity2d();
  void UpdatePeriodicityCommon();

  bool SetDefaultDriftMedium();
 
#endif
  __DEVICE__
  int Field(const double x, const double y, const double z, double& fx,
            double& fy, double& fz, int& iel,
#ifdef __GPUCOMPILE__
            const double* potentials, const int numPotentials
#else
          const std::vector<double>& potentials
#endif
  ) const;
#ifndef __GPUCOMPILE__
  double Potential(const double x, const double y, const double z,
                   const std::vector<double>& potentials) const;
  static double Potential3(const std::array<double, 6>& v,
                           const std::array<double, 3>& t);
  static void Field3(const std::array<double, 6>& v,
                     const std::array<double, 3>& t, double jac[4][4],
                     const double det, double& ex, double& ey);
  static double Potential5(const std::array<double, 8>& v,
                           const std::array<double, 2>& t);
  static void Field5(const std::array<double, 8>& v,
                     const std::array<double, 2>& t, double jac[4][4],
                     const double det, double& ex, double& ey);
  static double Potential13(const std::array<double, 10>& v,
                            const std::array<double, 4>& t);
#endif
  __DEVICE__ static void Field13(
#ifndef __GPUCOMPILE__
      const std::array<double, 10>& v, const std::array<double, 4>& t,
#else
    const double v[10], const double t[4],
#endif
      double jac[4][4], const double det, double& ex, double& ey, double& ez);
#ifndef __GPUCOMPILE__
  int FindElement5(const double x, const double y, double& t1, double& t2,
                   double& t3, double& t4, double jac[4][4], double& det) const;
#endif
  __DEVICE__
  int FindElement13(const double x, const double y, const double z, double& t1,
                    double& t2, double& t3, double& t4, double jac[4][4],
                    double& det) const;
#ifndef __GPUCOMPILE__
  int FindElementCube(const double x, const double y, const double z,
                      double& t1, double& t2, double& t3, TMatrixD*& jac,
                      std::vector<TMatrixD*>& dN) const;
#endif

  __DEVICE__
  void MapCoordinates(double& xpos, double& ypos, double& zpos, bool& xmirrored,
                      bool& ymirrored, bool& zmirrored, double& rcoordinate,
                      double& rotation) const;
  __DEVICE__
  void UnmapFields(double& ex, double& ey, double& ez, const double xpos,
                   const double ypos, const double zpos, const bool xmirrored,
                   const bool ymirrored, const bool zmirrored,
                   const double rcoordinate, const double rotation) const;
 
#ifndef __GPUCOMPILE__
  static int ReadInteger(char* token, int def, bool& error);
  static double ReadDouble(char* token, double def, bool& error);
 
  virtual double GetElementVolume(const size_t i) const;
  virtual void GetAspectRatio(const size_t i, double& dmin, double& dmax) const;
 
  void PrintWarning(const std::string& header);
  void PrintNotReady(const std::string& header) const;
  void PrintCouldNotOpen(const std::string& header,
                         const std::string& filename) const;
  void PrintElement(const std::string& header, const double x, const double y,
                    const double z, const double t1, const double t2,
                    const double t3, const double t4, const size_t i,
                    const std::vector<double>& potential) const;
  void TimeInterpolation(const double t, double& f0, double& f1, int& i0,
                         int& i1);
#endif
 
 protected:
  bool m_checkMultipleElement = false;
 
  // Tetrahedral tree
  bool m_useTetrahedralTree = true;
#ifdef __GPUCOMPILE__
  GARFIELD_CLASS_NAME(TetrahedralTree) * m_octree = nullptr;
#else
std::unique_ptr<GARFIELD_CLASS_NAME(TetrahedralTree)> m_octree;
#endif
 
 protected:
#ifndef __GPUCOMPILE__
  int Coordinates3(const double x, const double y, double& t1, double& t2,
                   double& t3, double& t4, double jac[4][4], double& det,
                   const std::array<double, 8>& xn,
                   const std::array<double, 8>& yn) const;
  int Coordinates4(const double x, const double y, double& t1, double& t2,
                   double& t3, double& t4, double& det,
                   const std::array<double, 8>& xn,
                   const std::array<double, 8>& yn) const;
  int Coordinates5(const double x, const double y, double& t1, double& t2,
                   double& t3, double& t4, double jac[4][4], double& det,
                   const std::array<double, 8>& xn,
                   const std::array<double, 8>& yn) const;
#endif

  __DEVICE__
  void Coordinates12(const double x, const double y, const double z, double& t1,
                     double& t2, double& t3, double& t4,
#ifdef __GPUCOMPILE__
                     const double xn[10], const double yn[10],
                     const double zn[10], const double w[4][3]
#else
                   const std::array<double, 10>& xn,
                   const std::array<double, 10>& yn,
                   const std::array<double, 10>& zn,
                   const std::array<std::array<double, 3>, 4>& w
#endif
  ) const;

  __DEVICE__
  int Coordinates13(const double x, const double y, const double z, double& t1,
                    double& t2, double& t3, double& t4, double jac[4][4],
                    double& det,
#ifdef __GPUCOMPILE__
                    const double xn[10], const double yn[10],
                    const double zn[10], cuda_t** w
#else
                  const std::array<double, 10>& xn,
                  const std::array<double, 10>& yn,
                  const std::array<double, 10>& zn,
                  const std::array<std::array<double, 3>, 4>& w
#endif
  ) const;
 
#ifndef __GPUCOMPILE__
  int CoordinatesCube(const double x, const double y, const double z,
                      double& t1, double& t2, double& t3, TMatrixD*& jac,
                      std::vector<TMatrixD*>& dN, const Element& element) const;

  static void Jacobian3(const std::array<double, 8>& xn,
                        const std::array<double, 8>& yn, const double u,
                        const double v, const double w, double& det,
                        double jac[4][4]);
  static void Jacobian5(const std::array<double, 8>& xn,
                        const std::array<double, 8>& yn, const double u,
                        const double v, double& det, double jac[4][4]);
#endif

  __DEVICE__
  static void Jacobian13(
#ifdef __GPUCOMPILE__
      const double xn[10], const double yn[10], const double zn[10],
#else
    const std::array<double, 10>& xn, const std::array<double, 10>& yn,
    const std::array<double, 10>& zn,
#endif
      const double fourt0, const double fourt1, const double fourt2,
      const double fourt3, double& det, double jac[4][4]);
 
#ifndef __GPUCOMPILE__
  void JacobianCube(const Element& element, const double t1, const double t2,
                    const double t3, TMatrixD*& jac,
                    std::vector<TMatrixD*>& dN) const;
 
  static std::array<std::array<double, 3>, 4> Weights12(
      const std::array<double, 10>& xn, const std::array<double, 10>& yn,
      const std::array<double, 10>& zn);

  void CalculateElementBoundingBoxes();

  bool InitializeTetrahedralTree();
};
}  // namespace Garfield
 
#endif
#endif