Component for importing field maps computed by Elmer. More...

#include <ComponentElmer.hh>

Inheritance diagram for Garfield::ComponentElmer:

Public Member Functions
	ComponentElmer ()
	Default constructor.

	ComponentElmer (const std::string &header, const std::string &elist, const std::string &nlist, const std::string &mplist, const std::string &volt, const std::string &unit)
	Constructor with a set of field map files, see Initialise().

	~ComponentElmer ()
	Destructor.

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

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

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

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

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

bool	Initialise (const std::string &header="mesh.header", const std::string &elist="mesh.elements", const std::string &nlist="mesh.nodes", const std::string &mplist="dielectrics.dat", const std::string &volt="out.result", const std::string &unit="cm")

bool	SetWeightingField (std::string prnsol, std::string label)
	Import a list of voltages to be used as weighting field.

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.

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

double	GetElementVolume (const unsigned int i) override

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

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

virtual void	Reset ()=0
	Reset the component.

virtual void	UpdatePeriodicity ()=0
	Verify periodicities.

Additional Inherited Members
Protected Attributes inherited from Garfield::ComponentFieldMap
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.

Detailed Description

Component for importing field maps computed by Elmer.

Definition at line 12 of file ComponentElmer.hh.

Constructor & Destructor Documentation

◆ ComponentElmer() [1/2]

Garfield::ComponentElmer::ComponentElmer ( )

Default constructor.

Definition at line 28 of file ComponentElmer.cc.

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

◆ ComponentElmer() [2/2]

Garfield::ComponentElmer::ComponentElmer	(	const std::string &	header,
		const std::string &	elist,
		const std::string &	nlist,
		const std::string &	mplist,
		const std::string &	volt,
		const std::string &	unit
	)

Constructor with a set of field map files, see Initialise().

Definition at line 32 of file ComponentElmer.cc.

    : ComponentFieldMap() {
  m_className = "ComponentElmer";
  Initialise(header, elist, nlist, mplist, volt, unit);
}

◆ ~ComponentElmer()

Garfield::ComponentElmer::~ComponentElmer ( )

inline

Destructor.

Definition at line 21 of file ComponentElmer.hh.

21{}

Member Function Documentation

◆ ElectricField() [1/2]

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

overridevirtual

Calculate the drift field [V/cm] and potential [V] at (x, y, z).

Implements Garfield::ComponentBase.

Definition at line 541 of file ComponentElmer.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 = nullptr;
 
  // 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    Point (" << x << ", "
                << y << ", " << z << ") is not in the mesh.\n";
    }
    status = -6;
    return;
  }
 
  const Element& element = elements[imap];
  if (m_debug) {
    PrintElement("ElectricField", x, y, z, t1, t2, t3, t4, element, 10);
  }
  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]];
  // Shorthands.
  const double fourt1 = 4 * t1;
  const double fourt2 = 4 * t2;
  const double fourt3 = 4 * t3;
  const double fourt4 = 4 * t4;
  const double invdet = 1. / det;
  // 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 * (fourt1 - 1) * jac[0][1] + n1.v * (fourt2 - 1) * jac[1][1] +
         n2.v * (fourt3 - 1) * jac[2][1] + n3.v * (fourt4 - 1) * jac[3][1] +
         n4.v * (fourt2 * jac[0][1] + fourt1 * jac[1][1]) +
         n5.v * (fourt3 * jac[0][1] + fourt1 * jac[2][1]) +
         n6.v * (fourt4 * jac[0][1] + fourt1 * jac[3][1]) +
         n7.v * (fourt3 * jac[1][1] + fourt2 * jac[2][1]) +
         n8.v * (fourt4 * jac[1][1] + fourt2 * jac[3][1]) +
         n9.v * (fourt4 * jac[2][1] + fourt3 * jac[3][1])) *
       invdet;
  ey = -(n0.v * (fourt1 - 1) * jac[0][2] + n1.v * (fourt2 - 1) * jac[1][2] +
         n2.v * (fourt3 - 1) * jac[2][2] + n3.v * (fourt4 - 1) * jac[3][2] +
         n4.v * (fourt2 * jac[0][2] + fourt1 * jac[1][2]) +
         n5.v * (fourt3 * jac[0][2] + fourt1 * jac[2][2]) +
         n6.v * (fourt4 * jac[0][2] + fourt1 * jac[3][2]) +
         n7.v * (fourt3 * jac[1][2] + fourt2 * jac[2][2]) +
         n8.v * (fourt4 * jac[1][2] + fourt2 * jac[3][2]) +
         n9.v * (fourt4 * jac[2][2] + fourt3 * jac[3][2])) *
       invdet;
  ez = -(n0.v * (fourt1 - 1) * jac[0][3] + n1.v * (fourt2 - 1) * jac[1][3] +
         n2.v * (fourt3 - 1) * jac[2][3] + n3.v * (fourt4 - 1) * jac[3][3] +
         n4.v * (fourt2 * jac[0][3] + fourt1 * jac[1][3]) +
         n5.v * (fourt3 * jac[0][3] + fourt1 * jac[2][3]) +
         n6.v * (fourt4 * jac[0][3] + fourt1 * jac[3][3]) +
         n7.v * (fourt3 * jac[1][3] + fourt2 * jac[2][3]) +
         n8.v * (fourt4 * jac[1][3] + fourt2 * jac[3][3]) +
         n9.v * (fourt4 * jac[2][3] + fourt3 * jac[3][3])) *
       invdet;
 
  // Transform field to global coordinates
  UnmapFields(ex, ey, ez, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
 
  // Drift medium?
  const Material& mat = materials[element.matmap];
  if (m_debug) {
    std::cout << m_className << "::ElectricField:\n    Material "
              << element.matmap << ", drift flag " << mat.driftmedium << ".\n";
  }
  m = mat.medium;
  status = -5;
  if (mat.driftmedium && m && m->IsDriftable()) status = 0;
}

