Base class for components based on finite-element field maps. More...

#include <ComponentFieldMap.hh>

Inheritance diagram for Garfield::ComponentFieldMap:

Classes
struct	Element

struct	Material

struct	Node

Public Member Functions
	ComponentFieldMap ()
	Constructor.

virtual	~ComponentFieldMap ()
	Destructor.

virtual void	SetRange ()
	Calculate x, y, z, V and angular ranges.

void	PrintRange ()
	Show x, y, z, V and angular ranges.

bool	IsInBoundingBox (const double x, const double y, const double z) const

bool	GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
	Get the bounding box coordinates.

bool	GetVoltageRange (double &vmin, double &vmax) override
	Calculate the voltage range [V].

void	PrintMaterials ()
	List all currently defined materials.

void	DriftMedium (const unsigned int imat)
	Flag a field map material as a drift medium.

void	NotDriftMedium (const unsigned int imat)
	Flag a field map materials as a non-drift medium.

unsigned int	GetNumberOfMaterials () const
	Return the number of materials in the field map.

double	GetPermittivity (const unsigned int imat) const
	Return the permittivity of a field map material.

double	GetConductivity (const unsigned int imat) const
	Return the conductivity of a field map material.

void	SetMedium (const unsigned int imat, Medium *medium)
	Associate a field map material with a Medium class.

Medium *	GetMedium (const unsigned int i) const
	Return the Medium associated to a field map material.

unsigned int	GetNumberOfMedia () const

int	GetNumberOfElements () const
	Return the number of mesh elements.

bool	GetElement (const unsigned int i, double &vol, double &dmin, double &dmax)
	Return the volume and aspect ratio of a mesh element.

void	EnableCheckMapIndices ()

void	DisableCheckMapIndices ()

void	EnableDeleteBackgroundElements (const bool on=true)
	Option to eliminate mesh elements in conductors (default: on).

void	EnableTetrahedralTreeForElementSearch (const bool on=true)

virtual Medium *	GetMedium (const double x, const double y, const double z)
	Get the medium at a given location (x, y, z).

Public Member Functions inherited from Garfield::ComponentBase
	ComponentBase ()
	Constructor.

virtual	~ComponentBase ()
	Destructor.

virtual void	SetGeometry (GeometryBase *geo)
	Define the geometry.

virtual void	Clear ()
	Reset.

virtual Medium *	GetMedium (const double x, const double y, const double z)
	Get the medium at a given location (x, y, z).

virtual void	ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0

virtual void	ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
	Calculate the drift field [V/cm] and potential [V] at (x, y, z).

virtual bool	GetVoltageRange (double &vmin, double &vmax)=0
	Calculate the voltage range [V].

virtual void	WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)

virtual double	WeightingPotential (const double x, const double y, const double z, const std::string &label)

virtual void	DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)

virtual void	MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)

void	SetMagneticField (const double bx, const double by, const double bz)
	Set a constant magnetic field.

virtual bool	IsReady ()
	Ready for use?

virtual bool	GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
	Get the bounding box coordinates.

double	IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)

double	IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)

double	IntegrateFlux (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)

virtual bool	IsWireCrossed (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)

virtual bool	IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)

void	EnablePeriodicityX (const bool on=true)
	Enable simple periodicity in the direction.

void	DisablePeriodicityX ()

void	EnablePeriodicityY (const bool on=true)
	Enable simple periodicity in the direction.

void	DisablePeriodicityY ()

void	EnablePeriodicityZ (const bool on=true)
	Enable simple periodicity in the direction.

void	DisablePeriodicityZ ()

void	EnableMirrorPeriodicityX (const bool on=true)
	Enable mirror periodicity in the direction.

void	DisableMirrorPeriodicityX ()

void	EnableMirrorPeriodicityY (const bool on=true)
	Enable mirror periodicity in the direction.

void	DisableMirrorPeriodicityY ()

void	EnableMirrorPeriodicityZ (const bool on=true)
	Enable mirror periodicity in the direction.

void	DisableMirrorPeriodicityZ ()

void	EnableAxialPeriodicityX (const bool on=true)
	Enable axial periodicity in the direction.

void	DisableAxialPeriodicityX ()

void	EnableAxialPeriodicityY (const bool on=true)
	Enable axial periodicity in the direction.

void	DisableAxialPeriodicityY ()

void	EnableAxialPeriodicityZ (const bool on=true)
	Enable axial periodicity in the direction.

void	DisableAxialPeriodicityZ ()

void	EnableRotationSymmetryX (const bool on=true)
	Enable rotation symmetry around the axis.

void	DisableRotationSymmetryX ()

void	EnableRotationSymmetryY (const bool on=true)
	Enable rotation symmetry around the axis.

void	DisableRotationSymmetryY ()

void	EnableRotationSymmetryZ (const bool on=true)
	Enable rotation symmetry around the axis.

void	DisableRotationSymmetryZ ()

void	EnableDebugging ()
	Switch on debugging messages.

void	DisableDebugging ()
	Switch off debugging messages.

void	ActivateTraps ()
	Request trapping to be taken care of by the component (for TCAD).

void	DeactivateTraps ()

bool	IsTrapActive ()

void	ActivateVelocityMap ()
	Request velocity to be taken care of by the component (for TCAD).

void	DectivateVelocityMap ()

bool	IsVelocityActive ()

virtual bool	ElectronAttachment (const double, const double, const double, double &eta)
	Get the electron attachment coefficient.

virtual bool	HoleAttachment (const double, const double, const double, double &eta)
	Get the hole attachment coefficient.

virtual void	ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
	Get the electron drift velocity.

virtual void	HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
	Get the hole drift velocity.

virtual bool	GetElectronLifetime (const double, const double, const double, double &etau)

virtual bool	GetHoleLifetime (const double, const double, const double, double &htau)

Protected Member Functions
void	Reset () override
	Reset the component.

void	UpdatePeriodicity2d ()

void	UpdatePeriodicityCommon ()

int	FindElement5 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
	Find the element for a point in curved quadratic quadrilaterals.

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)
	Find the element for a point in curved quadratic tetrahedra.

int	FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD &jac, std::vector< TMatrixD > &dN)
	Find the element for a point in a cube.

void	MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
	Move (xpos, ypos, zpos) to field map coordinates.

void	UnmapFields (double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
	Move (ex, ey, ez) to global coordinates.

int	ReadInteger (char *token, int def, bool &error)

double	ReadDouble (char *token, double def, bool &error)

virtual double	GetElementVolume (const unsigned int i)=0

virtual void	GetAspectRatio (const unsigned int i, double &dmin, double &dmax)=0

void	PrintWarning (const std::string &header)

void	PrintNotReady (const std::string &header) 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 Element &element, const unsigned int n, const int iw=-1) const

Protected Member Functions inherited from Garfield::ComponentBase
virtual void	Reset ()=0
	Reset the component.

virtual void	UpdatePeriodicity ()=0
	Verify periodicities.

Protected Attributes
bool	m_is3d = true

int	nElements = -1

std::vector< Element >	elements

int	nNodes = -1

std::vector< Node >	nodes

unsigned int	m_nMaterials = 0

std::vector< Material >	materials

int	nWeightingFields = 0

std::vector< std::string >	wfields

std::vector< bool >	wfieldsOk

bool	hasBoundingBox = false

std::array< double, 3 >	m_minBoundingBox

std::array< double, 3 >	m_maxBoundingBox

std::array< double, 3 >	m_mapmin

std::array< double, 3 >	m_mapmax

std::array< double, 3 >	m_mapamin

std::array< double, 3 >	m_mapamax

std::array< double, 3 >	m_mapna

std::array< double, 3 >	m_cells

double	m_mapvmin

double	m_mapvmax

std::array< bool, 3 >	m_setang

bool	m_deleteBackground = true

bool	m_warning = false

unsigned int	m_nWarnings = 0

Protected Attributes inherited from Garfield::ComponentBase
std::string	m_className = "ComponentBase"
	Class name.

GeometryBase *	m_geometry = nullptr
	Pointer to the geometry.

