Component for importing and interpolating Comsol field maps. More...

#include <ComponentComsol.hh>

Inheritance diagram for Garfield::ComponentComsol:

Classes
struct	nodeCmp

Public Member Functions
	ComponentComsol ()

	ComponentComsol (std::string mesh, std::string mplist, std::string field)

	~ComponentComsol ()

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

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

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

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

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

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

bool	Initialise (std::string header="mesh.mphtxt", std::string mplist="dielectrics.dat", std::string field="field.txt")

bool	SetWeightingField (std::string file, std::string label)

Public Member Functions inherited from Garfield::ComponentFieldMap
	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.

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

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

virtual bool	GetVoltageRange (double &vmin, double &vmax)
	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.

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

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.

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

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

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

void	EnableCheckMapIndices ()

void	DisableCheckMapIndices ()

void	EnableDeleteBackgroundElements ()

void	DisableDeleteBackgroundElements ()

void	EnableTetrahedralTreeForElementSearch (const bool on=true)

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

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	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.

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)

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	UpdatePeriodicity ()
	Verify periodicities.

double	GetElementVolume (const unsigned int i)

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

Protected Member Functions inherited from Garfield::ComponentFieldMap
void	Reset ()
	Geometry checks.

virtual void	UpdatePeriodicity ()=0
	Verify periodicities.

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 unsigned int i, const unsigned int n, const int iw=-1) const

virtual void	Reset ()=0
	Geometry checks.

virtual void	UpdatePeriodicity ()=0
	Verify periodicities.

Additional Inherited Members
Protected Attributes inherited from Garfield::ComponentFieldMap
bool	m_is3d

int	nElements

std::vector< Element >	elements

int	nNodes

std::vector< Node >	nodes

unsigned int	m_nMaterials

std::vector< Material >	materials

int	nWeightingFields

std::vector< std::string >	wfields

std::vector< bool >	wfieldsOk

bool	hasBoundingBox

double	xMinBoundingBox

double	yMinBoundingBox

double	zMinBoundingBox

double	xMaxBoundingBox

double	yMaxBoundingBox

double	zMaxBoundingBox

double	mapxmin

double	mapymin

double	mapzmin

double	mapxmax

double	mapymax

double	mapzmax

double	mapxamin

double	mapyamin

double	mapzamin

double	mapxamax

double	mapyamax

double	mapzamax

double	mapvmin

double	mapvmax

bool	setangx

bool	setangy

bool	setangz

double	mapsx

double	mapsy

double	mapsz

double	cellsx

double	cellsy

double	cellsz

double	mapnxa

double	mapnya

double	mapnza

bool	m_deleteBackground

bool	m_warning

unsigned int	m_nWarnings

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

GeometryBase *	m_geometry
	Pointer to the geometry.

bool	m_ready
	Ready for use?

bool	m_activeTraps
	Does the component have traps?

bool	m_hasVelocityMap
	Does the component have velocity maps?

bool	m_xPeriodic
	Simple periodicity in x.

bool	m_yPeriodic
	Simple periodicity in y.

bool	m_zPeriodic
	Simple periodicity in z.

bool	m_xMirrorPeriodic
	Mirror periodicity in x.

bool	m_yMirrorPeriodic
	Mirror periodicity in y.

bool	m_zMirrorPeriodic
	Mirror periodicity in z.

bool	m_xAxiallyPeriodic
	Axial periodicity in x.

bool	m_yAxiallyPeriodic
	Axial periodicity in y.

bool	m_zAxiallyPeriodic
	Axial periodicity in z.

bool	m_xRotationSymmetry
	Rotation symmetry around x-axis.

bool	m_yRotationSymmetry
	Rotation symmetry around y-axis.

bool	m_zRotationSymmetry
	Rotation symmetry around z-axis.

double	m_bx0

double	m_by0

double	m_bz0

bool	m_debug
	Switch on/off debugging messages.

Detailed Description

Component for importing and interpolating Comsol field maps.