◆ ElectricField() [2/2]

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

overridevirtual

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

Definition at line 534 of file ComponentElmer.cc.

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

Referenced by ElectricField().

◆ GetAspectRatio()

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

overrideprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 869 of file ComponentElmer.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::ComponentElmer::GetElementVolume ( const unsigned int i )

overrideprotectedvirtual

Implements Garfield::ComponentFieldMap.

Definition at line 848 of file ComponentElmer.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::ComponentElmer::GetMedium	(	const double	x,
		const double	y,
		const double	z
	)

overridevirtual

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

Reimplemented from Garfield::ComponentBase.

Definition at line 806 of file ComponentElmer.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    Point (" << x << ", " << y
                << ", " << z << ") is 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    Point (" << x << ", " << y
                << ", " << z << ") has out of range material number " << imap
                << ".\n";
    }
    return nullptr;
  }
 
  if (m_debug) PrintElement("GetMedium", x, y, z, t1, t2, t3, t4, element, 10);
 
  return materials[element.matmap].medium;
}

◆ Initialise()

bool Garfield::ComponentElmer::Initialise	(	const std::string &	header = `"mesh.header"`,
		const std::string &	elist = `"mesh.elements"`,
		const std::string &	nlist = `"mesh.nodes"`,
		const std::string &	mplist = `"dielectrics.dat"`,
		const std::string &	volt = `"out.result"`,
		const std::string &	unit = `"cm"`
	)

Import a field map from a set of files.

Parameters

header	name of the header file (contains the number of elements and nodes).
elist	name of the file that contains the list of mesh elements
nlist	name of the file that contains the list of mesh nodes
mplist	name of the file that contains the material properties
volt	output of the field solver (list of voltages)
unit	length unit to be used