bool	m_ready = false
	Ready for use?

bool	m_activeTraps = false
	Does the component have traps?

bool	m_hasVelocityMap = false
	Does the component have velocity maps?

std::array< bool, 3 >	m_periodic = {{false, false, false}}
	Simple periodicity in x, y, z.

std::array< bool, 3 >	m_mirrorPeriodic = {{false, false, false}}
	Mirror periodicity in x, y, z.

std::array< bool, 3 >	m_axiallyPeriodic = {{false, false, false}}
	Axial periodicity in x, y, z.

std::array< bool, 3 >	m_rotationSymmetric = {{false, false, false}}
	Rotation symmetry around x-axis, y-axis, z-axis.

double	m_bx0 = 0.

double	m_by0 = 0.

double	m_bz0 = 0.

bool	m_debug = false
	Switch on/off debugging messages.

Friends
class	ViewFEMesh

Detailed Description

Base class for components based on finite-element field maps.

Definition at line 13 of file ComponentFieldMap.hh.

Constructor & Destructor Documentation

◆ ComponentFieldMap()

Garfield::ComponentFieldMap::ComponentFieldMap ( )

Constructor.

Definition at line 14 of file ComponentFieldMap.cc.

                                     : ComponentBase() {
  m_className = "ComponentFieldMap";
}

◆ ~ComponentFieldMap()

Garfield::ComponentFieldMap::~ComponentFieldMap ( )

virtual

Destructor.

Definition at line 18 of file ComponentFieldMap.cc.

                                      {
  if (m_tetTree) delete m_tetTree;
}

Member Function Documentation

◆ DisableCheckMapIndices()

void Garfield::ComponentFieldMap::DisableCheckMapIndices ( )

inline

Definition at line 72 of file ComponentFieldMap.hh.

72{ m_checkMultipleElement = false; }

◆ DriftMedium()

void Garfield::ComponentFieldMap::DriftMedium ( const unsigned int imat )

Flag a field map material as a drift medium.

Definition at line 51 of file ComponentFieldMap.cc.

                                                           {
  // Do not proceed if not properly initialised.
  if (!m_ready) PrintNotReady("DriftMedium");
 
  // Check value
  if (imat >= m_nMaterials) {
    std::cerr << m_className << "::DriftMedium: Index out of range.\n";
    return;
  }
 
  // Make drift medium
  materials[imat].driftmedium = true;
}

◆ EnableCheckMapIndices()

void Garfield::ComponentFieldMap::EnableCheckMapIndices ( )

inline

Definition at line 68 of file ComponentFieldMap.hh.

                               {
    m_checkMultipleElement = true;
    m_lastElement = -1;
  }

◆ EnableDeleteBackgroundElements()

void Garfield::ComponentFieldMap::EnableDeleteBackgroundElements ( const bool on = true )

inline

Option to eliminate mesh elements in conductors (default: on).

Definition at line 74 of file ComponentFieldMap.hh.

                                                            { 
    m_deleteBackground = on; 
  }

◆ EnableTetrahedralTreeForElementSearch()

void Garfield::ComponentFieldMap::EnableTetrahedralTreeForElementSearch ( const bool on = true )

inline

Enable or disable the usage of the tetrahedral tree for searching the element in the mesh.

Definition at line 80 of file ComponentFieldMap.hh.

                                                                   {
    m_useTetrahedralTree = on;
  }

◆ FindElement13()

int Garfield::ComponentFieldMap::FindElement13	(	const double	x,
		const double	y,
		const double	z,
		double &	t1,
		double &	t2,
		double &	t3,
		double &	t4,
		double	jac[4][4],
		double &	det
	)

protected

Find the element for a point in curved quadratic tetrahedra.

Definition at line 298 of file ComponentFieldMap.cc.

                                                  {
  // Check if bounding boxes of elements have been computed
  if (!m_cacheElemBoundingBoxes) {
    std::cout << m_className << "::FindElement13:\n"
              << "    Caching the bounding boxes of all elements...";
    CalculateElementBoundingBoxes();
    std::cout << " done.\n";
    m_cacheElemBoundingBoxes = true;
  }
 
  // Backup
  double jacbak[4][4];
  double detbak = 1.;
  double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
  int imapbak = -1;
 
  // Initial values.
  t1 = t2 = t3 = t4 = 0.;
 
  // Check previously used element
  if (m_lastElement > -1 && !m_checkMultipleElement) {
    const Element& element = elements[m_lastElement];
    if (Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, element) == 0) {
      if (t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 && t3 <= +1 &&
          t4 >= 0 && t4 <= +1) {
        return m_lastElement;
      }
    }
  }
 
  // Tetra list in the block that contains the input 3D point.
  std::vector<int> tetList;
  if (m_useTetrahedralTree) {
    if (!m_isTreeInitialized) {
      if (!InitializeTetrahedralTree()) {
        std::cerr << m_className << "::FindElement13:\n";
        std::cerr << "    Tetrahedral tree initialization failed.\n";
        return -1;
      }
    }
    tetList = m_tetTree->GetTetListInBlock(Vec3(x, y, z));
  }
  // Number of elements to scan.
  // With tetra tree disabled, all elements are scanned.
  const int numElemToSearch = m_useTetrahedralTree ? tetList.size() : nElements;
  // Verify the count of volumes that contain the point.
  int nfound = 0;
  int imap = -1;
 
  // Scan all elements
  for (int i = 0; i < numElemToSearch; i++) {
    const int idxToElemList = m_useTetrahedralTree ? tetList[i] : i;
    const Element& element = elements[idxToElemList];
    if (x < element.xmin || x > element.xmax || y < element.ymin ||
        y > element.ymax || z < element.zmin || z > element.zmax)
      continue;
    if (Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, element) != 0) {
      continue;
    }
    if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1 || t4 < 0 ||
        t4 > 1) {
      continue;
    }
    ++nfound;
    imap = idxToElemList;
    m_lastElement = idxToElemList;
    if (m_debug) {
      std::cout << m_className << "::FindElement13:\n";
      std::cout << "    Found matching element " << i << ".\n";
    }
    if (!m_checkMultipleElement) return idxToElemList;
    for (int j = 0; j < 4; ++j) {
      for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
    }
    detbak = det;
    t1bak = t1;
    t2bak = t2;
    t3bak = t3;
    t4bak = t4;
    imapbak = imap;
    if (m_debug) {
      PrintElement("FindElement13", x, y, z, t1, t2, t3, t4, element, 10);
    }
  }
 
  // In checking mode, verify the tetrahedron/triangle count.
  if (m_checkMultipleElement) {
    if (nfound < 1) {
      if (m_debug) {
        std::cout << m_className << "::FindElement13:\n";
        std::cout << "    No element matching point (" << x << ", " << y << ", "
                  << z << ") found.\n";
      }
      m_lastElement = -1;
      return -1;
    }
    if (nfound > 1) {
      std::cerr << m_className << "::FindElement13:\n";
      std::cerr << "    Found << " << nfound << " elements matching point ("
                << x << ", " << y << ", " << z << ").\n";
    }
    if (nfound > 0) {
      for (int j = 0; j < 4; ++j) {
        for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
      }
      det = detbak;
      t1 = t1bak;
      t2 = t2bak;
      t3 = t3bak;
      t4 = t4bak;
      imap = imapbak;
      m_lastElement = imap;
      return imap;
    }
  }
 
  if (m_debug) {
    std::cout << m_className << "::FindElement13:\n";
    std::cout << "    No element matching point (" << x << ", " << y << ", "
              << z << ") found.\n";
  }
  return -1;
}

◆ FindElement5()

int Garfield::ComponentFieldMap::FindElement5	(	const double	x,
		const double	y,
		const double	z,
		double &	t1,
		double &	t2,
		double &	t3,
		double &	t4,
		double	jac[4][4],
		double &	det
	)

protected

Find the element for a point in curved quadratic quadrilaterals.