Definition at line 9 of file ComponentComsol.hh.

Constructor & Destructor Documentation

◆ ComponentComsol() [1/2]

Garfield::ComponentComsol::ComponentComsol ( )

Definition at line 14 of file ComponentComsol.cc.

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

◆ ComponentComsol() [2/2]

Garfield::ComponentComsol::ComponentComsol	(	std::string	mesh,
		std::string	mplist,
		std::string	field
	)

Definition at line 19 of file ComponentComsol.cc.

    : ComponentFieldMap() {
 
  m_className = "ComponentComsol";
  Initialise(mesh, mplist, field);
}

◆ ~ComponentComsol()

Garfield::ComponentComsol::~ComponentComsol ( )

inline

Definition at line 16 of file ComponentComsol.hh.

16{}

Member Function Documentation

◆ ElectricField() [1/2]

void Garfield::ComponentComsol::ElectricField	(	const double	x,
		const double	y,
		const double	z,
		double &	ex,
		double &	ey,
		double &	ez,
		double &	v,
		Medium *&	m,
		int &	status
	)

virtual

Implements Garfield::ComponentFieldMap.

Definition at line 344 of file ComponentComsol.cc.

                                                 {
 
  // Copy the coordinates
  double x = xin, y = yin, z = zin;
 
  // Map the coordinates onto field map coordinates
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  // Initial values
  ex = ey = ez = volt = 0.;
  status = 0;
  m = NULL;
 
  // Do not proceed if not properly initialised.
  if (!m_ready) {
    status = -10;
    PrintNotReady("ElectricField");
    return;
  }
 
  if (m_warning) PrintWarning("ElectricField");
 
  // Find the element that contains this point
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
  if (imap < 0) {
    if (m_debug) {
      std::cout << m_className << "::ElectricField:\n";
      std::cout << "    Point (" << x << ", " << y << ", " << z
                << " not in the mesh.\n";
    }
    status = -6;
    return;
  }
 
  if (m_debug) {
    PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, imap, 10);
  }
 
  const Element& element = elements[imap];
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const Node& n4 = nodes[element.emap[4]];
  const Node& n5 = nodes[element.emap[5]];
  const Node& n6 = nodes[element.emap[6]];
  const Node& n7 = nodes[element.emap[7]];
  const Node& n8 = nodes[element.emap[8]];
  const Node& n9 = nodes[element.emap[9]];
  // Tetrahedral field
  volt = n0.v * t1 * (2 * t1 - 1) + n1.v * t2 * (2 * t2 - 1) +
         n2.v * t3 * (2 * t3 - 1) + n3.v * t4 * (2 * t4 - 1) +
         4 * n4.v * t1 * t2 + 4 * n5.v * t1 * t3 + 4 * n6.v * t1 * t4 +
         4 * n7.v * t2 * t3 + 4 * n8.v * t2 * t4 + 4 * n9.v * t3 * t4;
  ex = -(n0.v * (4 * t1 - 1) * jac[0][1] + n1.v * (4 * t2 - 1) * jac[1][1] +
         n2.v * (4 * t3 - 1) * jac[2][1] + n3.v * (4 * t4 - 1) * jac[3][1] +
         n4.v * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
         n5.v * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
         n6.v * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
         n7.v * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
         n8.v * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
         n9.v * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
       det;
  ey = -(n0.v * (4 * t1 - 1) * jac[0][2] + n1.v * (4 * t2 - 1) * jac[1][2] +
         n2.v * (4 * t3 - 1) * jac[2][2] + n3.v * (4 * t4 - 1) * jac[3][2] +
         n4.v * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
         n5.v * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
         n6.v * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
         n7.v * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
         n8.v * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
         n9.v * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
       det;
  ez = -(n0.v * (4 * t1 - 1) * jac[0][3] + n1.v * (4 * t2 - 1) * jac[1][3] +
         n2.v * (4 * t3 - 1) * jac[2][3] + n3.v * (4 * t4 - 1) * jac[3][3] +
         n4.v * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
         n5.v * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
         n6.v * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
         n7.v * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
         n8.v * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
         n9.v * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
       det;
 
  // Transform field to global coordinates
  UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
  //  std::cout << "ef @(" << xin << ", " << yin << ", " << zin << ") = " <<
  // volt << "\n";
 
  // Drift medium?
  if (m_debug) {
    std::cout << m_className << "::ElectricField:\n";
    std::cout << "    Material " << element.matmap << ", drift flag "
              << materials[element.matmap].driftmedium << "\n";
  }
  m = materials[element.matmap].medium;
  status = -5;
  if (materials[element.matmap].driftmedium) {
    if (m && m->IsDriftable()) status = 0;
  }
}

◆ ElectricField() [2/2]

void Garfield::ComponentComsol::ElectricField	(	const double	x,
		const double	y,
		const double	z,
		double &	ex,
		double &	ey,
		double &	ez,
		Medium *&	m,
		int &	status
	)

virtual

Calculate the drift field at given point.

Parameters

x,y,z	coordinates [cm].
ex,ey,ez	components of the electric field [V/cm].
m	pointer to the medium at this location.
status	status flag

Status flags:

        0: Inside an active medium
      > 0: Inside a wire of type X
-4 ... -1: On the side of a plane where no wires are
       -5: Inside the mesh but not in an active medium
       -6: Outside the mesh
      -10: Unknown potential type (should not occur)
    other: Other cases (should not occur)

Implements Garfield::ComponentFieldMap.

Definition at line 336 of file ComponentComsol.cc.

                                                                         {
 
  double v = 0.;
  ElectricField(x, y, z, ex, ey, ez, v, m, status);
}

Referenced by ElectricField().

◆ GetAspectRatio()

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

protectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 677 of file ComponentComsol.cc.

                                                   {
 
  if (i >= elements.size()) {
    dmin = dmax = 0.;
    return;
  }
 
  const Element& element = elements[i];
  const int np = 4;
  // Loop over all pairs of vertices.
  for (int j = 0; j < np - 1; ++j) {
    const Node& nj = nodes[element.emap[j]];
    for (int k = j + 1; k < np; ++k) {
      const Node& nk = nodes[element.emap[k]];
      // Compute distance.
      const double dx = nj.x - nk.x;
      const double dy = nj.y - nk.y;
      const double dz = nj.z - nk.z;
      const double dist = sqrt(dx * dx + dy * dy + dz * dz);
      if (k == 1) {
        dmin = dmax = dist;
      } else {
        if (dist < dmin) dmin = dist;
        if (dist > dmax) dmax = dist;
      }
    }
  }
}

◆ GetElementVolume()

double Garfield::ComponentComsol::GetElementVolume ( const unsigned int i )

protectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 655 of file ComponentComsol.cc.

                                                             {
 
  if (i >= elements.size()) return 0.;
  const Element& element = elements[i];
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
 
  // Uses formula V = |a (dot) b x c|/6
  // with a => "3", b => "1", c => "2" and origin "0"
  const double vol =
      fabs((n3.x - n0.x) *
               ((n1.y - n0.y) * (n2.z - n0.z) - (n2.y - n0.y) * (n1.z - n0.z)) +
           (n3.y - n0.y) *
               ((n1.z - n0.z) * (n2.x - n0.x) - (n2.z - n0.z) * (n1.x - n0.x)) +
           (n3.z - n0.z) * ((n1.x - n0.x) * (n2.y - n0.y) -
                            (n3.x - n0.x) * (n1.y - n0.y))) /
      6.;
  return vol;
}

◆ GetMedium()

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

virtual

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

Implements Garfield::ComponentFieldMap.

Definition at line 609 of file ComponentComsol.cc.

                                                     {
 
  // Copy the coordinates
  double x = xin, y = yin, z = zin;
 
  // Map the coordinates onto field map coordinates
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  // Do not proceed if not properly initialised.
  if (!m_ready) {
    PrintNotReady("GetMedium");
    return nullptr;
  }
  if (m_warning) PrintWarning("GetMedium");
 
  // Find the element that contains this point
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
  if (imap < 0) {
    if (m_debug) {
      std::cout << m_className << "::GetMedium:\n";
      std::cout << "    Point (" << x << ", " << y << ", " << z
                << ") not in the mesh.\n";
    }
    return nullptr;
  }
  const Element& element = elements[imap];
  if (element.matmap >= m_nMaterials) {
    if (m_debug) {
      std::cerr << m_className << "::GetMedium:\n";
      std::cerr << "    Point (" << x << ", " << y
                << ") has out of range material number " << imap << ".\n";
    }
    return nullptr;
  }
 
  if (m_debug) {
    PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, imap, 10);
  }
 
  return materials[element.matmap].medium;
}

◆ Initialise()

bool Garfield::ComponentComsol::Initialise	(	std::string	header = `"mesh.mphtxt"`,
		std::string	mplist = `"dielectrics.dat"`,
		std::string	field = `"field.txt"`
	)

Definition at line 38 of file ComponentComsol.cc.

                                                  {
 
  m_ready = false;
  m_warning = false;
  m_nWarnings = 0;
 
  double unit = 100.0;  // m
 
  std::string line;
 
  // Open the materials file.
  materials.clear();
  std::ifstream fmplist;
  fmplist.open(mplist.c_str(), std::ios::in);
  if (fmplist.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open result file " << mplist
              << " for reading.\n";
    return false;
  }
  fmplist >> m_nMaterials;
  for (unsigned int i = 0; i < m_nMaterials; ++i) {
    Material newMaterial;
    newMaterial.driftmedium = true;
    newMaterial.medium = nullptr;
    newMaterial.ohm = -1;
    fmplist >> newMaterial.eps;
    materials.push_back(newMaterial);
  }
  {
    // add default material
    Material newMaterial;
    newMaterial.driftmedium = false;
    newMaterial.medium = nullptr;
    newMaterial.eps = newMaterial.ohm = -1;
    materials.push_back(newMaterial);
    m_nMaterials++;
  }
  std::map<int, int> domain2material;
  int d2msize;
  fmplist >> d2msize;
  for (int i = 0; i < d2msize; ++i) {
    int domain;
    fmplist >> domain;
    fmplist >> domain2material[domain];
  }
  fmplist.close();
 
  nodes.clear();
  std::ifstream fmesh;
  fmesh.open(mesh.c_str(), std::ios::in);
  if (fmesh.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open nodes file " << mesh << " for reading.\n";
    return false;
  }
 
  do {
    std::getline(fmesh, line);
  } while (!ends_with(line, "# number of mesh points"));
  nNodes = readInt(line);
 
  std::cout << m_className << "::Initialise:\n";
  std::cout << "    Read " << nNodes << " nodes from file " << mesh << ".\n";
  do {
    std::getline(fmesh, line);
  } while (line != "# Mesh point coordinates");
  double minx = 1e100, miny = 1e100, minz = 1e100, maxx = -1e100, maxy = -1e100,
         maxz = -1e100;
  for (int i = 0; i < nNodes; ++i) {
    Node newNode;
    fmesh >> newNode.x >> newNode.y >> newNode.z;
    newNode.x *= unit;
    newNode.y *= unit;
    newNode.z *= unit;
    nodes.push_back(newNode);
    minx = std::min(minx, newNode.x);
    maxx = std::max(maxx, newNode.x);
    miny = std::min(miny, newNode.y);
    maxy = std::max(maxy, newNode.y);
    minz = std::min(minz, newNode.z);
    maxz = std::max(maxz, newNode.z);
  }
  std::cout << minx << " < x < " << maxx << "\n";
  std::cout << miny << " < y < " << maxy << "\n";
  std::cout << minz << " < z < " << maxz << "\n";
 
  do {
    std::getline(fmesh, line);
  } while (line != "4 tet2 # type name");
  do {
    std::getline(fmesh, line);
  } while (!ends_with(line, "# number of elements"));
  nElements = readInt(line);
  elements.clear();
  std::cout << m_className << "::Initialise:\n";
  std::cout << "    Read " << nElements << " elements from file " << mesh
            << ".\n";
  std::getline(fmesh, line);
  // elements 6 & 7 are swapped due to differences in COMSOL and ANSYS
  // representation
  int perm[10] = {0, 1, 2, 3, 4, 5, 7, 6, 8, 9};
  for (int i = 0; i < nElements; ++i) {
    Element newElement;
    newElement.degenerate = false;
    for (int j = 0; j < 10; ++j) {
      fmesh >> newElement.emap[perm[j]];
    }
    elements.push_back(newElement);
  }
 
  do {
    std::getline(fmesh, line);
  } while (line != "# Geometric entity indices");
  for (int i = 0; i < nElements; ++i) {
    int domain;
    fmesh >> domain;
    elements[i].matmap = domain2material.count(domain) ? domain2material[domain]
                                                       : m_nMaterials - 1;
  }
  fmesh.close();
 
  std::map<Node, std::vector<int>, nodeCmp> nodeIdx;
  for (int i = 0; i < nNodes; ++i) {
    nodeIdx[nodes[i]].push_back(i);
  }
  std::cout << "Map size: " << nodeIdx.size() << std::endl;
 
  std::ifstream ffield;
  ffield.open(field.c_str(), std::ios::in);
  if (ffield.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open field potentials file " << field
              << " for reading.\n";
    return false;
  }
  do {
    std::getline(ffield, line);
  } while (line.substr(0, 81) !=
           "% x                       y                        z               "
           "         V (V)");
  {
    std::istringstream sline(line);
    std::string token;
    sline >> token;  // %
    sline >> token;  // x
    sline >> token;  // y
    sline >> token;  // z
    sline >> token;  // V
    sline >> token;  // (V)
    while (sline >> token) {
      std::cout << m_className << "::Initialise:\n";
      std::cout << "    Reading data for weighting field " << token << ".\n";
      nWeightingFields++;
      wfields.push_back(token);
      wfieldsOk.push_back(true);
      sline >> token;  // (V)
    }
  }
  for (int i = 0; i < nNodes; ++i) {
    Node tmp;
    ffield >> tmp.x >> tmp.y >> tmp.z >> tmp.v;
    tmp.x *= unit;
    tmp.y *= unit;
    tmp.z *= unit;
    for (int j = 0; j < nWeightingFields; ++j) {
      double w;
      ffield >> w;
      tmp.w.push_back(w);
    }
    int closest = -1;
    double closestDist = 1;
    const unsigned int nIdx = nodeIdx[tmp].size();
    // for (int j : nodeIdx[tmp]) {
    for (unsigned int k = 0; k < nIdx; ++k) {
      int j = nodeIdx[tmp][k];
      double dist = (tmp.x - nodes[j].x) * (tmp.x - nodes[j].x) +
                    (tmp.y - nodes[j].y) * (tmp.y - nodes[j].y) +
                    (tmp.z - nodes[j].z) * (tmp.z - nodes[j].z);
      if (dist < closestDist) {
        closestDist = dist;
        closest = j;
      }
    }
    if (closest == -1) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Could not match the node from field potentials file: "
                << tmp.x << " " << tmp.y << " " << tmp.z << "\n.";
      return false;
    }
    nodes[closest].v = tmp.v;
    nodes[closest].w = tmp.w;
  }
 
  m_ready = true;
 
  //  for (int i = 0; i < nNodes; ++i) {
  //    double ex, ey, ez, v;
  //    Medium* m;
  //    int status;
  //    ElectricField(nodes[i].x, nodes[i].y, nodes[i].z, ex, ey, ez, v, m,
  // status);
  //    std::cout << "Field at " << nodes[i].x << " " << nodes[i].y << " " <<
  // nodes[i].z << ": " << ex << " " << ey << " " << ez << " " << v << "\n";
  //  }
 
  // Establish the ranges.
  SetRange();
  UpdatePeriodicity();
  return true;
}