Definition at line 42 of file ComponentElmer.cc.

                                                       {
  const std::string hdr = m_className + "::Initialise:";
  m_debug = false;
  m_ready = false;
  m_warning = false;
  m_nWarnings = 0;
 
  // Keep track of the success.
  bool ok = true;
 
  // Buffer for reading
  const int size = 100;
  char line[size];
 
  // Open the header.
  std::ifstream fheader;
  fheader.open(header.c_str(), std::ios::in);
  if (fheader.fail()) {
    PrintErrorOpeningFile(hdr, "header", header);
  }
 
  // Temporary variables for use in file reading
  char* token = NULL;
  bool readerror = false;
  bool readstop = false;
  int il = 0;
 
  // Read the header to get the number of nodes and elements.
  fheader.getline(line, size, '\n');
  token = strtok(line, " ");
  nNodes = ReadInteger(token, 0, readerror);
  token = strtok(NULL, " ");
  nElements = ReadInteger(token, 0, readerror);
  std::cout << hdr << "\n    Read " << nNodes << " nodes and " << nElements
            << " elements from file " << header << ".\n";
  if (readerror) {
    PrintErrorReadingFile(hdr, header, il);
    fheader.close();
    return false;
  }
 
  // Close the header file.
  fheader.close();
 
  // Open the nodes list.
  std::ifstream fnodes;
  fnodes.open(nlist.c_str(), std::ios::in);
  if (fnodes.fail()) {
    PrintErrorOpeningFile(hdr, "nodes", nlist);
  }
 
  // Check the value of the unit.
  double funit;
  if (strcmp(unit.c_str(), "mum") == 0 || strcmp(unit.c_str(), "micron") == 0 ||
      strcmp(unit.c_str(), "micrometer") == 0) {
    funit = 0.0001;
  } else if (strcmp(unit.c_str(), "mm") == 0 ||
             strcmp(unit.c_str(), "millimeter") == 0) {
    funit = 0.1;
  } else if (strcmp(unit.c_str(), "cm") == 0 ||
             strcmp(unit.c_str(), "centimeter") == 0) {
    funit = 1.0;
  } else if (strcmp(unit.c_str(), "m") == 0 ||
             strcmp(unit.c_str(), "meter") == 0) {
    funit = 100.0;
  } else {
    std::cerr << hdr << " Unknown length unit " << unit << ".\n";
    ok = false;
    funit = 1.0;
  }
  if (m_debug) std::cout << hdr << " Unit scaling factor = " << funit << ".\n";
 
  // Read the nodes from the file.
  for (il = 0; il < nNodes; il++) {
    // Get a line from the nodes file.
    fnodes.getline(line, size, '\n');
 
    // Ignore the first two characters.
    token = strtok(line, " ");
    token = strtok(NULL, " ");
 
    // Get the node coordinates.
    token = strtok(NULL, " ");
    double xnode = ReadDouble(token, -1, readerror);
    token = strtok(NULL, " ");
    double ynode = ReadDouble(token, -1, readerror);
    token = strtok(NULL, " ");
    double znode = ReadDouble(token, -1, readerror);
    if (readerror) {
      PrintErrorReadingFile(hdr, nlist, il);
      fnodes.close();
      return false;
    }
 
    // Set up and create a new node.
    Node newNode;
    newNode.w.clear();
    newNode.x = xnode * funit;
    newNode.y = ynode * funit;
    newNode.z = znode * funit;
    nodes.push_back(std::move(newNode));
  }
 
  // Close the nodes file.
  fnodes.close();
 
  // Open the potential file.
  std::ifstream fvolt;
  fvolt.open(volt.c_str(), std::ios::in);
  if (fvolt.fail()) {
    PrintErrorOpeningFile(hdr, "potentials", volt);
  }
 
  // Reset the line counter.
  il = 1;
 
  // Read past the header.
  while (!readstop && fvolt.getline(line, size, '\n')) {
    token = strtok(line, " ");
    if (strcmp(token, "Perm:") == 0) readstop = true;
    il++;
  }
 
  // Should have stopped: if not, print error message.
  if (!readstop) {
    std::cerr << hdr << "\n    Error reading past header of potentials file "
              << volt << ".\n";
    fvolt.close();
    return false;
  }
 
  // Read past the permutation map (number of lines = nNodes).
  for (int tl = 0; tl < nNodes; tl++) {
    fvolt.getline(line, size, '\n');
    il++;
  }
 
  // Read the potentials.
  for (int tl = 0; tl < nNodes; tl++) {
    fvolt.getline(line, size, '\n');
    token = strtok(line, " ");
    double v = ReadDouble(token, -1, readerror);
    if (readerror) {
      PrintErrorReadingFile(hdr, volt, il);
      fvolt.close();
      return false;
    }
    // Place the voltage in its appropriate node.
    nodes[tl].v = v;
  }
 
  // Close the potentials file.
  fvolt.close();
 
  // Open the materials file.
  std::ifstream fmplist;
  fmplist.open(mplist.c_str(), std::ios::in);
  if (fmplist.fail()) {
    PrintErrorOpeningFile(hdr, "materials", mplist);
  }
 
  // Read the dielectric constants from the materials file.
  fmplist.getline(line, size, '\n');
  token = strtok(line, " ");
  if (readerror) {
    std::cerr << hdr << "\n    Error reading number of materials from "
              << mplist << ".\n";
    fmplist.close();
    ok = false;
    return false;
  }
  m_nMaterials = ReadInteger(token, 0, readerror);
  materials.resize(m_nMaterials);
  for (unsigned int i = 0; i < m_nMaterials; ++i) {
    materials[i].ohm = -1;
    materials[i].eps = -1;
    materials[i].medium = nullptr;
  }
  for (il = 2; il < ((int)m_nMaterials + 2); il++) {
    fmplist.getline(line, size, '\n');
    token = strtok(line, " ");
    ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    double dc = ReadDouble(token, -1.0, readerror);
    if (readerror) {
      PrintErrorReadingFile(hdr, mplist, il);
      fmplist.close();
      ok = false;
      return false;
    }
    materials[il - 2].eps = dc;
    std::cout << hdr << "\n    Set material " << il - 2 << " of "
              << m_nMaterials << " to eps " << dc << ".\n";
  }
 
  // Close the materials file.
  fmplist.close();
 
  // Find the lowest epsilon, check for eps = 0, set default drift media.
  double epsmin = -1.;
  unsigned int iepsmin = 0;
  for (unsigned int imat = 0; imat < m_nMaterials; ++imat) {
    if (materials[imat].eps < 0) continue;
    if (materials[imat].eps == 0) {
      std::cerr << hdr << "\n    Material " << imat
                << " has been assigned a permittivity equal to zero in\n    "
                << mplist << ".\n";
      ok = false;
    } else if (epsmin < 0. || epsmin > materials[imat].eps) {
      epsmin = materials[imat].eps;
      iepsmin = imat;
    }
  }
 
  if (epsmin < 0.) {
    std::cerr << hdr << "\n    No material with positive permittivity found \n"
              << "    in material list " << mplist << ".\n";
    ok = false;
  } else {
    for (unsigned int imat = 0; imat < m_nMaterials; ++imat) {
      materials[imat].driftmedium = imat == iepsmin ? true : false;
    }
  }
 
  // Open the elements file.
  std::ifstream felems;
  felems.open(elist.c_str(), std::ios::in);
  if (felems.fail()) {
    PrintErrorOpeningFile(hdr, "elements", elist);
  }
 
  // Read the elements and their material indices.
  elements.clear();
  int highestnode = 0;
  Element newElement;
  for (il = 0; il < nElements; il++) {
    // Get a line
    felems.getline(line, size, '\n');
 
    // Split into tokens.
    token = strtok(line, " ");
    // Read the 2nd-order element
    // Note: Ordering of Elmer elements can be described in the
    // ElmerSolver manual (appendix D. at the time of this comment)
    // If the order read below is compared to the shape functions used
    // eg. in ElectricField, the order is wrong, but note at the
    // end of this function the order of elements 5,6,7 will change to
    // 7,5,6 when actually recorded in newElement.emap to correct for this
    token = strtok(NULL, " ");
    int imat = ReadInteger(token, -1, readerror) - 1;
    token = strtok(NULL, " ");
    token = strtok(NULL, " ");
    int in0 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in1 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in2 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in3 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in4 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in5 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in6 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in7 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in8 = ReadInteger(token, -1, readerror);
    token = strtok(NULL, " ");
    int in9 = ReadInteger(token, -1, readerror);
 
    if (m_debug && il < 10) {
      std::cout << "    Read nodes " << in0 << ", " << in1 << ", " << in2
                << ", " << in3 << ", ... from element " << il + 1 << " of "
                << nElements << " with mat " << imat << ".\n";
    }
 
    // Check synchronisation.
    if (readerror) {
      PrintErrorReadingFile(hdr, elist, il);
      felems.close();
      ok = false;
      return false;
    }
 
    // Check the material number and ensure that epsilon is non-negative.
    if (imat < 0 || imat > (int)m_nMaterials) {
      std::cerr << hdr << "\n    Out-of-range material number on file " << elist
                << " (line " << il << ").\n";
      std::cerr << "    Element: " << il << ", material: " << imat << "\n";
      std::cerr << "    nodes: (" << in0 << ", " << in1 << ", " << in2 << ", "
                << in3 << ", " << in4 << ", " << in5 << ", " << in6 << ", "
                << in7 << ", " << in8 << ", " << in9 << ")\n";
      ok = false;
    }
    if (materials[imat].eps < 0) {
      std::cerr << hdr << "\n    Element " << il << " in element list " << elist
                << "\n    uses material " << imat
                << " which has not been assigned a positive permittivity in "
                << mplist << ".\n";
      ok = false;
    }
 
    // Check the node numbers.
    if (in0 < 1 || in1 < 1 || in2 < 1 || in3 < 1 || in4 < 1 || in5 < 1 ||
        in6 < 1 || in7 < 1 || in8 < 1 || in9 < 1) {
      std::cerr << hdr << "\n    Found a node number < 1 on file " << elist
                << " (line " << il << ").\n    Element: " << il
                << ", material: " << imat << "\n    nodes: (" << in0 << ", "
                << in1 << ", " << in2 << ", " << in3 << ", " << in4 << ", "
                << in5 << ", " << in6 << ", " << in7 << ", " << in8 << ", "
                << in9 << ")\n";
      ok = false;
    }
    if (in0 > highestnode) highestnode = in0;
    if (in1 > highestnode) highestnode = in1;
    if (in2 > highestnode) highestnode = in2;
    if (in3 > highestnode) highestnode = in3;
    if (in4 > highestnode) highestnode = in4;
    if (in5 > highestnode) highestnode = in5;
    if (in6 > highestnode) highestnode = in6;
    if (in7 > highestnode) highestnode = in7;
    if (in8 > highestnode) highestnode = in8;
    if (in9 > highestnode) highestnode = in9;
 
    // These elements must not be degenerate.
    if (in0 == in1 || in0 == in2 || in0 == in3 || in0 == in4 || in0 == in5 ||
        in0 == in6 || in0 == in7 || in0 == in8 || in0 == in9 || in1 == in2 ||
        in1 == in3 || in1 == in4 || in1 == in5 || in1 == in6 || in1 == in7 ||
        in1 == in8 || in1 == in9 || in2 == in3 || in2 == in4 || in2 == in5 ||
        in2 == in6 || in2 == in7 || in2 == in8 || in2 == in9 || in3 == in4 ||
        in3 == in5 || in3 == in6 || in3 == in7 || in3 == in8 || in3 == in9 ||
        in4 == in5 || in4 == in6 || in4 == in7 || in4 == in8 || in4 == in9 ||
        in5 == in6 || in5 == in7 || in5 == in8 || in5 == in9 || in6 == in7 ||
        in6 == in8 || in6 == in9 || in7 == in8 || in7 == in9 || in8 == in9) {
      std::cerr << hdr << "\n    Element " << il << " of file " << elist
                << " is degenerate,\n"
                << "    no such elements are allowed in this type of map.\n";
      ok = false;
    }
 
    newElement.degenerate = false;
 
    // Store the material reference.
    newElement.matmap = imat;
 
    // Node references
    newElement.emap[0] = in0 - 1;
    newElement.emap[1] = in1 - 1;
    newElement.emap[2] = in2 - 1;
    newElement.emap[3] = in3 - 1;
    newElement.emap[4] = in4 - 1;
    newElement.emap[7] = in5 - 1;
    newElement.emap[5] = in6 - 1;
    newElement.emap[6] = in7 - 1;
    newElement.emap[8] = in8 - 1;
    newElement.emap[9] = in9 - 1;
    elements.push_back(newElement);
  }
 
  // Close the elements file.
  felems.close();
 
  // Set the ready flag.
  if (ok) {
    m_ready = true;
  } else {
    std::cerr << hdr << "\n    Field map could not be "
              << "read and cannot be interpolated.\n";
    return false;
  }
 
  std::cout << hdr << " Finished.\n";
 
  // Remove weighting fields (if any).
  wfields.clear();
  wfieldsOk.clear();
  nWeightingFields = 0;
 
  // Establish the ranges.
  SetRange();
  UpdatePeriodicity();
  return true;
}