Definition at line 140 of file ComponentFieldMap.cc.

                                                 {
  // Check if bounding boxes of elements have been computed
  if (!m_cacheElemBoundingBoxes) {
    std::cout << m_className << "::FindElement5:\n"
              << "    Caching the bounding boxes of all elements...";
    CalculateElementBoundingBoxes();
    std::cout << " done.\n";
    m_cacheElemBoundingBoxes = true;
  }
 
  // Tetra list in the block that contains the input 3D point.
  std::vector<int> tetList;
  if (m_useTetrahedralTree) {
    if (!m_isTreeInitialized) {
      if (!InitializeTetrahedralTree()) {
        std::cerr << m_className << "::FindElement5:\n";
        std::cerr << "    Tetrahedral tree initialization failed.\n";
        return -1;
      }
    }
    tetList = m_tetTree->GetTetListInBlock(Vec3(x, y, z));
  }
  // Backup
  double jacbak[4][4], detbak = 1.;
  double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
  int imapbak = -1;
 
  // Initial values.
  t1 = t2 = t3 = t4 = 0;
 
  // Check previously used element
  if (m_lastElement > -1 && !m_checkMultipleElement) {
    const Element& element = elements[m_lastElement];
    if (element.degenerate) {
      if (Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, element) == 0) {
        if (t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 && t3 <= +1) {
          return m_lastElement;
        }
      }
    } else {
      if (Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, element) == 0) {
        if (t1 >= -1 && t1 <= +1 && t2 >= -1 && t2 <= +1) return m_lastElement;
      }
    }
  }
 
  // Verify the count of volumes that contain the point.
  int nfound = 0;
  int imap = -1;
 
  // Number of elements to scan.
  // With tetra tree disabled, all elements are scanned.
  const int numElemToSearch = m_useTetrahedralTree ? tetList.size() : nElements;
  for (int i = 0; i < numElemToSearch; ++i) {
    const int idxToElemList = m_useTetrahedralTree ? tetList[i] : i;
    const Element& element = elements[idxToElemList];
    if (x < element.xmin || x > element.xmax || y < element.ymin ||
        y > element.ymax || z < element.zmin || z > element.zmax)
      continue;
    if (element.degenerate) {
      // Degenerate element
      if (Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, element) != 0) {
        continue;
      }
      if (t1 < 0 || t1 > 1 || t2 < 0 || t2 > 1 || t3 < 0 || t3 > 1) continue;
      ++nfound;
      imap = idxToElemList;
      m_lastElement = idxToElemList;
      if (m_debug) {
        std::cout << m_className << "::FindElement5:\n";
        std::cout << "    Found matching degenerate element " << idxToElemList
                  << ".\n";
      }
      if (!m_checkMultipleElement) return idxToElemList;
      for (int j = 0; j < 4; ++j) {
        for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
      }
      detbak = det;
      t1bak = t1;
      t2bak = t2;
      t3bak = t3;
      t4bak = t4;
      imapbak = imap;
      if (m_debug) {
        PrintElement("FindElement5", x, y, z, t1, t2, t3, t4, element, 6);
      }
    } else {
      // Non-degenerate element
      if (Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, element) != 0) {
        continue;
      }
      if (t1 < -1 || t1 > 1 || t2 < -1 || t2 > 1) continue;
      ++nfound;
      imap = idxToElemList;
      m_lastElement = idxToElemList;
      if (m_debug) {
        std::cout << m_className << "::FindElement5:\n";
        std::cout << "    Found matching non-degenerate element "
                  << idxToElemList << ".\n";
      }
      if (!m_checkMultipleElement) return idxToElemList;
      for (int j = 0; j < 4; ++j) {
        for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
      }
      detbak = det;
      t1bak = t1;
      t2bak = t2;
      t3bak = t3;
      t4bak = t4;
      imapbak = imap;
      if (m_debug) {
        PrintElement("FindElement5", x, y, z, t1, t2, t3, t4, element, 8);
      }
    }
  }
 
  // In checking mode, verify the tetrahedron/triangle count.
  if (m_checkMultipleElement) {
    if (nfound < 1) {
      if (m_debug) {
        std::cout << m_className << "::FindElement5:\n";
        std::cout << "    No element matching point (" << x << ", " << y
                  << ") found.\n";
      }
      m_lastElement = -1;
      return -1;
    }
    if (nfound > 1) {
      std::cout << m_className << "::FindElement5:\n";
      std::cout << "    Found " << nfound << " elements matching point (" << x
                << ", " << y << ").\n";
    }
    if (nfound > 0) {
      for (int j = 0; j < 4; ++j) {
        for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
      }
      det = detbak;
      t1 = t1bak;
      t2 = t2bak;
      t3 = t3bak;
      t4 = t4bak;
      imap = imapbak;
      m_lastElement = imap;
      return imap;
    }
  }
 
  if (m_debug) {
    std::cout << m_className << "::FindElement5:\n";
    std::cout << "    No element matching point (" << x << ", " << y
              << ") found.\n";
  }
  return -1;
}

Referenced by Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys121::WeightingField(), and Garfield::ComponentAnsys121::WeightingPotential().

◆ FindElementCube()

int Garfield::ComponentFieldMap::FindElementCube	(	const double	x,
		const double	y,
		const double	z,
		double &	t1,
		double &	t2,
		double &	t3,
		TMatrixD *&	jac,
		std::vector< TMatrixD * > &	dN
	)

protected

Find the element for a point in a cube.

Definition at line 425 of file ComponentFieldMap.cc.

                                                                 {
  int imap = -1;
  if (m_lastElement >= 0) {
    const Element& element = elements[m_lastElement];
    const Node& n3 = nodes[element.emap[3]];
    if (x >= n3.x && y >= n3.y && z >= n3.z) {
      const Node& n0 = nodes[element.emap[0]];
      const Node& n2 = nodes[element.emap[2]];
      const Node& n7 = nodes[element.emap[7]];
      if (x < n0.x && y < n2.y && z < n7.z) {
        imap = m_lastElement;
      }
    }
  }
 
  // Default element loop
  if (imap == -1) {
    for (int i = 0; i < nElements; ++i) {
      const Element& element = elements[i];
      const Node& n3 = nodes[element.emap[3]];
      if (x < n3.x || y < n3.y || z < n3.z) continue;
      const Node& n0 = nodes[element.emap[0]];
      const Node& n2 = nodes[element.emap[2]];
      const Node& n7 = nodes[element.emap[7]];
      if (x < n0.x && y < n2.y && z < n7.z) {
        imap = i;
        break;
      }
    }
  }
 
  if (imap < 0) {
    if (m_debug) {
      std::cout << m_className << "::FindElementCube:\n";
      std::cout << "    Point (" << x << "," << y << "," << z
                << ") not in the mesh, it is background or PEC.\n";
      const Node& first0 = nodes[elements.front().emap[0]];
      const Node& first2 = nodes[elements.front().emap[2]];
      const Node& first3 = nodes[elements.front().emap[3]];
      const Node& first7 = nodes[elements.front().emap[7]];
      std::cout << "    First node (" << first3.x << "," << first3.y << ","
                << first3.z << ") in the mesh.\n";
      std::cout << "  dx= " << (first0.x - first3.x)
                << ", dy= " << (first2.y - first3.y)
                << ", dz= " << (first7.z - first3.z) << "\n";
      const Node& last0 = nodes[elements.back().emap[0]];
      const Node& last2 = nodes[elements.back().emap[2]];
      const Node& last3 = nodes[elements.back().emap[3]];
      const Node& last5 = nodes[elements.back().emap[5]];
      const Node& last7 = nodes[elements.back().emap[7]];
      std::cout << "    Last node (" << last5.x << "," << last5.y << ","
                << last5.z << ") in the mesh.\n";
      std::cout << "  dx= " << (last0.x - last3.x)
                << ", dy= " << (last2.y - last3.y)
                << ", dz= " << (last7.z - last3.z) << "\n";
    }
    return -1;
  }
  CoordinatesCube(x, y, z, t1, t2, t3, jac, dN, elements[imap]);
  if (m_debug) {
    PrintElement("FindElementCube", x, y, z, t1, t2, t3, 0., elements[imap], 8);
  }
  return imap;
}

