Garfield::ComponentNeBem3d Class Reference

Interface to neBEM. More...

#include <ComponentNeBem3d.hh>

Inheritance diagram for Garfield::ComponentNeBem3d:

Public Member Functions
	ComponentNeBem3d ()
	Constructor.
	~ComponentNeBem3d ()
	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).
bool	GetVoltageRange (double &vmin, double &vmax) override
	Calculate the voltage range [V].
unsigned int	GetNumberOfPrimitives () const
bool	GetPrimitive (const unsigned int i, double &a, double &b, double &c, std::vector< double > &xv, std::vector< double > &yv, std::vector< double > &zv, int &interface, double &v, double &q, double &lambda) const
unsigned int	GetNumberOfElements () const
bool	GetElement (const unsigned int i, std::vector< double > &xv, std::vector< double > &yv, std::vector< double > &zv, int &interface, double &bc, double &lambda) const
void	Reset () override
	Reset the component.
void	UpdatePeriodicity () override
	Verify periodicities.
bool	Initialise ()
void	SetTargetElementSize (const double length)
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)

Additional Inherited Members
virtual void	Reset ()=0
	Reset the component.
virtual void	UpdatePeriodicity ()=0
	Verify periodicities.
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

Interface to neBEM.

Definition at line 9 of file ComponentNeBem3d.hh.

Constructor & Destructor Documentation

ComponentNeBem3d()

Garfield::ComponentNeBem3d::ComponentNeBem3d

(

)

Constructor.

Definition at line 377 of file ComponentNeBem3d.cc.

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

~ComponentNeBem3d()

Garfield::ComponentNeBem3d::~ComponentNeBem3d ( )

inline

Destructor.

Definition at line 14 of file ComponentNeBem3d.hh.

{}

Member Function Documentation

ElectricField() [1/2]

void Garfield::ComponentNeBem3d::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 381 of file ComponentNeBem3d.cc.

                                                  {
  ex = ey = ez = v = 0.;
  status = 0;
  // Check if the requested point is inside a medium
  m = GetMedium(x, y, z);
  if (!m) {
    status = -6;
    return;
  }
 
  if (!m_ready) {
    if (!Initialise()) {
      std::cerr << m_className << "::ElectricField: Initialisation failed.\n";
      status = -11;
      return;
    }
    m_ready = true;
  }
}

ElectricField() [2/2]

void Garfield::ComponentNeBem3d::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 404 of file ComponentNeBem3d.cc.

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

Referenced by ElectricField().

GetElement()

bool Garfield::ComponentNeBem3d::GetElement	(	const unsigned int	i,
		std::vector< double > &	xv,
		std::vector< double > &	yv,
		std::vector< double > &	zv,
		int &	interface,
		double &	bc,
		double &	lambda
	)		const

Definition at line 2499 of file ComponentNeBem3d.cc.

                                                                    {
 
  if (i >= m_elements.size()) {
    std::cerr << m_className << "::GetElement: Index out of range.\n";
    return false;
  }
  const auto& element = m_elements[i];
  xv = element.xv;
  yv = element.yv;
  zv = element.zv;
  interface = element.interface;
  bc = element.bc;
  lambda = element.lambda;
  return true;
}

GetNumberOfElements()

unsigned int Garfield::ComponentNeBem3d::GetNumberOfElements ( ) const

inline

Definition at line 28 of file ComponentNeBem3d.hh.

{ return m_elements.size(); }

GetNumberOfPrimitives()

unsigned int Garfield::ComponentNeBem3d::GetNumberOfPrimitives ( ) const

inline

Definition at line 23 of file ComponentNeBem3d.hh.

{ return m_primitives.size(); }

GetPrimitive()

bool Garfield::ComponentNeBem3d::GetPrimitive	(	const unsigned int	i,
		double &	a,
		double &	b,
		double &	c,
		std::vector< double > &	xv,
		std::vector< double > &	yv,
		std::vector< double > &	zv,
		int &	interface,
		double &	v,
		double &	q,
		double &	lambda
	)		const