Referenced by ComponentElmer().

◆ SetWeightingField()

bool Garfield::ComponentElmer::SetWeightingField	(	std::string	prnsol,
		std::string	label
	)

Import a list of voltages to be used as weighting field.

Definition at line 433 of file ComponentElmer.cc.

                                                                       {
  const std::string hdr = m_className + "::SetWeightingField:";
  if (!m_ready) {
    PrintNotReady("SetWeightingField");
    std::cerr << "    Weighting field cannot be added.\n";
    return false;
  }
 
  // Keep track of the success.
  bool ok = true;
 
  // Open the voltage list.
  std::ifstream fwvolt;
  fwvolt.open(wvolt.c_str(), std::ios::in);
  if (fwvolt.fail()) {
    PrintErrorOpeningFile(hdr, "potential", wvolt);
    return false;
  }
 
  // Check if a weighting field with the same label already 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 = nNodes; j--;) {
      nodes[j].w.resize(nWeightingFields);
    }
  } else {
    std::cout << hdr << "\n    Replacing existing weighting field " << label
              << ".\n";
  }
  wfields[iw] = label;
  wfieldsOk[iw] = false;
 
  // Temporary variables for use in file reading
  const int size = 100;
  char line[size];
  char* token = NULL;
  bool readerror = false;
  bool readstop = false;
  int il = 1;
 
  // Read past the header.
  while (!readstop && fwvolt.getline(line, size, '\n')) {
    token = strtok(line, " ");
    if (strcmp(token, "Perm:") == 0) readstop = true;
    il++;
  }
 
  // Should have stopped: if not, print error message.
  if (!readstop) {
    std::cerr << hdr << "\n    Error reading past header of potentials file "
              << wvolt << ".\n";
    fwvolt.close();
    ok = false;
    return false;
  }
 
  // Read past the permutation map (number of lines = nNodes).
  for (int tl = 0; tl < nNodes; tl++) {
    fwvolt.getline(line, size, '\n');
    il++;
  }
 
  // Read the potentials.
  for (int tl = 0; tl < nNodes; tl++) {
    double v;
    fwvolt.getline(line, size, '\n');
    token = strtok(line, " ");
    v = ReadDouble(token, -1, readerror);
    if (readerror) {
      PrintErrorReadingFile(hdr, wvolt, il);
      fwvolt.close();
      ok = false;
      return false;
    }
    // Place the weighting potential at its appropriate node and index.
    nodes[tl].w[iw] = v;
  }
 
  // Close the potentials file.
  fwvolt.close();
  std::cout << hdr << "\n    Read potentials from file " << wvolt << ".\n";
 
  // Set the ready flag.
  wfieldsOk[iw] = ok;
  if (!ok) {
    std::cerr << hdr << "\n    Field map could not "
              << "be read and cannot be interpolated.\n";
    return false;
  }
  return true;
}