◆ GetAspectRatio()

virtual void Garfield::ComponentFieldMap::GetAspectRatio	(	const unsigned int	i,
		double &	dmin,
		double &	dmax
	)

protectedpure virtual

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentComsol, Garfield::ComponentCST, and Garfield::ComponentElmer.

Referenced by GetElement().

◆ GetBoundingBox()

bool Garfield::ComponentFieldMap::GetBoundingBox	(	double &	xmin,
		double &	ymin,
		double &	zmin,
		double &	xmax,
		double &	ymax,
		double &	zmax
	)

overridevirtual

Get the bounding box coordinates.

Reimplemented from Garfield::ComponentBase.

Definition at line 2060 of file ComponentFieldMap.cc.

                                                     {
  if (!m_ready) return false;
 
  xmin = m_minBoundingBox[0];
  xmax = m_maxBoundingBox[0];
  ymin = m_minBoundingBox[1];
  ymax = m_maxBoundingBox[1];
  zmin = m_minBoundingBox[2];
  zmax = m_maxBoundingBox[2];
  return true;
}

◆ GetConductivity()

double Garfield::ComponentFieldMap::GetConductivity ( const unsigned int imat ) const

Return the conductivity of a field map material.

Definition at line 88 of file ComponentFieldMap.cc.

                                                                       {
  if (imat >= m_nMaterials) {
    std::cerr << m_className << "::GetConductivity: Index out of range.\n";
    return -1.;
  }
 
  return materials[imat].ohm;
}

◆ GetElement()

bool Garfield::ComponentFieldMap::GetElement	(	const unsigned int	i,
		double &	vol,
		double &	dmin,
		double &	dmax
	)

Return the volume and aspect ratio of a mesh element.

Definition at line 127 of file ComponentFieldMap.cc.

                                                               {
  if ((int)i >= nElements) {
    std::cerr << m_className << "::GetElement:\n";
    std::cerr << "    Element index (" << i << ") out of range.\n";
    return false;
  }
 
  vol = GetElementVolume(i);
  GetAspectRatio(i, dmin, dmax);
  return true;
}

◆ GetElementVolume()

virtual double Garfield::ComponentFieldMap::GetElementVolume ( const unsigned int i )

protectedpure virtual

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentComsol, Garfield::ComponentCST, and Garfield::ComponentElmer.

Referenced by GetElement().

◆ GetMedium() [1/2]

Medium * Garfield::ComponentBase::GetMedium	(	const double	x,
		const double	y,
		const double	z
	)

virtual

Get the medium at a given location (x, y, z).

Reimplemented from Garfield::ComponentBase.

Definition at line 26 of file ComponentBase.cc.

                                                 {
  if (!m_geometry) return nullptr;
  return m_geometry->GetMedium(x, y, z);
}

◆ GetMedium() [2/2]

Medium * Garfield::ComponentFieldMap::GetMedium ( const unsigned int i ) const

Return the Medium associated to a field map material.

Definition at line 117 of file ComponentFieldMap.cc.

                                                                  {
  if (imat >= m_nMaterials) {
    std::cerr << m_className << "::GetMedium:\n"
              << "    Material index " << imat << " is out of range.\n";
    return nullptr;
  }
 
  return materials[imat].medium;
}

◆ GetNumberOfElements()

int Garfield::ComponentFieldMap::GetNumberOfElements ( ) const

inline

Return the number of mesh elements.

Definition at line 62 of file ComponentFieldMap.hh.

62{ return nElements; }

◆ GetNumberOfMaterials()

unsigned int Garfield::ComponentFieldMap::GetNumberOfMaterials ( ) const

inline

Return the number of materials in the field map.

Definition at line 48 of file ComponentFieldMap.hh.

48{ return m_nMaterials; }

◆ GetNumberOfMedia()

unsigned int Garfield::ComponentFieldMap::GetNumberOfMedia ( ) const

inline

Definition at line 59 of file ComponentFieldMap.hh.

59{ return m_nMaterials; }

◆ GetPermittivity()

double Garfield::ComponentFieldMap::GetPermittivity ( const unsigned int imat ) const

Return the permittivity of a field map material.

Definition at line 79 of file ComponentFieldMap.cc.

                                                                       {
  if (imat >= m_nMaterials) {
    std::cerr << m_className << "::GetPermittivity: Index out of range.\n";
    return -1.;
  }
 
  return materials[imat].eps;
}

◆ GetVoltageRange()

bool Garfield::ComponentFieldMap::GetVoltageRange	(	double &	vmin,
		double &	vmax
	)

inlineoverridevirtual

Calculate the voltage range [V].

Implements Garfield::ComponentBase.

Definition at line 35 of file ComponentFieldMap.hh.

                                                            {
    vmin = m_mapvmin;
    vmax = m_mapvmax;
    return true;
  }

Referenced by Garfield::ViewFEMesh::Plot().

◆ IsInBoundingBox()

bool Garfield::ComponentFieldMap::IsInBoundingBox	(	const double	x,
		const double	y,
		const double	z
	)		const

inline

Definition at line 25 of file ComponentFieldMap.hh.

                                                                             {
    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];
  }

◆ MapCoordinates()

void Garfield::ComponentFieldMap::MapCoordinates	(	double &	xpos,
		double &	ypos,
		double &	zpos,
		bool &	xmirrored,
		bool &	ymirrored,
		bool &	zmirrored,
		double &	rcoordinate,
		double &	rotation
	)		const

protected

Move (xpos, ypos, zpos) to field map coordinates.