Referenced by ComponentComsol().

◆ IsInBoundingBox()

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

inlinevirtual

Reimplemented from Garfield::ComponentFieldMap.

Definition at line 33 of file ComponentComsol.hh.

                                                     {
    return x >= xMinBoundingBox && x <= xMaxBoundingBox &&
           y >= yMinBoundingBox && y <= yMaxBoundingBox &&
           z >= zMinBoundingBox && y <= zMaxBoundingBox;
  }

◆ SetWeightingField()

bool Garfield::ComponentComsol::SetWeightingField	(	std::string	file,
		std::string	label
	)

Definition at line 251 of file ComponentComsol.cc.

                                                                        {
  double unit = 100.0;  // m;
 
  if (!m_ready) {
    std::cerr << m_className << "::SetWeightingField:\n";
    std::cerr << "    No valid field map is present.\n";
    std::cerr << "    Weighting field cannot be added.\n";
    return false;
  }
 
  // Open the voltage list.
  std::ifstream ffield;
  ffield.open(field.c_str(), std::ios::in);
  if (ffield.fail()) {
    std::cerr << m_className << "::Initialise:\n";
    std::cerr << "    Could not open field potentials file " << field
              << " for reading.\n";
    return false;
  }
 
  // Check if a weighting field with the same label alm_ready exists.
  int iw = nWeightingFields;
  for (int i = nWeightingFields; i--;) {
    if (wfields[i] == label) {
      iw = i;
      break;
    }
  }
  if (iw == nWeightingFields) {
    ++nWeightingFields;
    wfields.resize(nWeightingFields);
    wfieldsOk.resize(nWeightingFields);
    for (int j = 0; j < nNodes; ++j) {
      nodes[j].w.resize(nWeightingFields);
    }
  } else {
    std::cout << m_className << "::SetWeightingField:\n";
    std::cout << "    Replacing existing weighting field " << label << ".\n";
  }
  wfields[iw] = label;
  wfieldsOk[iw] = false;
  std::map<Node, std::vector<int>, nodeCmp> nodeIdx;
  for (int i = 0; i < nNodes; ++i) {
    nodeIdx[nodes[i]].push_back(i);
  }
  std::cout << "Map size: " << nodeIdx.size() << std::endl;
 
  std::string line;
  do {
    std::getline(ffield, line);
  } while (line !=
           "% x                       y                        z               "
           "         V (V)");
  for (int i = 0; i < nNodes; ++i) {
    Node tmp;
    ffield >> tmp.x >> tmp.y >> tmp.z >> tmp.v;
    tmp.x *= unit;
    tmp.y *= unit;
    tmp.z *= unit;
    int closest = -1;
    double closestDist = 1;
    const unsigned int nIdx = nodeIdx[tmp].size();
    // for (int j : nodeIdx[tmp]) {
    for (unsigned int k = 0; k < nIdx; ++k) {
      int j = nodeIdx[tmp][k];
      double dist = (tmp.x - nodes[j].x) * (tmp.x - nodes[j].x) +
                    (tmp.y - nodes[j].y) * (tmp.y - nodes[j].y) +
                    (tmp.z - nodes[j].z) * (tmp.z - nodes[j].z);
      if (dist < closestDist) {
        closestDist = dist;
        closest = j;
      }
    }
    if (closest == -1) {
      std::cerr << m_className << "::Initialise:\n";
      std::cerr << "    Could not match the node from field potentials file: "
                << tmp.x << " " << tmp.y << " " << tmp.z << "\n.";
      return false;
    }
    nodes[closest].w[iw] = tmp.v;
  }
 
  return true;
}