Definition at line 2475 of file ComponentNeBem3d.cc.

                                                          {
  if (i >= m_primitives.size()) {
    std::cerr << m_className << "::GetPrimitive: Index out of range.\n";
    return false;
  }
  const auto& primitive = m_primitives[i];
  a = primitive.a;
  b = primitive.b;
  c = primitive.c;
  xv = primitive.xv;
  yv = primitive.yv;
  zv = primitive.zv;
  interface = primitive.interface;
  v = primitive.v;
  q = primitive.q;
  lambda = primitive.lambda;
  return true;
}

GetVoltageRange()

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

overridevirtual

Calculate the voltage range [V].

Implements Garfield::ComponentBase.

Definition at line 411 of file ComponentNeBem3d.cc.

                                                                 {
  // Voltage and other bc have to come from the solids
  vmin = vmax = 0;
  return true;
}

Initialise()

bool Garfield::ComponentNeBem3d::Initialise

(

)

Retrieve surface panels, remove contacts and cut polygons to rectangles and right-angle triangles.

Definition at line 427 of file ComponentNeBem3d.cc.

                                  {
  // Reset the lists.
  m_primitives.clear();
  m_elements.clear();
 
  if (!m_geometry) {
    std::cerr << m_className << "::Initialise: Geometry not set.\n";
    return false;
  }
  // Be sure we won't have intersections with the bounding box.
  // TODO! Loop over the solids and call PLACYE, PLACHE, PLABXE, PLASPE, ...
  // Loop over the solids.
  std::map<int, Solid*> solids;
  std::map<int, Solid::BoundaryCondition> bc;
  std::map<int, double> volt;
  std::map<int, double> eps;
  std::map<int, double> charge;
  const unsigned int nSolids = m_geometry->GetNumberOfSolids();
  std::vector<Panel> panelsIn;
  for (unsigned int i = 0; i < nSolids; ++i) {
    Medium* medium = nullptr;
    const auto solid = m_geometry->GetSolid(i, medium);
    if (!solid) continue;
    // Get the panels.
    solid->SolidPanels(panelsIn);
    // Get the boundary condition.
    const auto id = solid->GetId();
    solids[id] = solid;
    bc[id] = solid->GetBoundaryConditionType();
    volt[id] = solid->GetBoundaryPotential();
    charge[id] = solid->GetBoundaryChargeDensity();
    if (!medium) {
      eps[id] = 1.;
    } else {
      eps[id] = medium->GetDielectricConstant();
    }
  }
  // Apply cuts.
  // CALL CELSCT('APPLY')
  // Reduce to basic periodic copy.
  // CALL BEMBAS
 
  // Find contact panels and split into primitives.
 
  // *---------------------------------------------------------------------
  // * PLABEM - Prepares panels for BEM applications: removes the contacts
  // *          and cuts polygons to rectangles and right-angle triangles.
  // *---------------------------------------------------------------------
 
  // Establish tolerances.
  const double epsang = 1.e-6;  // BEMEPA
  const double epsxyz = 1.e-6;  // BEMEPD
  // CALL EPSSET('SET',EPSXYZ,EPSXYZ,EPSXYZ)
 
  const unsigned int nIn = panelsIn.size();
  if (m_debug) {
    std::cout << m_className << "::Initialise: Retrieved " 
              << nIn << " panels from the solids.\n";
  }
  // Keep track of which panels have been processed.
  std::vector<bool> mark(nIn, false);
  // Pick up panels which coincide potentially.
  for (unsigned int i = 0; i < nIn; ++i) {
    // Skip panels already done.
    if (mark[i]) continue;
    // Fetch panel parameters.
    const double a1 = panelsIn[i].a;
    const double b1 = panelsIn[i].b;
    const double c1 = panelsIn[i].c;
    const auto& xp1 = panelsIn[i].xv;
    const auto& yp1 = panelsIn[i].yv;
    const auto& zp1 = panelsIn[i].zv;
    const unsigned int np1 = xp1.size();
    // Establish its norm and offset.
    const double d1 = a1 * xp1[0] + b1 * yp1[0] + c1 * zp1[0];
    if (m_debug) {
      std::cout << "  Panel " << i << "\n    Norm vector: "
                << a1 << ", " << b1 << ", " << c1 << ", " << d1 << ".\n";
    }
    // Rotation matrix.
    std::array<std::array<double, 3>, 3> rot;
    if (fabs(c1) <= fabs(a1) && fabs(c1) <= fabs(b1)) {
      // Rotation: removing C
      rot[0][0] = b1 / sqrt(a1 * a1 + b1 * b1);
      rot[0][1] = -a1 / sqrt(a1 * a1 + b1 * b1);
      rot[0][2] = 0.;
    } else if (fabs(b1) <= fabs(a1) && fabs(b1) <= fabs(c1)) {
      // Rotation: removing B
      rot[0][0] = c1 / sqrt(a1 * a1 + c1 * c1);
      rot[0][1] = 0.;
      rot[0][2] = -a1 / sqrt(a1 * a1 + c1 * c1);
    } else {
      // Rotation: removing A
      rot[0][0] = 0.;
      rot[0][1] = c1 / sqrt(b1 * b1 + c1 * c1);
      rot[0][2] = -b1 / sqrt(b1 * b1 + c1 * c1);
    }
    rot[2][0] = a1;
    rot[2][1] = b1;
    rot[2][2] = c1;
    rot[1][0] = rot[2][1] * rot[0][2] - rot[2][2] * rot[0][1];
    rot[1][1] = rot[2][2] * rot[0][0] - rot[2][0] * rot[0][2];
    rot[1][2] = rot[2][0] * rot[0][1] - rot[2][1] * rot[0][0];
    // Rotate to the x, y plane.
    std::vector<double> xp(np1, 0.);
    std::vector<double> yp(np1, 0.);
    std::vector<double> zp(np1, 0.);
    double zm = 0.;
    for (unsigned int k = 0; k < np1; ++k) {
      xp[k] = rot[0][0] * xp1[k] + rot[0][1] * yp1[k] + rot[0][2] * zp1[k];
      yp[k] = rot[1][0] * xp1[k] + rot[1][1] * yp1[k] + rot[1][2] * zp1[k];
      zp[k] = rot[2][0] * xp1[k] + rot[2][1] * yp1[k] + rot[2][2] * zp1[k];
      zm += zp[k];
    }
    zm /= np1;
    // Store it.
    std::vector<Panel> newPanels;
    std::vector<int> vol1;
    std::vector<int> vol2;
    Panel panel1 = panelsIn[i];
    panel1.xv = xp;
    panel1.yv = yp;
    panel1.zv = zp;
    vol1.push_back(panel1.volume);
    vol2.push_back(-1);
    newPanels.push_back(std::move(panel1));
    // Pick up all matching planes.
    for (unsigned int j = i + 1; j < nIn; ++j) {
      if (mark[j]) continue;
      const double a2 = panelsIn[j].a;
      const double b2 = panelsIn[j].b;
      const double c2 = panelsIn[j].c;
      const auto& xp2 = panelsIn[j].xv;
      const auto& yp2 = panelsIn[j].yv;
      const auto& zp2 = panelsIn[j].zv;
      const unsigned int np2 = xp2.size();
      // See whether this matches the first.
      const double d2 = a2 * xp2[0] + b2 * yp2[0] + c2 * zp2[0];
      // Inner product.
      const double dot = a1 * a2 + b1 * b2 + c1 * c2;
      // Offset between the two planes.
      const double offset = d1 - d2 * dot;
      if (fabs(fabs(dot) - 1.) > epsang || fabs(offset) > epsxyz) continue;
      // Found a match.
      mark[j] = true;
      if (m_debug) std::cout << "    Match with panel " << j << ".\n";
      // Rotate this plane too.
      xp.assign(np2, 0.);
      yp.assign(np2, 0.);
      zp.assign(np2, 0.);
      zm = 0.;
      for (unsigned int k = 0; k < np2; ++k) {
        xp[k] = rot[0][0] * xp2[k] + rot[0][1] * yp2[k] + rot[0][2] * zp2[k];
        yp[k] = rot[1][0] * xp2[k] + rot[1][1] * yp2[k] + rot[1][2] * zp2[k];
        zp[k] = rot[2][0] * xp2[k] + rot[2][1] * yp2[k] + rot[2][2] * zp2[k];
        zm += zp[k];
      }
      zm /= np2;
      // Store it.
      Panel panel2 = panelsIn[j];
      panel2.xv = xp;
      panel2.yv = yp;
      panel2.zv = zp;
      vol1.push_back(panel2.volume);
      vol2.push_back(-1);
      newPanels.push_back(std::move(panel2));
    }
    std::vector<bool> obsolete(newPanels.size(), false);
    // Cut them as long as needed till no contacts remain.
    unsigned int jmin = 0;
    bool change = true;
    while (change) {
      change = false;
      const unsigned int n = newPanels.size();
      for (unsigned int j = 0; j < n; ++j) {
        if (obsolete[j] || j < jmin) continue;
        if (vol1[j] >= 0 && vol2[j] >= 0) continue;
        const auto& panelj = newPanels[j];
        for (unsigned int k = j + 1; k < n; ++k) {
          if (obsolete[k]) continue;
          if (vol1[k] >= 0 && vol2[k] >= 0) continue;
          const auto& panelk = newPanels[k];
          if (m_debug) std::cout << "    Cutting " << j << ", " << k << ".\n";
          // Separate contact and non-contact areas.
          std::vector<Panel> panelsOut;
          std::vector<int> itypo;
          EliminateOverlaps(panelj, panelk, panelsOut, itypo);
          const unsigned int nOut = panelsOut.size();
          if (nOut == 2) {
            // TODO: retrieve epsx, epsy from overlap finding?
            const double epsx = epsxyz;
            const double epsy = epsxyz;
            // If there are just 2 panels, see whether there is a new one.
            const bool equal1 = Equal(panelj, panelsOut[0], epsx, epsy);
            const bool equal2 = Equal(panelj, panelsOut[1], epsx, epsy);
            const bool equal3 = Equal(panelk, panelsOut[0], epsx, epsy);
            const bool equal4 = Equal(panelk, panelsOut[1], epsx, epsy);
            if ((equal1 || equal3) && (equal2 || equal4)) {
              if (m_debug) {
                std::cout << "    Original and new panels are identical.\n";
              }
            } else {
              change = true;
            }
          } else {
            change = true;
          }
          if (m_debug) std::cout << "    Change flag: " << change << ".\n";
          // If there is no change, keep the two panels and proceed.
          if (!change) continue;
          // Flag the existing panels as inactive.
          obsolete[j] = true;
          obsolete[k] = true;
 
          // Add the new panels.
          for (unsigned int l = 0; l < nOut; ++l) {
            if (itypo[l] == 1) {
              vol1.push_back(std::max(vol1[j], vol2[j]));
              vol2.push_back(-1);
            } else if (itypo[l] == 2) {
              vol1.push_back(std::max(vol1[k], vol2[k]));
              vol2.push_back(-1);
            } else {
              vol1.push_back(std::max(vol1[j], vol2[j]));
              vol2.push_back(std::max(vol1[k], vol2[k]));
            }
            newPanels.push_back(std::move(panelsOut[l]));
            obsolete.push_back(false);
          }
          jmin = j + 1;
          // Restart the loops.
          break;
        }
        if (change) break;
      }
    }
    // And rotate the panels back in place.
    const unsigned int nNew = newPanels.size();
    for (unsigned int j = 0; j < nNew; ++j) {
      if (obsolete[j]) continue;
      // Examine the boundary conditions.
      int interfaceType = 0;
      double potential = 0.;
      double lambda = 0.;
      double chargeDensity = 0.;
      if (m_debug) {
        std::cout << "    Volume 1: " << vol1[j] 
                  << ". Volume 2: " << vol2[j] << ".\n";
      }
      if (vol1[j] < 0 && vol2[j] < 0) {
        // Shouldn't happen.
        continue;
      } else if (vol1[j] < 0 || vol2[j] < 0) {
        // Interface between a solid and vacuum/background.
        const auto vol = vol1[j] < 0 ? vol2[j] : vol1[j];
        interfaceType = InterfaceType(bc[vol]);
        if (bc[vol] == Solid::Dielectric || 
            bc[vol] == Solid::DielectricCharge) {
          if (fabs(eps[vol] - 1.) < 1.e-6) {
            // Same epsilon on both sides. Skip.
            interfaceType = 0;
          } else {
            lambda = (eps[vol] - 1.) / (eps[vol] + 1.);
          }
        } else if (bc[vol] == Solid::Voltage) {
          potential = volt[vol];
        }
        if (bc[vol] == Solid::Charge || bc[vol] == Solid::DielectricCharge) {
          chargeDensity = charge[vol];
        }
      } else {
        const auto bc1 = bc[vol1[j]];
        const auto bc2 = bc[vol2[j]];
        if (bc1 == Solid::Voltage || bc1 == Solid::Charge ||
            bc1 == Solid::Float) {
          interfaceType = InterfaceType(bc1);
          // First volume is a conductor. Other volume must be a dielectric.
          if (bc2 == Solid::Dielectric || bc2 == Solid::DielectricCharge) {
            if (bc1 == Solid::Voltage) {
              potential = volt[vol1[j]];
            } else if (bc1 == Solid::Charge) {
              chargeDensity = charge[vol1[j]];
            }
          } else {
            interfaceType = -1;
          }
          if (bc1 == Solid::Voltage && bc2 == Solid::Voltage) {
            const double v1 = volt[vol1[j]];
            const double v2 = volt[vol2[j]];
            if (fabs(v1 - v2) < 1.e-6 * (1. + fabs(v1) + fabs(v2))) {
              interfaceType = 0;
            }
          }
        } else if (bc1 == Solid::Dielectric || 
                   bc1 == Solid::DielectricCharge) {
          interfaceType = InterfaceType(bc2);
          // First volume is a dielectric.
          if (bc2 == Solid::Voltage) {
            potential = volt[vol2[j]];
          } else if (bc2 == Solid::Charge) {
            chargeDensity = charge[vol2[j]];
          } else if (bc2 == Solid::Dielectric || 
                     bc2 == Solid::DielectricCharge) {
            const double eps1 = eps[vol1[j]];
            const double eps2 = eps[vol2[j]];
            if (fabs(eps1 - eps2) < 1.e-6 * (1. + fabs(eps1) + fabs(eps2))) {
              // Same epsilon. Skip.
              interfaceType = 0;
            } else {
              lambda = (eps1 - eps2) / (eps1 + eps2);
            }
          }
        }
      }
      if (m_debug) {
        if (interfaceType < 0) {
          std::cout << "    Conflicting boundary conditions. Skip.\n";
        } else if (interfaceType < 1) {
          std::cout << "    Trivial interface. Skip.\n";
        } else if (interfaceType > 5) {
          std::cout << "    Interface type " << interfaceType 
                    << " is not implemented. Skip.\n";
        }
      }
      if (interfaceType < 1 || interfaceType > 5) continue;
 
      std::vector<Panel> panelsOut;
      // Reduce to rectangles and right-angle triangles.
      if (m_debug) std::cout << "    Creating primitives.\n";
      MakePrimitives(newPanels[j], panelsOut);
      // Loop over the rectangles and triangles.
      for (auto& panel : panelsOut) {
        const auto& up = panel.xv;
        const auto& vp = panel.yv;
        const auto& wp = panel.zv;
        const unsigned int np = up.size();
        // Rotate.
        xp.assign(np, 0.);
        yp.assign(np, 0.);
        zp.assign(np, 0.);
        for (unsigned int k = 0; k < np; ++k) {
          xp[k] = rot[0][0] * up[k] + rot[1][0] * vp[k] + rot[2][0] * wp[k];
          yp[k] = rot[0][1] * up[k] + rot[1][1] * vp[k] + rot[2][1] * wp[k];
          zp[k] = rot[0][2] * up[k] + rot[1][2] * vp[k] + rot[2][2] * wp[k];
        }
        Primitive primitive;
        primitive.a = panel.a;
        primitive.b = panel.b;
        primitive.c = panel.c;
        primitive.xv = xp;
        primitive.yv = yp;
        primitive.zv = zp;
        primitive.v = potential;
        primitive.q = chargeDensity;
        primitive.lambda = lambda;
        primitive.interface = interfaceType;
        // Set the requested discretization level (target element size).
        primitive.elementSize = -1.;
        if (solids.find(vol1[j]) != solids.end()) {
          const auto solid = solids[vol1[j]];
          if (solid) {
            primitive.elementSize = solid->GetDiscretisationLevel(panel);
          }
        }
        m_primitives.push_back(std::move(primitive));
      }
    }
  }
 
  if (m_debug) {
    std::cout << m_className << "::Initialise:\n"
              << "    Created " << m_primitives.size() << " primitives.\n";
  }
 
  // Discretize the primitives.
  for (const auto& primitive : m_primitives) {
    const auto nVertices = primitive.xv.size();
    if (nVertices < 2 || nVertices > 4) continue;
    std::vector<Element> elements;
    // Get the target element size.
    double targetSize = primitive.elementSize;
    if (targetSize < MinDist) targetSize = m_targetElementSize;
    if (nVertices == 2) {
      DiscretizeWire(primitive, targetSize, elements);
    } else if (nVertices == 3) {
      DiscretizeTriangle(primitive, targetSize, elements);
    } else if (nVertices == 4) {
      DiscretizeRectangle(primitive, targetSize, elements);
    }
    for (auto& element: elements) {
      element.interface = primitive.interface;
      element.lambda = primitive.lambda;
      element.bc = primitive.v;
      element.assigned = primitive.q;
      element.solution = 0.;
    }
    m_elements.insert(m_elements.end(),
                      std::make_move_iterator(elements.begin()),
                      std::make_move_iterator(elements.end()));
  }
  return true;
}

Referenced by ElectricField().

Reset()

void Garfield::ComponentNeBem3d::Reset ( )

overridevirtual

Reset the component.

Implements Garfield::ComponentBase.

Definition at line 2519 of file ComponentNeBem3d.cc.

                             {
  m_primitives.clear();
  m_elements.clear();
  m_ready = false;
}

SetTargetElementSize()

void Garfield::ComponentNeBem3d::SetTargetElementSize

(

const double

length

)

Set the default value of the target linear size of the elements produced by neBEM's discretisation process.

Definition at line 417 of file ComponentNeBem3d.cc.

                                                               {
 
  if (length < MinDist) {
    std::cerr << m_className << "::SetTargetElementSize: Value must be > "
              << MinDist << ".\n";
    return; 
  }
  m_targetElementSize = length;
}

UpdatePeriodicity()

void Garfield::ComponentNeBem3d::UpdatePeriodicity ( )

overridevirtual

Verify periodicities.

Implements Garfield::ComponentBase.

Definition at line 2525 of file ComponentNeBem3d.cc.

                                         {
  std::cerr << m_className << "::UpdatePeriodicity:\n"
            << "    Periodicities are not supported.\n";
}

---

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

• ComponentNeBem3d.hh
• ComponentNeBem3d.cc