Definition at line 2074 of file ComponentFieldMap.cc.

                                                               {
  // Initial values
  rotation = 0;
 
  // If chamber is periodic, reduce to the cell volume.
  xmirrored = false;
  double auxr, auxphi;
  if (m_periodic[0]) {
    xpos = m_mapmin[0] + fmod(xpos - m_mapmin[0], m_mapmax[0] - m_mapmin[0]);
    if (xpos < m_mapmin[0]) xpos += m_mapmax[0] - m_mapmin[0];
  } else if (m_mirrorPeriodic[0]) {
    double xnew =
        m_mapmin[0] + fmod(xpos - m_mapmin[0], m_mapmax[0] - m_mapmin[0]);
    if (xnew < m_mapmin[0]) xnew += m_mapmax[0] - m_mapmin[0];
    int nx = int(floor(0.5 + (xnew - xpos) / (m_mapmax[0] - m_mapmin[0])));
    if (nx != 2 * (nx / 2)) {
      xnew = m_mapmin[0] + m_mapmax[0] - xnew;
      xmirrored = true;
    }
    xpos = xnew;
  }
  if (m_axiallyPeriodic[0] && (zpos != 0 || ypos != 0)) {
    auxr = sqrt(zpos * zpos + ypos * ypos);
    auxphi = atan2(zpos, ypos);
    rotation = (m_mapamax[0] - m_mapamin[0]) *
               floor(0.5 +
                     (auxphi - 0.5 * (m_mapamin[0] + m_mapamax[0])) /
                         (m_mapamax[0] - m_mapamin[0]));
    if (auxphi - rotation < m_mapamin[0])
      rotation = rotation - (m_mapamax[0] - m_mapamin[0]);
    if (auxphi - rotation > m_mapamax[0])
      rotation = rotation + (m_mapamax[0] - m_mapamin[0]);
    auxphi = auxphi - rotation;
    ypos = auxr * cos(auxphi);
    zpos = auxr * sin(auxphi);
  }
 
  ymirrored = false;
  if (m_periodic[1]) {
    ypos = m_mapmin[1] + fmod(ypos - m_mapmin[1], m_mapmax[1] - m_mapmin[1]);
    if (ypos < m_mapmin[1]) ypos += m_mapmax[1] - m_mapmin[1];
  } else if (m_mirrorPeriodic[1]) {
    double ynew =
        m_mapmin[1] + fmod(ypos - m_mapmin[1], m_mapmax[1] - m_mapmin[1]);
    if (ynew < m_mapmin[1]) ynew += m_mapmax[1] - m_mapmin[1];
    int ny = int(floor(0.5 + (ynew - ypos) / (m_mapmax[1] - m_mapmin[1])));
    if (ny != 2 * (ny / 2)) {
      ynew = m_mapmin[1] + m_mapmax[1] - ynew;
      ymirrored = true;
    }
    ypos = ynew;
  }
  if (m_axiallyPeriodic[1] && (xpos != 0 || zpos != 0)) {
    auxr = sqrt(xpos * xpos + zpos * zpos);
    auxphi = atan2(xpos, zpos);
    rotation = (m_mapamax[1] - m_mapamin[1]) *
               floor(0.5 +
                     (auxphi - 0.5 * (m_mapamin[1] + m_mapamax[1])) /
                         (m_mapamax[1] - m_mapamin[1]));
    if (auxphi - rotation < m_mapamin[1])
      rotation = rotation - (m_mapamax[1] - m_mapamin[1]);
    if (auxphi - rotation > m_mapamax[1])
      rotation = rotation + (m_mapamax[1] - m_mapamin[1]);
    auxphi = auxphi - rotation;
    zpos = auxr * cos(auxphi);
    xpos = auxr * sin(auxphi);
  }
 
  zmirrored = false;
  if (m_periodic[2]) {
    zpos = m_mapmin[2] + fmod(zpos - m_mapmin[2], m_mapmax[2] - m_mapmin[2]);
    if (zpos < m_mapmin[2]) zpos += m_mapmax[2] - m_mapmin[2];
  } else if (m_mirrorPeriodic[2]) {
    double znew =
        m_mapmin[2] + fmod(zpos - m_mapmin[2], m_mapmax[2] - m_mapmin[2]);
    if (znew < m_mapmin[2]) znew += m_mapmax[2] - m_mapmin[2];
    int nz = int(floor(0.5 + (znew - zpos) / (m_mapmax[2] - m_mapmin[2])));
    if (nz != 2 * (nz / 2)) {
      znew = m_mapmin[2] + m_mapmax[2] - znew;
      zmirrored = true;
    }
    zpos = znew;
  }
  if (m_axiallyPeriodic[2] && (ypos != 0 || xpos != 0)) {
    auxr = sqrt(ypos * ypos + xpos * xpos);
    auxphi = atan2(ypos, xpos);
    rotation = (m_mapamax[2] - m_mapamin[2]) *
               floor(0.5 +
                     (auxphi - 0.5 * (m_mapamin[2] + m_mapamax[2])) /
                         (m_mapamax[2] - m_mapamin[2]));
    if (auxphi - rotation < m_mapamin[2])
      rotation = rotation - (m_mapamax[2] - m_mapamin[2]);
    if (auxphi - rotation > m_mapamax[2])
      rotation = rotation + (m_mapamax[2] - m_mapamin[2]);
    auxphi = auxphi - rotation;
    xpos = auxr * cos(auxphi);
    ypos = auxr * sin(auxphi);
  }
 
  // If we have a rotationally symmetric field map, store coordinates.
  rcoordinate = 0;
  double zcoordinate = 0;
  if (m_rotationSymmetric[0]) {
    rcoordinate = sqrt(ypos * ypos + zpos * zpos);
    zcoordinate = xpos;
  } else if (m_rotationSymmetric[1]) {
    rcoordinate = sqrt(xpos * xpos + zpos * zpos);
    zcoordinate = ypos;
  } else if (m_rotationSymmetric[2]) {
    rcoordinate = sqrt(xpos * xpos + ypos * ypos);
    zcoordinate = zpos;
  }
 
  if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
      m_rotationSymmetric[2]) {
    xpos = rcoordinate;
    ypos = zcoordinate;
    zpos = 0;
  }
}

◆ NotDriftMedium()

void Garfield::ComponentFieldMap::NotDriftMedium ( const unsigned int imat )

Flag a field map materials as a non-drift medium.

Definition at line 65 of file ComponentFieldMap.cc.

                                                              {
  // Do not proceed if not properly initialised.
  if (!m_ready) PrintNotReady("NotDriftMedium");
 
  // Check value
  if (imat >= m_nMaterials) {
    std::cerr << m_className << "::NotDriftMedium: Index out of range.\n";
    return;
  }
 
  // Make drift medium
  materials[imat].driftmedium = false;
}

◆ PrintElement()

void Garfield::ComponentFieldMap::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 Element &	element,
		const unsigned int	n,
		const int	iw = `-1`
	)		const

protected

Definition at line 2396 of file ComponentFieldMap.cc.

                                                                               {
  std::cout << m_className << "::" << header << ":\n"
            << "    Global = (" << x << ", " << y << ", " << z << ")\n"
            << "    Local = (" << t1 << ", " << t2 << ", " << t3 << ", " << t4
            << ")\n";
  if (element.degenerate) std::cout << "    Element is degenerate.\n";
  std::cout << " Node             x            y            z            V\n";
  for (unsigned int ii = 0; ii < n; ++ii) {
    const Node& node = nodes[element.emap[ii]];
    const double v = iw < 0 ? node.v : node.w[iw];
    printf("      %-5d %12g %12g %12g %12g\n", element.emap[ii], node.x, node.y,
           node.z, v);
  }
}

◆ PrintMaterials()

void Garfield::ComponentFieldMap::PrintMaterials ( )

List all currently defined materials.

Definition at line 22 of file ComponentFieldMap.cc.

                                       {
  // Do not proceed if not properly initialised.
  if (!m_ready) PrintNotReady("PrintMaterials");
 
  if (materials.empty()) {
    std::cerr << m_className << "::PrintMaterials:\n"
              << "    No materials are currently defined.\n";
    return;
  }
 
  std::cout << m_className << "::PrintMaterials:\n"
            << "    Currently " << m_nMaterials << " materials are defined.\n"
            << "      Index Permittivity  Resistivity Notes\n";
  for (unsigned int i = 0; i < m_nMaterials; ++i) {
    printf("      %5d %12g %12g", i, materials[i].eps, materials[i].ohm);
    if (materials[i].medium) {
      std::string name = materials[i].medium->GetName();
      std::cout << " " << name;
      if (materials[i].medium->IsDriftable()) std::cout << ", drift medium";
      if (materials[i].medium->IsIonisable()) std::cout << ", ionisable";
    }
    if (materials[i].driftmedium) {
      std::cout << " (drift medium)\n";
    } else {
      std::cout << "\n";
    }
  }
}

Referenced by Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), and Garfield::ComponentAnsys123::Initialise().

◆ PrintNotReady()

void Garfield::ComponentFieldMap::PrintNotReady ( const std::string & header ) const

inlineprotected

Definition at line 199 of file ComponentFieldMap.hh.

                                                    {
    std::cerr << m_className << "::" << header << ":\n"
              << "    Field map not yet initialised.\n";
  }