◆ UpdatePeriodicity()

void Garfield::ComponentElmer::UpdatePeriodicity ( )

inlineoverrideprotectedvirtual

Verify periodicities.

Implements Garfield::ComponentBase.

Definition at line 57 of file ComponentElmer.hh.

57{ UpdatePeriodicityCommon(); }

Garfield::ComponentFieldMap::UpdatePeriodicityCommon

void UpdatePeriodicityCommon()

Definition: ComponentFieldMap.cc:1753

Referenced by Initialise().

◆ WeightingField()

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

overridevirtual

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

Reimplemented from Garfield::ComponentBase.

Definition at line 646 of file ComponentElmer.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;
 
  const Element& element = elements[imap];
  if (m_debug) {
    PrintElement("WeightingField", x, y, z, t1, t2, t3, t4, element, 10, iw);
  }
  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]];
  // Shorthands.
  const double fourt1 = 4 * t1;
  const double fourt2 = 4 * t2;
  const double fourt3 = 4 * t3;
  const double fourt4 = 4 * t4;
  const double invdet = 1. / det;
  // Tetrahedral field
  wx = -(n0.w[iw] * (fourt1 - 1) * jac[0][1] +
         n1.w[iw] * (fourt2 - 1) * jac[1][1] +
         n2.w[iw] * (fourt3 - 1) * jac[2][1] +
         n3.w[iw] * (fourt4 - 1) * jac[3][1] +
         n4.w[iw] * (fourt2 * jac[0][1] + fourt1 * jac[1][1]) +
         n5.w[iw] * (fourt3 * jac[0][1] + fourt1 * jac[2][1]) +
         n6.w[iw] * (fourt4 * jac[0][1] + fourt1 * jac[3][1]) +
         n7.w[iw] * (fourt3 * jac[1][1] + fourt2 * jac[2][1]) +
         n8.w[iw] * (fourt4 * jac[1][1] + fourt2 * jac[3][1]) +
         n9.w[iw] * (fourt4 * jac[2][1] + fourt3 * jac[3][1])) *
       invdet;
  wy = -(n0.w[iw] * (fourt1 - 1) * jac[0][2] +
         n1.w[iw] * (fourt2 - 1) * jac[1][2] +
         n2.w[iw] * (fourt3 - 1) * jac[2][2] +
         n3.w[iw] * (fourt4 - 1) * jac[3][2] +
         n4.w[iw] * (fourt2 * jac[0][2] + fourt1 * jac[1][2]) +
         n5.w[iw] * (fourt3 * jac[0][2] + fourt1 * jac[2][2]) +
         n6.w[iw] * (fourt4 * jac[0][2] + fourt1 * jac[3][2]) +
         n7.w[iw] * (fourt3 * jac[1][2] + fourt2 * jac[2][2]) +
         n8.w[iw] * (fourt4 * jac[1][2] + fourt2 * jac[3][2]) +
         n9.w[iw] * (fourt4 * jac[2][2] + fourt3 * jac[3][2])) *
       invdet;
  wz = -(n0.w[iw] * (fourt1 - 1) * jac[0][3] +
         n1.w[iw] * (fourt2 - 1) * jac[1][3] +
         n2.w[iw] * (fourt3 - 1) * jac[2][3] +
         n3.w[iw] * (fourt4 - 1) * jac[3][3] +
         n4.w[iw] * (fourt2 * jac[0][3] + fourt1 * jac[1][3]) +
         n5.w[iw] * (fourt3 * jac[0][3] + fourt1 * jac[2][3]) +
         n6.w[iw] * (fourt4 * jac[0][3] + fourt1 * jac[3][3]) +
         n7.w[iw] * (fourt3 * jac[1][3] + fourt2 * jac[2][3]) +
         n8.w[iw] * (fourt4 * jac[1][3] + fourt2 * jac[3][3]) +
         n9.w[iw] * (fourt4 * jac[2][3] + fourt3 * jac[3][3])) *
       invdet;
 
  // Transform field to global coordinates
  UnmapFields(wx, wy, wz, x, y, z, xmirr, ymirr, zmirr, rcoordinate, rotation);
}

◆ WeightingPotential()

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

overridevirtual

Calculate the weighting potential at a given point.

Parameters

x,y,z	coordinates [cm].
label	name of the electrode.

Returns: weighting potential [dimensionless].

Reimplemented from Garfield::ComponentBase.

Definition at line 746 of file ComponentElmer.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.;
 
  const Element& element = elements[imap];
  if (m_debug) {
    PrintElement("WeightingPotential", x, y, z, t1, t2, t3, t4, element, 10,
                 iw);
  }
  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:

Public Member Functions

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ComponentElmer() [1/2]

◆ ComponentElmer() [2/2]

◆ ~ComponentElmer()

Member Function Documentation

◆ ElectricField() [1/2]

◆ ElectricField() [2/2]

◆ GetAspectRatio()

◆ GetElementVolume()

◆ GetMedium()

◆ Initialise()

◆ SetWeightingField()

◆ UpdatePeriodicity()

◆ WeightingField()

◆ WeightingPotential()