◆ UpdatePeriodicity()

void Garfield::ComponentComsol::UpdatePeriodicity ( )

inlineprotectedvirtual

Verify periodicities.

Implements Garfield::ComponentFieldMap.

Definition at line 48 of file ComponentComsol.hh.

48{ UpdatePeriodicityCommon(); }

Garfield::ComponentFieldMap::UpdatePeriodicityCommon

void UpdatePeriodicityCommon()

Definition: ComponentFieldMap.cc:2008

Referenced by Initialise().

◆ WeightingField()

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

virtual

Calculate the weighting field at a given point and for a given electrode.

Parameters

x,y,z	coordinates [cm].
wx,wy,wz	components of the weighting field [1/cm].
label	name of the electrode

Implements Garfield::ComponentFieldMap.

Definition at line 450 of file ComponentComsol.cc.

                                                                         {
 
  // Initial values
  wx = wy = wz = 0;
 
  // Do not proceed if not properly initialised.
  if (!m_ready) return;
 
  // Look for the label.
  int iw = 0;
  bool found = false;
  for (int i = nWeightingFields; i--;) {
    if (wfields[i] == label) {
      iw = i;
      found = true;
      break;
    }
  }
 
  // Do not proceed if the requested weighting field does not exist.
  if (!found) return;
  // Check if the weighting field is properly initialised.
  if (!wfieldsOk[iw]) return;
 
  // Copy the coordinates.
  double x = xin, y = yin, z = zin;
 
  // Map the coordinates onto field map coordinates
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  if (m_warning) PrintWarning("WeightingField");
 
  // Find the element that contains this point.
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
  // Check if the point is in the mesh.
  if (imap < 0) return;
 
  if (m_debug) {
    PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, imap, 10, iw);
  }
 
  const Element& element = elements[imap];
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const Node& n4 = nodes[element.emap[4]];
  const Node& n5 = nodes[element.emap[5]];
  const Node& n6 = nodes[element.emap[6]];
  const Node& n7 = nodes[element.emap[7]];
  const Node& n8 = nodes[element.emap[8]];
  const Node& n9 = nodes[element.emap[9]];
  // Tetrahedral field
  wx = -(n0.w[iw] * (4 * t1 - 1) * jac[0][1] +
         n1.w[iw] * (4 * t2 - 1) * jac[1][1] +
         n2.w[iw] * (4 * t3 - 1) * jac[2][1] +
         n3.w[iw] * (4 * t4 - 1) * jac[3][1] +
         n4.w[iw] * (4 * t2 * jac[0][1] + 4 * t1 * jac[1][1]) +
         n5.w[iw] * (4 * t3 * jac[0][1] + 4 * t1 * jac[2][1]) +
         n6.w[iw] * (4 * t4 * jac[0][1] + 4 * t1 * jac[3][1]) +
         n7.w[iw] * (4 * t3 * jac[1][1] + 4 * t2 * jac[2][1]) +
         n8.w[iw] * (4 * t4 * jac[1][1] + 4 * t2 * jac[3][1]) +
         n9.w[iw] * (4 * t4 * jac[2][1] + 4 * t3 * jac[3][1])) /
       det;
 
  wy = -(n0.w[iw] * (4 * t1 - 1) * jac[0][2] +
         n1.w[iw] * (4 * t2 - 1) * jac[1][2] +
         n2.w[iw] * (4 * t3 - 1) * jac[2][2] +
         n3.w[iw] * (4 * t4 - 1) * jac[3][2] +
         n4.w[iw] * (4 * t2 * jac[0][2] + 4 * t1 * jac[1][2]) +
         n5.w[iw] * (4 * t3 * jac[0][2] + 4 * t1 * jac[2][2]) +
         n6.w[iw] * (4 * t4 * jac[0][2] + 4 * t1 * jac[3][2]) +
         n7.w[iw] * (4 * t3 * jac[1][2] + 4 * t2 * jac[2][2]) +
         n8.w[iw] * (4 * t4 * jac[1][2] + 4 * t2 * jac[3][2]) +
         n9.w[iw] * (4 * t4 * jac[2][2] + 4 * t3 * jac[3][2])) /
       det;
 
  wz = -(n0.w[iw] * (4 * t1 - 1) * jac[0][3] +
         n1.w[iw] * (4 * t2 - 1) * jac[1][3] +
         n2.w[iw] * (4 * t3 - 1) * jac[2][3] +
         n3.w[iw] * (4 * t4 - 1) * jac[3][3] +
         n4.w[iw] * (4 * t2 * jac[0][3] + 4 * t1 * jac[1][3]) +
         n5.w[iw] * (4 * t3 * jac[0][3] + 4 * t1 * jac[2][3]) +
         n6.w[iw] * (4 * t4 * jac[0][3] + 4 * t1 * jac[3][3]) +
         n7.w[iw] * (4 * t3 * jac[1][3] + 4 * t2 * jac[2][3]) +
         n8.w[iw] * (4 * t4 * jac[1][3] + 4 * t2 * jac[3][3]) +
         n9.w[iw] * (4 * t4 * jac[2][3] + 4 * t3 * jac[3][3])) /
       det;
 
  // Transform field to global coordinates
  UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
}