Referenced by DriftMedium(), Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentComsol::ElectricField(), Garfield::ComponentElmer::ElectricField(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentComsol::GetMedium(), Garfield::ComponentElmer::GetMedium(), NotDriftMedium(), PrintMaterials(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), and Garfield::ComponentElmer::SetWeightingField().

◆ PrintRange()

void Garfield::ComponentFieldMap::PrintRange ( )

Show x, y, z, V and angular ranges.

Definition at line 2031 of file ComponentFieldMap.cc.

                                   {
  std::cout << m_className << "::PrintRange:\n";
  std::cout << "        Dimensions of the elementary block\n";
  printf("            %15g < x < %-15g cm,\n", m_mapmin[0], m_mapmax[0]);
  printf("            %15g < y < %-15g cm,\n", m_mapmin[1], m_mapmax[1]);
  printf("            %15g < z < %-15g cm,\n", m_mapmin[2], m_mapmax[2]);
  printf("            %15g < V < %-15g V.\n", m_mapvmin, m_mapvmax);
 
  std::cout << "        Periodicities\n";
  const std::array<std::string, 3> axes = {{"x", "y", "z"}};
  for (unsigned int i = 0; i < 3; ++i) {
    std::cout << "            " << axes[i] << ":";
    if (m_periodic[i]) {
      std::cout << " simple with length " << m_cells[i] << " cm";
    }
    if (m_mirrorPeriodic[i]) {
      std::cout << " mirror with length " << m_cells[i] << " cm";
    }
    if (m_axiallyPeriodic[i]) {
      std::cout << " axial " << int(0.5 + m_mapna[i]) << "-fold repetition";
    }
    if (m_rotationSymmetric[i]) std::cout << " rotational symmetry";
    if (!(m_periodic[i] || m_mirrorPeriodic[i] || m_axiallyPeriodic[i] ||
          m_rotationSymmetric[i]))
      std::cout << " none";
    std::cout << "\n";
  }
}

Referenced by Garfield::ComponentCST::Coordinate2Index(), SetRange(), and UpdatePeriodicityCommon().

◆ PrintWarning()

void Garfield::ComponentFieldMap::PrintWarning ( const std::string & header )

inlineprotected

Definition at line 193 of file ComponentFieldMap.hh.

                                             {
    if (!m_warning || m_nWarnings > 10) return;
    std::cerr << m_className << "::" << header << ":\n"
              << "    Warnings have been issued for this field map.\n";
    ++m_nWarnings;
  }

◆ ReadDouble()

double Garfield::ComponentFieldMap::ReadDouble	(	char *	token,
		double	def,
		bool &	error
	)

protected

Definition at line 2282 of file ComponentFieldMap.cc.

                                                                         {
  if (!token) {
    error = true;
    return def;
  }
  return atof(token);
}

Referenced by Garfield::ComponentElmer::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), and Garfield::ComponentCST::SetWeightingField().

◆ ReadInteger()

int Garfield::ComponentFieldMap::ReadInteger	(	char *	token,
		int	def,
		bool &	error
	)

protected

Definition at line 2273 of file ComponentFieldMap.cc.

                                                                    {
  if (!token) {
    error = true;
    return def;
  }
 
  return atoi(token);
}

Referenced by Garfield::ComponentElmer::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), and Garfield::ComponentCST::SetWeightingField().

◆ Reset()

void Garfield::ComponentFieldMap::Reset ( )

inlineoverrideprotectedvirtual

Reset the component.

Implements Garfield::ComponentBase.

Definition at line 157 of file ComponentFieldMap.hh.

157{};

◆ SetMedium()

void Garfield::ComponentFieldMap::SetMedium	(	const unsigned int	imat,
		Medium *	medium
	)

Associate a field map material with a Medium class.

Definition at line 97 of file ComponentFieldMap.cc.

                                                                    {
  if (imat >= m_nMaterials) {
    std::cerr << m_className << "::SetMedium:\n";
    std::cerr << "    Material index " << imat << " is out of range.\n";
    return;
  }
 
  if (!m) {
    std::cerr << m_className << "::SetMedium:    Null pointer.\n";
    return;
  }
 
  if (m_debug) {
    std::cout << m_className << "::SetMedium:\n    Associated material " << imat
              << " with medium " << m->GetName() << ".\n";
  }
 
  materials[imat].medium = m;
}

Referenced by main().

◆ SetRange()

void Garfield::ComponentFieldMap::SetRange ( )

virtual

Calculate x, y, z, V and angular ranges.

Reimplemented in Garfield::ComponentCST.

Definition at line 1934 of file ComponentFieldMap.cc.

                                 {
  // Initial values
  m_mapmin.fill(0.);
  m_mapmax.fill(0.);
  m_mapamin.fill(0.);
  m_mapamax.fill(0.);
  m_mapvmin = m_mapvmax = 0.;
  m_setang.fill(false);
 
  // Make sure the required data is available.
  if (!m_ready || nNodes < 1) {
    std::cerr << m_className << "::SetRange:\n";
    std::cerr << "    Field map not yet set.\n";
    return;
  }
  if (nNodes < 1) {
    std::cerr << m_className << "::SetRange:\n";
    std::cerr << "    Number of nodes < 1.\n";
    return;
  }
 
  // Loop over the nodes.
  m_mapmin[0] = m_mapmax[0] = nodes[0].x;
  m_mapmin[1] = m_mapmax[1] = nodes[0].y;
  m_mapmin[2] = m_mapmax[2] = nodes[0].z;
  m_mapvmin = m_mapvmax = nodes[0].v;
 
  for (const auto& node : nodes) {
    const std::array<double, 3> pos = {{node.x, node.y, node.z}};
    for (unsigned int i = 0; i < 3; ++i) {
      m_mapmin[i] = std::min(m_mapmin[i], pos[i]);
      m_mapmax[i] = std::max(m_mapmax[i], pos[i]);
    }
    m_mapvmin = std::min(m_mapvmin, node.v);
    m_mapvmax = std::max(m_mapvmax, node.v);
 
    if (node.y != 0 || node.z != 0) {
      const double ang = atan2(node.z, node.y);
      if (m_setang[0]) {
        m_mapamin[0] = std::min(m_mapamin[0], ang);
        m_mapamax[0] = std::max(m_mapamax[0], ang);
      } else {
        m_mapamin[0] = m_mapamax[0] = ang;
        m_setang[0] = true;
      }
    }
 
    if (node.z != 0 || node.x != 0) {
      const double ang = atan2(node.x, node.z);
      if (m_setang[1]) {
        m_mapamin[1] = std::min(m_mapamin[1], ang);
        m_mapamax[1] = std::max(m_mapamax[1], ang);
      } else {
        m_mapamin[1] = m_mapamax[1] = ang;
        m_setang[1] = true;
      }
    }
 
    if (node.x != 0 || node.y != 0) {
      const double ang = atan2(node.y, node.x);
      if (m_setang[2]) {
        m_mapamin[2] = std::min(m_mapamin[2], ang);
        m_mapamax[2] = std::max(m_mapamax[2], ang);
      } else {
        m_mapamin[2] = m_mapamax[2] = ang;
        m_setang[2] = true;
      }
    }
  }
 
  // Fix the angular ranges.
  for (unsigned int i = 0; i < 3; ++i) {
    if (m_mapamax[i] - m_mapamin[i] > Pi) {
      const double aux = m_mapamin[i];
      m_mapamin[i] = m_mapamax[i];
      m_mapamax[i] = aux + TwoPi;
    }
  }
 
  // Set provisional cell dimensions.
  m_minBoundingBox[0] = m_mapmin[0];
  m_maxBoundingBox[0] = m_mapmax[0];
  m_minBoundingBox[1] = m_mapmin[1];
  m_maxBoundingBox[1] = m_mapmax[1];
  if (m_is3d) {
    m_minBoundingBox[2] = m_mapmin[2];
    m_maxBoundingBox[2] = m_mapmax[2];
  } else {
    m_mapmin[2] = m_minBoundingBox[2];
    m_mapmax[2] = m_maxBoundingBox[2];
  }
  hasBoundingBox = true;
 
  // Display the range if requested.
  if (m_debug) PrintRange();
}

Referenced by Garfield::ComponentElmer::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), and Garfield::ComponentComsol::Initialise().

◆ UnmapFields()

void Garfield::ComponentFieldMap::UnmapFields	(	double &	ex,
		double &	ey,
		double &	ez,
		double &	xpos,
		double &	ypos,
		double &	zpos,
		bool &	xmirrored,
		bool &	ymirrored,
		bool &	zmirrored,
		double &	rcoordinate,
		double &	rotation
	)		const

protected

Move (ex, ey, ez) to global coordinates.

Definition at line 2198 of file ComponentFieldMap.cc.

                                                            {
  // Apply mirror imaging.
  if (xmirrored) ex = -ex;
  if (ymirrored) ey = -ey;
  if (zmirrored) ez = -ez;
 
  // Rotate the field.
  double er, theta;
  if (m_axiallyPeriodic[0]) {
    er = sqrt(ey * ey + ez * ez);
    theta = atan2(ez, ey);
    theta += rotation;
    ey = er * cos(theta);
    ez = er * sin(theta);
  }
  if (m_axiallyPeriodic[1]) {
    er = sqrt(ez * ez + ex * ex);
    theta = atan2(ex, ez);
    theta += rotation;
    ez = er * cos(theta);
    ex = er * sin(theta);
  }
  if (m_axiallyPeriodic[2]) {
    er = sqrt(ex * ex + ey * ey);
    theta = atan2(ey, ex);
    theta += rotation;
    ex = er * cos(theta);
    ey = er * sin(theta);
  }
 
  // Take care of symmetry.
  double eaxis;
  er = ex;
  eaxis = ey;
 
  // Rotational symmetry
  if (m_rotationSymmetric[0]) {
    if (rcoordinate <= 0) {
      ex = eaxis;
      ey = 0;
      ez = 0;
    } else {
      ex = eaxis;
      ey = er * ypos / rcoordinate;
      ez = er * zpos / rcoordinate;
    }
  }
  if (m_rotationSymmetric[1]) {
    if (rcoordinate <= 0) {
      ex = 0;
      ey = eaxis;
      ez = 0;
    } else {
      ex = er * xpos / rcoordinate;
      ey = eaxis;
      ez = er * zpos / rcoordinate;
    }
  }
  if (m_rotationSymmetric[2]) {
    if (rcoordinate <= 0) {
      ex = 0;
      ey = 0;
      ez = eaxis;
    } else {
      ex = er * xpos / rcoordinate;
      ey = er * ypos / rcoordinate;
      ez = eaxis;
    }
  }
}

Referenced by Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentComsol::ElectricField(), Garfield::ComponentElmer::ElectricField(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentComsol::WeightingField(), and Garfield::ComponentElmer::WeightingField().

◆ UpdatePeriodicity2d()

void Garfield::ComponentFieldMap::UpdatePeriodicity2d ( )

protected

Definition at line 1905 of file ComponentFieldMap.cc.

                                            {
  // Check the required data is available.
  if (!m_ready) {
    std::cerr << m_className << "::UpdatePeriodicity2d:\n";
    std::cerr << "    No valid field map available.\n";
    return;
  }
 
  // No z-periodicity in 2d
  if (m_periodic[2] || m_mirrorPeriodic[2]) {
    std::cerr << m_className << "::UpdatePeriodicity2d:\n";
    std::cerr << "    Simple or mirror periodicity along z\n";
    std::cerr << "    requested for a 2d map; reset.\n";
    m_periodic[2] = false;
    m_mirrorPeriodic[2] = false;
    m_warning = true;
  }
 
  // Only z-axial periodicity in 2d maps
  if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1]) {
    std::cerr << m_className << "::UpdatePeriodicity2d:\n";
    std::cerr << "    Axial symmetry has been requested \n";
    std::cerr << "    around x or y for a 2D map; reset.\n";
    m_axiallyPeriodic[0] = false;
    m_axiallyPeriodic[1] = false;
    m_warning = true;
  }
}

Referenced by Garfield::ComponentAnsys121::UpdatePeriodicity(), and Garfield::ComponentCST::UpdatePeriodicity().

◆ UpdatePeriodicityCommon()

void Garfield::ComponentFieldMap::UpdatePeriodicityCommon ( )

protected

Definition at line 1753 of file ComponentFieldMap.cc.

                                                {
  // Check the required data is available.
  if (!m_ready) {
    std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
    std::cerr << "    No valid field map available.\n";
    return;
  }
 
  for (unsigned int i = 0; i < 3; ++i) {
    // No regular and mirror periodicity at the same time.
    if (m_periodic[i] && m_mirrorPeriodic[i]) {
      std::cerr << m_className << "::UpdatePeriodicityCommon:\n"
                << "    Both simple and mirror periodicity requested. Reset.\n";
      m_periodic[i] = false;
      m_mirrorPeriodic[i] = false;
      m_warning = true;
    }
    // In case of axial periodicity,
    // the range must be an integral part of two pi.
    if (m_axiallyPeriodic[i]) {
      if (m_mapamin[i] >= m_mapamax[i]) {
        m_mapna[i] = 0;
      } else {
        m_mapna[i] = TwoPi / (m_mapamax[i] - m_mapamin[i]);
      }
      if (fabs(m_mapna[i] - int(0.5 + m_mapna[i])) > 0.001 ||
          m_mapna[i] < 1.5) {
        std::cerr
            << m_className << "::UpdatePeriodicityCommon:\n"
            << "    Axial symmetry has been requested but the map\n"
            << "    does not cover an integral fraction of 2 pi. Reset.\n";
        m_axiallyPeriodic[i] = false;
        m_warning = true;
      }
    }
  }
 
  // Not more than 1 rotational symmetry
  if ((m_rotationSymmetric[0] && m_rotationSymmetric[1]) ||
      (m_rotationSymmetric[0] && m_rotationSymmetric[2]) ||
      (m_rotationSymmetric[1] && m_rotationSymmetric[2])) {
    std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
    std::cerr << "    Only 1 rotational symmetry allowed; reset.\n";
    m_rotationSymmetric.fill(false);
    m_warning = true;
  }
 
  // No rotational symmetry as well as axial periodicity
  if ((m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
       m_rotationSymmetric[2]) &&
      (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2])) {
    std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
    std::cerr << "    Not allowed to combine rotational symmetry\n";
    std::cerr << "    and axial periodicity; reset.\n";
    m_axiallyPeriodic.fill(false);
    m_rotationSymmetric.fill(false);
    m_warning = true;
  }
 
  // In case of rotational symmetry, the x-range should not straddle 0.
  if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
      m_rotationSymmetric[2]) {
    if (m_mapmin[0] * m_mapmax[0] < 0) {
      std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
      std::cerr << "    Rotational symmetry requested, \n";
      std::cerr << "    but x-range straddles 0; reset.\n";
      m_rotationSymmetric.fill(false);
      m_warning = true;
    }
  }
 
  // Recompute the cell ranges.
  for (unsigned int i = 0; i < 3; ++i) {
    m_minBoundingBox[i] = m_mapmin[i];
    m_maxBoundingBox[i] = m_mapmax[i];
    m_cells[i] = fabs(m_mapmax[i] - m_mapmin[i]);
  }
  if (m_rotationSymmetric[0]) {
    m_minBoundingBox[0] = m_mapmin[1];
    m_maxBoundingBox[0] = m_mapmax[1];
    m_minBoundingBox[1] = -std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_maxBoundingBox[1] = +std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_minBoundingBox[2] = -std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_maxBoundingBox[2] = +std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
  } else if (m_rotationSymmetric[1]) {
    m_minBoundingBox[0] = -std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_maxBoundingBox[0] = +std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_minBoundingBox[1] = m_mapmin[1];
    m_maxBoundingBox[1] = m_mapmax[1];
    m_minBoundingBox[2] = -std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_maxBoundingBox[2] = +std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
  } else if (m_rotationSymmetric[2]) {
    m_minBoundingBox[0] = -std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_maxBoundingBox[0] = +std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_minBoundingBox[1] = -std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_maxBoundingBox[1] = +std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0]));
    m_minBoundingBox[2] = m_mapmin[1];
    m_maxBoundingBox[2] = m_mapmax[1];
  }
 
  if (m_axiallyPeriodic[0]) {
    m_minBoundingBox[1] =
        -std::max(std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
    m_maxBoundingBox[1] =
        +std::max(std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
    m_minBoundingBox[2] =
        -std::max(std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
    m_maxBoundingBox[2] =
        +std::max(std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
  } else if (m_axiallyPeriodic[1]) {
    m_minBoundingBox[0] =
        -std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
    m_maxBoundingBox[0] =
        +std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
    m_minBoundingBox[2] =
        -std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
    m_maxBoundingBox[2] =
        +std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[2]), fabs(m_mapmax[2])));
  } else if (m_axiallyPeriodic[2]) {
    m_minBoundingBox[0] =
        -std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])));
    m_maxBoundingBox[0] =
        +std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])));
    m_minBoundingBox[1] =
        -std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])));
    m_maxBoundingBox[1] =
        +std::max(std::max(fabs(m_mapmin[0]), fabs(m_mapmax[0])),
                  std::max(fabs(m_mapmin[1]), fabs(m_mapmax[1])));
  }
 
  for (unsigned int i = 0; i < 3; ++i) {
    if (m_periodic[i] || m_mirrorPeriodic[i]) {
      m_minBoundingBox[i] = -INFINITY;
      m_maxBoundingBox[i] = +INFINITY;
    }
  }
 
  // Display the range if requested.
  if (m_debug) PrintRange();
}