◆ WeightingPotential()

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

virtual

Implements Garfield::ComponentFieldMap.

Definition at line 548 of file ComponentComsol.cc.

                                                                   {
 
  // Do not proceed if not properly initialised.
  if (!m_ready) return 0.;
 
  // Look for the label.
  int iw = 0;
  bool found = false;
  for (int i = nWeightingFields; i--;) {
    if (wfields[i] == label) {
      iw = i;
      found = true;
      break;
    }
  }
 
  // Do not proceed if the requested weighting field does not exist.
  if (!found) return 0.;
  // Check if the weighting field is properly initialised.
  if (!wfieldsOk[iw]) return 0.;
 
  // Copy the coordinates.
  double x = xin, y = yin, z = zin;
 
  // Map the coordinates onto field map coordinates.
  bool xmirr, ymirr, zmirr;
  double rcoordinate, rotation;
  MapCoordinates(x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  if (m_warning) PrintWarning("WeightingPotential");
 
  // Find the element that contains this point.
  double t1, t2, t3, t4, jac[4][4], det;
  const int imap = FindElement13(x, y, z, t1, t2, t3, t4, jac, det);
  if (imap < 0) return 0.;
 
  if (m_debug) {
    PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, imap, 10, iw);
  }
 
  const Element& element = elements[imap];
  const Node& n0 = nodes[element.emap[0]];
  const Node& n1 = nodes[element.emap[1]];
  const Node& n2 = nodes[element.emap[2]];
  const Node& n3 = nodes[element.emap[3]];
  const Node& n4 = nodes[element.emap[4]];
  const Node& n5 = nodes[element.emap[5]];
  const Node& n6 = nodes[element.emap[6]];
  const Node& n7 = nodes[element.emap[7]];
  const Node& n8 = nodes[element.emap[8]];
  const Node& n9 = nodes[element.emap[9]];
  // Tetrahedral field
  return n0.w[iw] * t1 * (2 * t1 - 1) + n1.w[iw] * t2 * (2 * t2 - 1) +
         n2.w[iw] * t3 * (2 * t3 - 1) + n3.w[iw] * t4 * (2 * t4 - 1) +
         4 * n4.w[iw] * t1 * t2 + 4 * n5.w[iw] * t1 * t3 +
         4 * n6.w[iw] * t1 * t4 + 4 * n7.w[iw] * t2 * t3 +
         4 * n8.w[iw] * t2 * t4 + 4 * n9.w[iw] * t3 * t4;
}

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

Classes

Public Member Functions

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ComponentComsol() [1/2]

◆ ComponentComsol() [2/2]

◆ ~ComponentComsol()

Member Function Documentation

◆ ElectricField() [1/2]

◆ ElectricField() [2/2]

◆ GetAspectRatio()

◆ GetElementVolume()

◆ GetMedium()

◆ Initialise()

◆ IsInBoundingBox()

◆ SetWeightingField()

◆ UpdatePeriodicity()

◆ WeightingField()

◆ WeightingPotential()