Referenced by Garfield::ComponentAnsys121::UpdatePeriodicity(), Garfield::ComponentAnsys123::UpdatePeriodicity(), Garfield::ComponentComsol::UpdatePeriodicity(), Garfield::ComponentCST::UpdatePeriodicity(), and Garfield::ComponentElmer::UpdatePeriodicity().

Friends And Related Function Documentation

◆ ViewFEMesh

friend class ViewFEMesh

friend

Definition at line 84 of file ComponentFieldMap.hh.

Member Data Documentation

◆ elements

std::vector<Element> Garfield::ComponentFieldMap::elements

protected

Definition at line 100 of file ComponentFieldMap.hh.

◆ hasBoundingBox

bool Garfield::ComponentFieldMap::hasBoundingBox = false

protected

Definition at line 132 of file ComponentFieldMap.hh.

Referenced by SetRange(), and Garfield::ComponentCST::SetRange().

◆ m_cells

std::array<double, 3> Garfield::ComponentFieldMap::m_cells

protected

Definition at line 142 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ m_deleteBackground

bool Garfield::ComponentFieldMap::m_deleteBackground = true

protected

Definition at line 150 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::ComponentCST(), EnableDeleteBackgroundElements(), Garfield::ComponentAnsys121::Initialise(), and Garfield::ComponentAnsys123::Initialise().

◆ m_is3d

bool Garfield::ComponentFieldMap::m_is3d = true

protected

Definition at line 87 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), SetRange(), and Garfield::ComponentCST::SetRange().

◆ m_mapamax

std::array<double, 3> Garfield::ComponentFieldMap::m_mapamax

protected

Definition at line 140 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ m_mapamin

std::array<double, 3> Garfield::ComponentFieldMap::m_mapamin

protected

Definition at line 139 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ m_mapmax

std::array<double, 3> Garfield::ComponentFieldMap::m_mapmax

protected

Definition at line 138 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), SetRange(), Garfield::ComponentCST::SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_mapmin

std::array<double, 3> Garfield::ComponentFieldMap::m_mapmin

protected

Definition at line 137 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), SetRange(), Garfield::ComponentCST::SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_mapna

std::array<double, 3> Garfield::ComponentFieldMap::m_mapna

protected

Definition at line 141 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ m_mapvmax

double Garfield::ComponentFieldMap::m_mapvmax

protected

Definition at line 144 of file ComponentFieldMap.hh.

Referenced by GetVoltageRange(), PrintRange(), SetRange(), and Garfield::ComponentCST::SetRange().

◆ m_mapvmin

double Garfield::ComponentFieldMap::m_mapvmin

protected

Definition at line 144 of file ComponentFieldMap.hh.

Referenced by GetVoltageRange(), PrintRange(), SetRange(), and Garfield::ComponentCST::SetRange().

◆ m_maxBoundingBox

std::array<double, 3> Garfield::ComponentFieldMap::m_maxBoundingBox

protected

Definition at line 134 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::ComponentCST(), Garfield::ComponentAnsys121::ElectricField(), GetBoundingBox(), Garfield::ComponentAnsys121::GetMedium(), IsInBoundingBox(), SetRange(), Garfield::ComponentCST::SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentCST::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_minBoundingBox

std::array<double, 3> Garfield::ComponentFieldMap::m_minBoundingBox

protected

Definition at line 133 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::ComponentCST(), Garfield::ComponentAnsys121::ElectricField(), GetBoundingBox(), Garfield::ComponentAnsys121::GetMedium(), IsInBoundingBox(), SetRange(), Garfield::ComponentCST::SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentCST::SetRangeZ(), and UpdatePeriodicityCommon().

◆ m_nMaterials

unsigned int Garfield::ComponentFieldMap::m_nMaterials = 0

protected

◆ m_nWarnings

unsigned int Garfield::ComponentFieldMap::m_nWarnings = 0

protected

Definition at line 154 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentElmer::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), and PrintWarning().

◆ m_setang

std::array<bool, 3> Garfield::ComponentFieldMap::m_setang

protected

Definition at line 146 of file ComponentFieldMap.hh.

Referenced by SetRange().

◆ m_warning

bool Garfield::ComponentFieldMap::m_warning = false

protected

Definition at line 153 of file ComponentFieldMap.hh.

◆ materials

std::vector<Material> Garfield::ComponentFieldMap::materials

protected

◆ nElements

int Garfield::ComponentFieldMap::nElements = -1

protected

Definition at line 90 of file ComponentFieldMap.hh.

Referenced by FindElement13(), FindElement5(), FindElementCube(), Garfield::ComponentCST::GetAspectRatio(), GetElement(), Garfield::ComponentCST::GetElementVolume(), GetNumberOfElements(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), and Garfield::ComponentComsol::Initialise().

◆ nNodes

int Garfield::ComponentFieldMap::nNodes = -1

protected

Definition at line 103 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentElmer::Initialise(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentComsol::Initialise(), SetRange(), Garfield::ComponentComsol::SetWeightingField(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), and Garfield::ComponentCST::SetWeightingField().

◆ nodes

std::vector<Node> Garfield::ComponentFieldMap::nodes

protected

Definition at line 112 of file ComponentFieldMap.hh.

◆ nWeightingFields

int Garfield::ComponentFieldMap::nWeightingFields = 0

protected

Definition at line 127 of file ComponentFieldMap.hh.

◆ wfields

std::vector<std::string> Garfield::ComponentFieldMap::wfields

protected

Definition at line 128 of file ComponentFieldMap.hh.

◆ wfieldsOk

std::vector<bool> Garfield::ComponentFieldMap::wfieldsOk

protected

Definition at line 129 of file ComponentFieldMap.hh.

The documentation for this class was generated from the following files:

Classes

Public Member Functions

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

◆ ComponentFieldMap()

◆ ~ComponentFieldMap()

Member Function Documentation

◆ DisableCheckMapIndices()

◆ DriftMedium()

◆ EnableCheckMapIndices()

◆ EnableDeleteBackgroundElements()

◆ EnableTetrahedralTreeForElementSearch()

◆ FindElement13()

◆ FindElement5()

◆ FindElementCube()

◆ GetAspectRatio()

◆ GetBoundingBox()

◆ GetConductivity()

◆ GetElement()

◆ GetElementVolume()

◆ GetMedium() [1/2]

◆ GetMedium() [2/2]

◆ GetNumberOfElements()

◆ GetNumberOfMaterials()

◆ GetNumberOfMedia()

◆ GetPermittivity()

◆ GetVoltageRange()

◆ IsInBoundingBox()

◆ MapCoordinates()

◆ NotDriftMedium()

◆ PrintElement()

◆ PrintMaterials()

◆ PrintNotReady()

◆ PrintRange()

◆ PrintWarning()

◆ ReadDouble()

◆ ReadInteger()

◆ Reset()

◆ SetMedium()

◆ SetRange()

◆ UnmapFields()

◆ UpdatePeriodicity2d()

◆ UpdatePeriodicityCommon()

Friends And Related Function Documentation

◆ ViewFEMesh

Member Data Documentation

◆ elements

◆ hasBoundingBox

◆ m_cells

◆ m_deleteBackground

◆ m_is3d

◆ m_mapamax

◆ m_mapamin

◆ m_mapmax

◆ m_mapmin

◆ m_mapna

◆ m_mapvmax

◆ m_mapvmin

◆ m_maxBoundingBox

◆ m_minBoundingBox

◆ m_nMaterials

◆ m_nWarnings

◆ m_setang

◆ m_warning

◆ materials

◆ nElements

◆ nNodes

◆ nodes

◆ nWeightingFields

◆ wfields

◆ wfieldsOk