ComponentTcad2d.cc

Go to the documentation of this file.

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
 
#include "Garfield/ComponentTcad2d.hh"
#include "Garfield/GarfieldConstants.hh"
#include "Garfield/Utilities.hh"
 
namespace Garfield {
 
ComponentTcad2d::ComponentTcad2d() : ComponentBase() {
  m_className = "ComponentTcad2d";
 
  m_regions.reserve(10);
  m_vertices.reserve(10000);
  m_elements.reserve(10000);
}
 
bool ComponentTcad2d::SetDonor(const unsigned int donorNumber,
                               const double eXsec, const double hXsec,
                               const double conc) {
  if (donorNumber >= m_donors.size()) {
    std::cerr << m_className << "::SetDonor: Index out of range.\n";
    return false;
  }
  m_donors[donorNumber].xsece = eXsec;
  m_donors[donorNumber].xsech = hXsec;
  m_donors[donorNumber].conc = conc;
 
  m_validTraps = CheckTraps();
  return true;
}
 
bool ComponentTcad2d::SetAcceptor(const unsigned int acceptorNumber,
                                  const double eXsec, const double hXsec,
                                  const double conc) {
  if (acceptorNumber >= m_acceptors.size()) {
    std::cerr << m_className << "::SetAcceptor: Index out of range.\n";
    return false;
  }
  m_acceptors[acceptorNumber].xsece = eXsec;
  m_acceptors[acceptorNumber].xsech = hXsec;
  m_acceptors[acceptorNumber].conc = conc;
 
  m_validTraps = CheckTraps();
  return true;
}
 
bool ComponentTcad2d::ElectronAttachment(const double x, const double y,
                                         const double z, double& eta) {
  eta = 0.;
  if (!m_validTraps) {
    std::cerr << m_className << "::ElectronAttachment:\n"
              << "    Trap cross-sections or concentrations are invalid.\n";
    return false;
  }
 
  if (m_donors.empty() && m_acceptors.empty()) {
    if (m_debug) {
      std::cerr << m_className << "::ElectronAttachment: No traps defined.\n";
    }
    return true;
  }
 
  const unsigned int nAcceptors = m_acceptors.size();
  for (unsigned int i = 0; i < nAcceptors; ++i) {
    const Defect& acceptor = m_acceptors[i];
    double f = 0.;
    GetAcceptorOccupation(x, y, z, i, f);
    eta += acceptor.conc * acceptor.xsece * (1. - f);
  }
  const unsigned int nDonors = m_donors.size();
  for (unsigned int i = 0; i < nDonors; ++i) {
    const Defect& donor = m_donors[i];
    double f = 0.;
    GetDonorOccupation(x, y, z, i, f);
    eta += donor.conc * donor.xsece * f;
  }
  return true;
}
 
bool ComponentTcad2d::HoleAttachment(const double x, const double y,
                                     const double z, double& eta) {
  eta = 0.;
  if (!m_validTraps) {
    std::cerr << m_className << "::HoleAttachment:\n"
              << "    Trap cross-sections or concentrations are invalid.\n";
    return false;
  }
 
  if (m_donors.empty() && m_acceptors.empty()) {
    if (m_debug) {
      std::cerr << m_className << "::HoleAttachment: No traps defined.\n";
    }
    return true;
  }
 
  const unsigned int nAcceptors = m_acceptors.size();
  for (unsigned int i = 0; i < nAcceptors; ++i) {
    const Defect& acceptor = m_acceptors[i];
    double f = 0.;
    GetAcceptorOccupation(x, y, z, i, f);
    eta += acceptor.conc * acceptor.xsech * f;
  }
  const unsigned int nDonors = m_donors.size();
  for (unsigned int i = 0; i < nDonors; ++i) {
    const Defect& donor = m_donors[i];
    double f = 0.;
    GetDonorOccupation(x, y, z, i, f);
    eta += donor.conc * donor.xsech * (1. - f);
  }
  return true;
}
 
void ComponentTcad2d::WeightingField(const double xin, const double yin,
                                     const double zin, double& wx, double& wy,
                                     double& wz, const std::string& /*label*/) {
  wx = wy = wz = 0.;
  if (m_wf.empty()) {
    std::cerr << m_className << "::WeightingField: Not available.\n";
    return;
  }
  // In case of periodicity, reduce to the cell volume.
  double x = xin, y = yin, z = zin;
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) return;
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const auto& f = m_wf[element.vertex[j]];
    wx += w[j] * f[0];
    wy += w[j] * f[1];
  }
  if (xmirr) wx = -wx;
  if (ymirr) wy = -wy;
}
 
double ComponentTcad2d::WeightingPotential(const double xin, const double yin,
                                           const double zin, 
                                           const std::string& /*label*/) {
 
  if (m_wp.empty()) {
    std::cerr << m_className << "::WeightingPotential: Not available.\n";
    return 0.;
  }
  // In case of periodicity, reduce to the cell volume.
  double x = xin, y = yin, z = zin;
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return 0.;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) return 0.;
 
  double v = 0.;
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    v += w[j] * m_wp[element.vertex[j]];
  }
  return v;
}
 
void ComponentTcad2d::ElectricField(const double xin, const double yin,
                                    const double zin, double& ex, double& ey,
                                    double& ez, double& p, Medium*& m,
                                    int& status) {
  // Assume this will work.
  status = 0;
  // Initialise.
  ex = ey = ez = p = 0.;
  m = nullptr;
 
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::ElectricField:\n"
              << "    Field map is not available for interpolation.\n";
    status = -10;
    return;
  }
 
  // In case of periodicity, reduce to the cell volume.
  double x = xin, y = yin, z = zin;
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) {
    status = -6;
    return;
  }
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    status = -6;
    return;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    ex += w[j] * vj.ex;
    ey += w[j] * vj.ey;
    p += w[j] * vj.p;
  }
  if (xmirr) ex = -ex;
  if (ymirr) ey = -ey;
  m = m_regions[element.region].medium;
  if (!m_regions[element.region].drift || !m) status = -5;
  m_lastElement = i;
}
 
void ComponentTcad2d::ElectronVelocity(const double xin, const double yin,
                                       const double zin, double& vx, double& vy,
                                       double& vz, Medium*& m, int& status) {
  // Assume this will work.
  status = 0;
  // Initialise.
  vx = vy = vz = 0.;
  m = nullptr;
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::ElectronVelocity:\n"
              << "    Field map is not available for interpolation.\n";
    status = -10;
    return;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  if (!InsideBoundingBox(x, y, z)) {
    status = -6;
    return;
  }
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    status = -6;
    return;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    vx += w[j] * vj.eVx;
    vy += w[j] * vj.eVy;
  }
  if (xmirr) vx = -vx;
  if (ymirr) vy = -vy;
  m = m_regions[element.region].medium;
  if (!m_regions[element.region].drift || !m) status = -5;
  m_lastElement = i;
}
 
void ComponentTcad2d::HoleVelocity(const double xin, const double yin,
                                   const double zin, double& vx, double& vy,
                                   double& vz, Medium*& m, int& status) {
 
  // Assume this will work.
  status = 0;
  // Initialise.
  vx = vy = vz = 0.;
  m = nullptr;
 
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::HoleVelocity:\n"
              << "    Field map is not available for interpolation.\n";
    status = -10;
    return;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  if (!InsideBoundingBox(x, y, z)) {
    status = -6;
    return;
  }
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    status = -6;
    return;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    vx += w[j] * vj.hVx;
    vy += w[j] * vj.hVy;
  }
  if (xmirr) vx = -vx;
  if (ymirr) vy = -vy;
  m = m_regions[element.region].medium;
  if (!m_regions[element.region].drift || !m) status = -5;
  m_lastElement = i;
}
 
Medium* ComponentTcad2d::GetMedium(const double xin, const double yin,
                                   const double zin) {
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::GetMedium:\n"
              << "    Field map not available for interpolation.\n";
    return nullptr;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return nullptr;
 
  // Shape functions
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    return nullptr;
  }
  m_lastElement = i;
  const Element& element = m_elements[i];
  return m_regions[element.region].medium;
}
 
bool ComponentTcad2d::GetElectronLifetime(const double xin, const double yin,
                                          const double zin, double& tau) {
  tau = 0.;
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::GetElectronLifetime:\n"
              << "    Field map is not available for interpolation.\n";
    return false;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return false;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    return false;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    tau += w[j] * vj.eTau;
  }
  m_lastElement = i;
  return true;
}
bool ComponentTcad2d::GetHoleLifetime(const double xin, const double yin,
                                      const double zin, double& tau) {
  tau = 0.;
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::GetHoleLifetime:\n"
              << "    Field map is not available for interpolation.\n";
    return false;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return false;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    return false;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    tau += w[j] * vj.hTau;
  }
  m_lastElement = i;
  return true;
}
 
bool ComponentTcad2d::GetMobility(const double xin, const double yin,
                                  const double zin, double& emob,
                                  double& hmob) {
  emob = hmob = 0.;
 
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::GetMobility:\n"
              << "    Field map is not available for interpolation.\n";
    return false;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return false;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    return false;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    emob += w[j] * vj.emob;
    hmob += w[j] * vj.hmob;
  }
  m_lastElement = i;
  return true;
}
 
bool ComponentTcad2d::GetDonorOccupation(const double xin, const double yin,
                                         const double zin,
                                         const unsigned int donorNumber,
                                         double& f) {
  f = 0.;
  if (donorNumber >= m_donors.size()) {
    std::cerr << m_className << "::GetDonorOccupation:\n"
              << "    Donor " << donorNumber << " does not exist.\n";
    return false;
  }
 
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::GetDonorOccupation:\n"
              << "    Field map is not available for interpolation.\n";
    return false;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return false;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    return false;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    f += w[j] * vj.donorOcc[donorNumber];
  }
  m_lastElement = i;
  return true;
}
 
bool ComponentTcad2d::GetAcceptorOccupation(const double xin, const double yin,
                                            const double zin,
                                            const unsigned int acceptorNumber,
                                            double& f) {
  f = 0.;
  if (acceptorNumber >= m_acceptors.size()) {
    std::cerr << m_className << "::GetAcceptorOccupation:\n"
              << "    Acceptor " << acceptorNumber << " does not exist.\n";
    return false;
  }
 
  // Make sure the field map has been loaded.
  if (!m_ready) {
    std::cerr << m_className << "::GetAcceptorOccupation:\n"
              << "    Field map is not available for interpolation.\n";
    return false;
  }
 
  double x = xin, y = yin, z = zin;
  // In case of periodicity, reduce to the cell volume.
  bool xmirr = false, ymirr = false;
  MapCoordinates(x, y, xmirr, ymirr);
  // Check if the point is inside the bounding box.
  if (!InsideBoundingBox(x, y, z)) return false;
 
  std::array<double, nMaxVertices> w;
  const unsigned int i = FindElement(x, y, w);
  if (i >= m_elements.size()) {
    // Point is outside the mesh.
    return false;
  }
 
  const Element& element = m_elements[i];
  const unsigned int nVertices = element.type + 1;
  for (unsigned int j = 0; j < nVertices; ++j) {
    const Vertex& vj = m_vertices[element.vertex[j]];
    f += w[j] * vj.acceptorOcc[acceptorNumber];
  }
  m_lastElement = i;
  return true;
}
 
bool ComponentTcad2d::Initialise(const std::string& gridfilename,
                                 const std::string& datafilename) {
  m_ready = false;
 
  m_hasPotential = m_hasField = false;
  m_hasElectronMobility = m_hasHoleMobility = false;
  m_validTraps = false;
  m_donors.clear();
  m_acceptors.clear();
 
  // Import mesh data from .grd file.
  if (!LoadGrid(gridfilename)) {
    std::cerr << m_className << "::Initialise:\n"
              << "    Importing mesh data failed.\n";
    return false;
  }
 
  // Import electric field, potential, mobilities and
  // trap occupation values from .dat file.
  if (!LoadData(datafilename)) {
    std::cerr << m_className << "::Initialise:\n"
              << "    Importing electric field and potential failed.\n";
    return false;
  }
 
  // Find min./max. coordinates and potentials.
  m_xMaxBB = m_xMinBB = m_vertices[m_elements[0].vertex[0]].x;
  m_yMaxBB = m_yMinBB = m_vertices[m_elements[0].vertex[0]].y;
  m_pMax = m_pMin = m_vertices[m_elements[0].vertex[0]].p;
  for (auto& element : m_elements) {
    const Vertex& v0 = m_vertices[element.vertex[0]];
    const Vertex& v1 = m_vertices[element.vertex[1]];
    double xmin = std::min(v0.x, v1.x);
    double xmax = std::max(v0.x, v1.x);
    double ymin = std::min(v0.y, v1.y);
    double ymax = std::max(v0.y, v1.y);
    m_pMin = std::min(m_pMin, std::min(v0.p, v1.p));
    m_pMax = std::max(m_pMax, std::max(v0.p, v1.p));
    if (element.type > 1) {
      const Vertex& v2 = m_vertices[element.vertex[2]];
      xmin = std::min(xmin, v2.x);
      xmax = std::max(xmax, v2.x);
      ymin = std::min(ymin, v2.y);
      ymax = std::max(ymax, v2.y);
      m_pMin = std::min(m_pMin, v2.p);
      m_pMax = std::max(m_pMax, v2.p);
    }
    if (element.type > 2) {
      const Vertex& v3 = m_vertices[element.vertex[3]];
      xmin = std::min(xmin, v3.x);
      xmax = std::max(xmax, v3.x);
      ymin = std::min(ymin, v3.y);
      ymax = std::max(ymax, v3.y);
      m_pMin = std::min(m_pMin, v3.p);
      m_pMax = std::max(m_pMax, v3.p);
    }
    constexpr double tol = 1.e-6;
    element.xmin = xmin - tol;
    element.xmax = xmax + tol;
    element.ymin = ymin - tol;
    element.ymax = ymax + tol;
    m_xMinBB = std::min(m_xMinBB, xmin);
    m_xMaxBB = std::max(m_xMaxBB, xmax);
    m_yMinBB = std::min(m_yMinBB, ymin);
    m_yMaxBB = std::max(m_yMaxBB, ymax);
  }
 
  std::cout << m_className << "::Initialise:\n"
            << "    Available data:\n";
  if (m_hasPotential) {
    std::cout << "      Electrostatic potential\n";
  }
  if (m_hasField) {
    std::cout << "      Electric field\n";
  }
  if (m_hasElectronMobility) {
    std::cout << "      Electron mobility\n";
  }
  if (m_hasHoleMobility) {
    std::cout << "      Hole mobility\n";
  }
  if (m_hasElectronVelocity) {
    std::cout << "      Electron velocity\n";
  }
  if (m_hasHoleVelocity) {
    std::cout << "      Hole velocity\n";
  }
  if (m_hasElectronLifetime) {
    std::cout << "      Electron lifetimes\n";
  }
  if (m_hasHoleLifetime) {
    std::cout << "      Hole lifetimes\n";
  }
  if (!m_donors.empty()) {
    std::cout << "      " << m_donors.size() << " donor-type traps\n";
  }
  if (!m_acceptors.empty()) {
    std::cout << "      " << m_acceptors.size() << " acceptor-type traps\n";
  }
  std::cout << "    Bounding box:\n"
            << "      " << m_xMinBB << " < x [cm] < " << m_xMaxBB << "\n"
            << "      " << m_yMinBB << " < y [cm] < " << m_yMaxBB << "\n"
            << "    Voltage range:\n"
            << "      " << m_pMin << " < V < " << m_pMax << "\n";
 
  bool ok = true;
 
  // Count the number of elements belonging to a region.
  const int nRegions = m_regions.size();
  std::vector<int> nElementsRegion(nRegions, 0);
 
  // Count the different element shapes.
  unsigned int nLines = 0;
  unsigned int nTriangles = 0;
  unsigned int nRectangles = 0;
  unsigned int nOtherShapes = 0;
 
  // Check if there are elements which are not part of any region.
  unsigned int nLoose = 0;
  std::vector<int> looseElements;
 
  // Check if there are degenerate elements.
  unsigned int nDegenerate = 0;
  std::vector<int> degenerateElements;
 
  const unsigned int nElements = m_elements.size();
  for (unsigned int i = 0; i < nElements; ++i) {
    const Element& element = m_elements[i];
    if (element.type == 1) {
      ++nLines;
      if (element.vertex[0] == element.vertex[1]) {
        degenerateElements.push_back(i);
        ++nDegenerate;
      }
    } else if (element.type == 2) {
      ++nTriangles;
      if (element.vertex[0] == element.vertex[1] ||
          element.vertex[1] == element.vertex[2] ||
          element.vertex[2] == element.vertex[0]) {
        degenerateElements.push_back(i);
        ++nDegenerate;
      }
    } else if (element.type == 3) {
      ++nRectangles;
      if (element.vertex[0] == element.vertex[1] ||
          element.vertex[0] == element.vertex[2] ||
          element.vertex[0] == element.vertex[3] ||
          element.vertex[1] == element.vertex[2] ||
          element.vertex[1] == element.vertex[3] ||
          element.vertex[2] == element.vertex[3]) {
        degenerateElements.push_back(i);
        ++nDegenerate;
      }
    } else {
      // Other shapes should not occur, since they were excluded in LoadGrid.
      ++nOtherShapes;
    }
    if (element.region >= 0 && element.region < nRegions) {
      ++nElementsRegion[element.region];
    } else {
      looseElements.push_back(i);
      ++nLoose;
    }
  }
 
  if (nDegenerate > 0) {
    std::cerr << m_className << "::Initialise:\n"
              << "    The following elements are degenerate:\n";
    for (unsigned int i = 0; i < nDegenerate; ++i) {
      std::cerr << "      " << degenerateElements[i] << "\n";
    }
    ok = false;
  }
 
  if (nLoose > 0) {
    std::cerr << m_className << "::Initialise:\n"
              << "    The following elements are not part of any region:\n";
    for (unsigned int i = 0; i < nLoose; ++i) {
      std::cerr << "      " << looseElements[i] << "\n";
    }
    ok = false;
  }
 
  std::cout << m_className << "::Initialise:\n"
            << "    Number of regions: " << nRegions << "\n";
  for (int i = 0; i < nRegions; ++i) {
    std::cout << "      " << i << ": " << m_regions[i].name << ", "
              << nElementsRegion[i] << " elements\n";
  }
 
  std::cout << "    Number of elements: " << nElements << "\n";
  if (nLines > 0) {
    std::cout << "      " << nLines << " lines\n";
  }
  if (nTriangles > 0) {
    std::cout << "      " << nTriangles << " triangles\n";
  }
  if (nRectangles > 0) {
    std::cout << "      " << nRectangles << " rectangles\n";
  }
  if (nOtherShapes > 0) {
    std::cerr << "      " << nOtherShapes << " elements of unknown type.\n"
              << "      Program bug!\n";
    m_ready = false;
    Cleanup();
    return false;
  }
 
  std::cout << "    Number of vertices: " << m_vertices.size() << "\n";
 
  // Find adjacent elements.
  std::cout << m_className << "::Initialise:\n"
            << "    Looking for neighbouring elements. Be patient...\n";
  FindNeighbours();
 
  if (!ok) {
    m_ready = false;
    Cleanup();
    return false;
  }
 
  m_ready = true;
  UpdatePeriodicity();
  std::cout << m_className << "::Initialise: Initialisation finished.\n";
  return true;
}
 
bool ComponentTcad2d::SetWeightingField(const std::string& datfile1,
                                        const std::string& datfile2,
                                        const double dv) {
 
  if (!m_ready) {
    std::cerr << m_className << "::SetWeightingField:\n"
              << "    Mesh is not available. Call Initialise first.\n";
    return false;
  }
  if (dv < Small) {
     std::cerr << m_className << "::SetWeightingField:\n"
               << "    Voltage difference must be > 0.\n";
     return false;
  }
  const double s = 1. / dv;
 
  m_wf.clear();
  m_wp.clear();
  // Load first the field/potential at nominal bias.
  std::vector<std::array<double, 2> > wf1;
  std::vector<double> wp1;
  if (!LoadWeightingField(datfile1, wf1, wp1)) {
    std::cerr << m_className << "::SetWeightingField:\n"
              << "    Could not import data from " << datfile1 << ".\n";
    return false;
  }
  // Then load the field/potential for the configuration with the potential 
  // at the electrode to be read out increased by small voltage dv. 
  std::vector<std::array<double, 2> > wf2;
  std::vector<double> wy2;
  std::vector<double> wp2;
  if (!LoadWeightingField(datfile2, wf2, wp2)) {
    std::cerr << m_className << "::SetWeightingField:\n"
              << "    Could not import data from " << datfile2 << ".\n";
    return false;
  }
  const unsigned int nVertices = m_vertices.size();
  bool foundField = true;
  if (wf1.size() != nVertices || wf2.size() != nVertices) {
    foundField = false;
    std::cerr << m_className << "::SetWeightingField:\n"
              << "    Could not load electric field values.\n";
  }
  bool foundPotential = true;
  if (wp1.size() != nVertices || wp2.size() != nVertices) {
    foundPotential = false;
    std::cerr << m_className << "::SetWeightingField:\n"
              << "    Could not load electrostatic potentials.\n";
  }
  if (!foundField && !foundPotential) return false;
  if (foundField) {
    m_wf.assign(nVertices, {0., 0.});
    for (unsigned int i = 0; i < nVertices; ++i) {
      m_wf[i][0] = (wf2[i][0] - wf1[i][0]) * s; 
      m_wf[i][1] = (wf2[i][1] - wf1[i][1]) * s; 
    }
  }
  if (foundPotential) {
    m_wp.assign(nVertices, 0.);
    for (unsigned int i = 0; i < nVertices; ++i) {
      m_wp[i] = (wp2[i] - wp1[i]) * s; 
    }
  }
  return true;
}
 
bool ComponentTcad2d::GetBoundingBox(double& xmin, double& ymin, double& zmin,
                                     double& xmax, double& ymax, double& zmax) {
  if (!m_ready) return false;
  if (m_periodic[0] || m_mirrorPeriodic[0]) {
    xmin = -INFINITY;
    xmax = +INFINITY;
  } else {
    xmin = m_xMinBB;
    xmax = m_xMaxBB;
  }
 
  if (m_periodic[1] || m_mirrorPeriodic[1]) {
    ymin = -INFINITY;
    ymax = +INFINITY;
  } else {
    ymin = m_yMinBB;
    ymax = m_yMaxBB;
  }
 
  if (m_hasRangeZ) {
    zmin = m_zMinBB;
    zmax = m_zMaxBB;
  }
  return true;
}
 
void ComponentTcad2d::SetRangeZ(const double zmin, const double zmax) {
  if (fabs(zmax - zmin) <= 0.) {
    std::cerr << m_className << "::SetRangeZ: Zero range is not permitted.\n";
    return;
  }
  m_zMinBB = std::min(zmin, zmax);
  m_zMaxBB = std::max(zmin, zmax);
  m_hasRangeZ = true;
}
 
bool ComponentTcad2d::GetVoltageRange(double& vmin, double& vmax) {
  if (!m_ready) return false;
  vmin = m_pMin;
  vmax = m_pMax;
  return true;
}
 
void ComponentTcad2d::PrintRegions() const {
  // Do not proceed if not properly initialised.
  if (!m_ready) {
    std::cerr << m_className << "::PrintRegions:\n"
              << "    Field map not yet initialised.\n";
    return;
  }
 
  if (m_regions.empty()) {
    std::cerr << m_className << "::PrintRegions:\n"
              << "    No regions are currently defined.\n";
    return;
  }
 
  const unsigned int nRegions = m_regions.size();
  std::cout << m_className << "::PrintRegions:\n"
            << "    Currently " << nRegions << " regions are defined.\n"
            << "      Index  Name       Medium\n";
  for (unsigned int i = 0; i < nRegions; ++i) {
    std::cout << "      " << i << "  " << m_regions[i].name;
    if (!m_regions[i].medium) {
      std::cout << "      none  ";
    } else {
      std::cout << "      " << m_regions[i].medium->GetName();
    }
    if (m_regions[i].drift) {
      std::cout << " (active region)\n";
    } else {
      std::cout << "\n";
    }
  }
}
 
void ComponentTcad2d::GetRegion(const unsigned int i, std::string& name,
                                bool& active) const {
  if (i >= m_regions.size()) {
    std::cerr << m_className << "::GetRegion: Index out of range.\n";
    return;
  }
  name = m_regions[i].name;
  active = m_regions[i].drift;
}
 
void ComponentTcad2d::SetDriftRegion(const unsigned int i) {
  if (i >= m_regions.size()) {
    std::cerr << m_className << "::SetDriftRegion: Index out of range.\n";
    return;
  }
  m_regions[i].drift = true;
}
 
void ComponentTcad2d::UnsetDriftRegion(const unsigned int i) {
  if (i >= m_regions.size()) {
    std::cerr << m_className << "::UnsetDriftRegion: Index out of range.\n";
    return;
  }
  m_regions[i].drift = false;
}
 
void ComponentTcad2d::SetMedium(const unsigned int i, Medium* medium) {
  if (i >= m_regions.size()) {
    std::cerr << m_className << "::SetMedium: Index out of range.\n";
    return;
  }
 
  if (!medium) {
    std::cerr << m_className << "::SetMedium: Null pointer.\n";
    return;
  }
 
  m_regions[i].medium = medium;
}
 
Medium* ComponentTcad2d::GetMedium(const unsigned int i) const {
  if (i >= m_regions.size()) {
    std::cerr << m_className << "::GetMedium: Index out of range.\n";
    return nullptr;
  }
 
  return m_regions[i].medium;
}
 
bool ComponentTcad2d::GetElement(const unsigned int i, double& vol,
                                 double& dmin, double& dmax, int& type) const {
  if (i >= m_elements.size()) {
    std::cerr << m_className << "::GetElement: Index out of range.\n";
    return false;
  }
 
  const Element& element = m_elements[i];
  if (element.type == 1) {
    const Vertex& v0 = m_vertices[element.vertex[0]];
    const Vertex& v1 = m_vertices[element.vertex[1]];
    const double dx = v1.x - v0.x;
    const double dy = v1.y - v0.y;
    const double d = sqrt(dx * dx + dy * dy);
    dmin = dmax = vol = d;
  } else if (m_elements[i].type == 2) {
    const Vertex& v0 = m_vertices[element.vertex[0]];
    const Vertex& v1 = m_vertices[element.vertex[1]];
    const Vertex& v2 = m_vertices[element.vertex[2]];
    vol = 0.5 *
          fabs((v2.x - v0.x) * (v1.y - v0.y) - (v2.y - v0.y) * (v1.x - v0.x));
    const double a = sqrt(pow(v1.x - v0.x, 2) + pow(v1.y - v0.y, 2));
    const double b = sqrt(pow(v2.x - v0.x, 2) + pow(v2.y - v0.y, 2));
    const double c = sqrt(pow(v1.x - v2.x, 2) + pow(v1.y - v2.y, 2));
    dmin = std::min(std::min(a, b), c);
    dmax = std::max(std::max(a, b), c);
  } else if (m_elements[i].type == 3) {
    const Vertex& v0 = m_vertices[element.vertex[0]];
    const Vertex& v1 = m_vertices[element.vertex[1]];
    const Vertex& v3 = m_vertices[element.vertex[3]];
    const double a = sqrt(pow(v1.x - v0.x, 2) + pow(v1.y - v0.y, 2));
    const double b = sqrt(pow(v3.x - v0.x, 2) + pow(v3.y - v0.y, 2));
    vol = a * b;
    dmin = std::min(a, b);
    dmax = sqrt(a * a + b * b);
  } else {
    std::cerr << m_className << "::GetElement:\n"
              << "    Unexpected element type (" << type << ")\n";
    return false;
  }
  return true;
}
 
bool ComponentTcad2d::GetElement(const unsigned int i, double& vol,
                                 double& dmin, double& dmax, int& type,
                                 int& node1, int& node2, int& node3, int& node4,
                                 int& reg) const {
  if (!GetElement(i, vol, dmin, dmax, type)) return false;
  const Element& element = m_elements[i];
  node1 = element.vertex[0];
  node2 = element.vertex[1];
  node3 = element.vertex[2];
  node4 = element.vertex[3];
  reg = element.region;
  return true;
}
 
bool ComponentTcad2d::GetNode(const unsigned int i, double& x, double& y,
                              double& v, double& ex, double& ey) const {
  if (i >= m_vertices.size()) {
    std::cerr << m_className << "::GetNode: Index out of range.\n";
    return false;
  }
 
  x = m_vertices[i].x;
  y = m_vertices[i].y;
  v = m_vertices[i].p;
  ex = m_vertices[i].ex;
  ey = m_vertices[i].ey;
  return true;
}
 
bool ComponentTcad2d::LoadData(const std::string& datafilename) {
  std::ifstream datafile;
  datafile.open(datafilename.c_str(), std::ios::in);
  if (!datafile) {
    std::cerr << m_className << "::LoadData:\n"
              << "    Could not open file " << datafilename << ".\n";
    return false;
  }
 
  const unsigned int nVertices = m_vertices.size();
  std::vector<unsigned int> fillCount(nVertices, 0);
  for (unsigned int i = 0; i < nVertices; ++i) {
    m_vertices[i].p = 0.;
    m_vertices[i].ex = 0.;
    m_vertices[i].ey = 0.;
    m_vertices[i].emob = 0.;
    m_vertices[i].hmob = 0.;
    m_vertices[i].donorOcc.clear();
    m_vertices[i].acceptorOcc.clear();
  }
  m_donors.clear();
  m_acceptors.clear();
 
  while (!datafile.fail()) {
    // Read one line and strip white space from the beginning of the line.
    std::string line;
    std::getline(datafile, line);
    ltrim(line);
    // Find data section.
    if (line.substr(0, 8) != "function") continue;
    // Read type of data set.
    const std::string::size_type pEq = line.find('=');
    if (pEq == std::string::npos) {
      // No "=" found.
      std::cerr << m_className << "::LoadData:\n"
                << "    Error reading file " << datafilename << ".\n"
                << "    Line:\n    " << line << "\n";
      datafile.close();
      Cleanup();
      return false;
    }
    line = line.substr(pEq + 1);
    std::string dataset;
    std::istringstream data;
    data.str(line);
    data >> dataset;
    if (dataset == "ElectrostaticPotential") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasPotential = true;
    } else if (dataset == "ElectricField") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasField = true;
    } else if (dataset == "eDriftVelocity") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasElectronVelocity = true;
    } else if (dataset == "hDriftVelocity") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasHoleVelocity = true;
    } else if (dataset == "eMobility") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasElectronMobility = true;
    } else if (dataset == "hMobility") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasHoleMobility = true;
    } else if (dataset == "eLifetime") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasElectronLifetime = true;
    } else if (dataset == "hLifetime") {
      if (!ReadDataset(datafile, dataset)) return false;
      m_hasHoleLifetime = true;
    } else if (dataset.substr(0, 14) == "TrapOccupation" &&
               dataset.substr(17, 2) == "Do") {
      if (!ReadDataset(datafile, dataset)) return false;
      Defect donor;
      donor.xsece = -1.;
      donor.xsech = -1.;
      donor.conc = -1.;
      m_donors.push_back(donor);
    } else if (dataset.substr(0, 14) == "TrapOccupation" &&
               dataset.substr(17, 2) == "Ac") {
      if (!ReadDataset(datafile, dataset)) return false;
      Defect acceptor;
      acceptor.xsece = -1.;
      acceptor.xsech = -1.;
      acceptor.conc = -1.;
      m_acceptors.push_back(acceptor);
    }
  }
  if (datafile.fail() && !datafile.eof()) {
    std::cerr << m_className << "::LoadData:\n"
              << "    Error reading file " << datafilename << "\n";
    datafile.close();
    Cleanup();
    return false;
  }
 
  datafile.close();
 
  return true;
}
 
bool ComponentTcad2d::ReadDataset(std::ifstream& datafile,
                                  const std::string& dataset) {
  enum DataSet {
    ElectrostaticPotential,
    EField,
    eDriftVelocity,
    hDriftVelocity,
    eMobility,
    hMobility,
    eLifetime,
    hLifetime,
    DonorTrapOccupation,
    AcceptorTrapOccupation,
    Unknown
  };
  DataSet ds = Unknown;
  if (dataset == "ElectrostaticPotential") {
    ds = ElectrostaticPotential;
  } else if (dataset == "ElectricField") {
    ds = EField;
  } else if (dataset == "eDriftVelocity") {
    ds = eDriftVelocity;
  } else if (dataset == "hDriftVelocity") {
    ds = hDriftVelocity;
  } else if (dataset == "eMobility") {
    ds = eMobility;
  } else if (dataset == "hMobility") {
    ds = hMobility;
  } else if (dataset == "eLifetime") {
    ds = eLifetime;
  } else if (dataset == "hLifetime") {
    ds = hLifetime;
  } else if (dataset.substr(0, 14) == "TrapOccupation") {
    if (dataset.substr(17, 2) == "Do") {
      ds = DonorTrapOccupation;
    } else if (dataset.substr(17, 2) == "Ac") {
      ds = AcceptorTrapOccupation;
    }
  } else {
    std::cerr << m_className << "::ReadDataset:\n"
              << "    Unexpected dataset " << dataset << ".\n";
    return false;
  }
 
  bool isVector = false;
  if (ds == EField || ds == eDriftVelocity || ds == hDriftVelocity) {
    isVector = true;
  }
 
  if (!datafile.is_open()) return false;
  std::string line;
  std::getline(datafile, line);
  std::getline(datafile, line);
  std::getline(datafile, line);
  std::getline(datafile, line);
  // Get the region name (given in brackets).
  std::string::size_type bra = line.find('[');
  std::string::size_type ket = line.find(']');
  if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
    std::cerr << m_className << "::ReadDataset:\n"
              << "    Cannot extract region name.\n"
              << "    Line:\n    " << line << "\n";
    datafile.close();
    Cleanup();
    return false;
  }
  line = line.substr(bra + 1, ket - bra - 1);
  std::string name;
  std::istringstream data;
  data.str(line);
  data >> name;
  data.clear();
  // Check if the region name matches one from the mesh file.
  const int index = FindRegion(name);
  if (index == -1) {
    std::cerr << m_className << "::ReadDataset:\n"
              << "    Unknown region " << name << ".\n";
    datafile.close();
    Cleanup();
    return false;
  }
  // Get the number of values.
  std::getline(datafile, line);
  bra = line.find('(');
  ket = line.find(')');
  if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
    std::cerr << m_className << "::ReadDataset:\n"
              << "    Cannot extract number of values to be read.\n"
              << "    Line:\n    " << line << "\n";
    datafile.close();
    Cleanup();
    return false;
  }
  line = line.substr(bra + 1, ket - bra - 1);
  int nValues;
  data.str(line);
  data >> nValues;
  data.clear();
  if (isVector) nValues = nValues / 2;
  // Mark the vertices belonging to this region.
  const unsigned int nVertices = m_vertices.size();
  std::vector<bool> isInRegion(nVertices, false);
  const unsigned int nElements = m_elements.size();
  for (unsigned int j = 0; j < nElements; ++j) {
    if (m_elements[j].region != index) continue;
    for (int k = 0; k <= m_elements[j].type; ++k) {
      isInRegion[m_elements[j].vertex[k]] = true;
    }
  }
 
  unsigned int ivertex = 0;
  for (int j = 0; j < nValues; ++j) {
    // Read the next value.
    double val1, val2;
    if (isVector) {
      datafile >> val1 >> val2;
    } else {
      datafile >> val1;
    }
    // Find the next vertex belonging to the region.
    while (ivertex < nVertices) {
      if (isInRegion[ivertex]) break;
      ++ivertex;
    }
    // Check if there is a mismatch between the number of vertices
    // and the number of values.
    if (ivertex >= nVertices) {
      std::cerr << m_className << "::ReadDataset:\n"
                << "    Dataset " << dataset
                << " has more values than vertices in region " << name << "\n";
      datafile.close();
      Cleanup();
      return false;
    }
    switch (ds) {
      case ElectrostaticPotential:
        m_vertices[ivertex].p = val1;
        break;
      case EField:
        m_vertices[ivertex].ex = val1;
        m_vertices[ivertex].ey = val2;
        break;
      case eDriftVelocity:
        // Scale from cm/s to cm/ns.
        m_vertices[ivertex].eVx = val1 * 1.e-9;
        m_vertices[ivertex].eVy = val2 * 1.e-9;
        break;
      case hDriftVelocity:
        // Scale from cm/s to cm/ns.
        m_vertices[ivertex].hVx = val1 * 1.e-9;
        m_vertices[ivertex].hVy = val2 * 1.e-9;
        break;
      case eMobility:
        // Convert from cm2 / (V s) to cm2 / (V ns).
        m_vertices[ivertex].emob = val1 * 1.e-9;
        break;
      case hMobility:
        // Convert from cm2 / (V s) to cm2 / (V ns).
        m_vertices[ivertex].hmob = val1 * 1.e-9;
        break;
      case eLifetime:
        // Convert from s to ns.
        m_vertices[ivertex].eTau = val1 * 1.e9;
        break;
      case hLifetime:
        // Convert from s to ns.
        m_vertices[ivertex].hTau = val1 * 1.e9;
        break;
      case DonorTrapOccupation:
        m_vertices[ivertex].donorOcc.push_back(val1);
        break;
      case AcceptorTrapOccupation:
        m_vertices[ivertex].acceptorOcc.push_back(val1);
        break;
      default:
        std::cerr << m_className << "::ReadDataset:\n"
                  << "    Unexpected dataset (" << ds << "). Program bug!\n";
        datafile.close();
        Cleanup();
        return false;
    }
    ++ivertex;
  }
  return true;
}
 
bool ComponentTcad2d::LoadWeightingField(const std::string& datafilename,
    std::vector<std::array<double, 2> >& wf, std::vector<double>& wp) {
  std::ifstream datafile;
  datafile.open(datafilename.c_str(), std::ios::in);
  if (!datafile) {
    std::cerr << m_className << "::LoadWeightingField:\n"
              << "    Could not open file " << datafilename << ".\n";
    return false;
  }
 
  const unsigned int nVertices = m_vertices.size();
  bool ok = true;
  while (!datafile.fail()) {
    // Read one line and strip white space from the beginning of the line.
    std::string line;
    std::getline(datafile, line);
    ltrim(line);
    // Find data section.
    if (line.substr(0, 8) != "function") continue;
    // Read type of data set.
    const std::string::size_type pEq = line.find('=');
    if (pEq == std::string::npos) {
      // No "=" found.
      std::cerr << m_className << "::LoadWeightingField:\n"
                << "    Error reading file " << datafilename << ".\n"
                << "    Line:\n    " << line << "\n";
      datafile.close();
      return false;
    }
    line = line.substr(pEq + 1);
    std::string dataset;
    std::istringstream data;
    data.str(line);
    data >> dataset;
    data.clear();
    if (dataset != "ElectrostaticPotential" && dataset != "ElectricField") {
      continue;
    }
    bool field = false;
    if (dataset == "ElectricField") {
      wf.assign(nVertices, {0., 0.});
      field = true;
    } else {
      wp.assign(nVertices, 0.);
    }
    std::getline(datafile, line);
    std::getline(datafile, line);
    std::getline(datafile, line);
    std::getline(datafile, line);
    // Get the region name (given in brackets).
    auto bra = line.find('[');
    auto ket = line.find(']');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      std::cerr << m_className << "::LoadWeightingField:\n"
                << "    Cannot extract region name.\n"
                << "    Line:\n    " << line << "\n";
      ok = false;
      break;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    std::string name;
    data.str(line);
    data >> name;
    data.clear();
    // Check if the region name matches one from the mesh file.
    const int index = FindRegion(name);
    if (index == -1) {
      std::cerr << m_className << "::LoadWeightingField:\n"
                << "    Unknown region " << name << ".\n";
      ok = false;
      break;
    }
    // Get the number of values.
    std::getline(datafile, line);
    bra = line.find('(');
    ket = line.find(')');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      std::cerr << m_className << "::LoadWeightingField:\n"
                << "    Cannot extract number of values to be read.\n"
                << "    Line:\n    " << line << "\n";
      ok = false;
      break;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    int nValues;
    data.str(line);
    data >> nValues;
    data.clear();
    if (field) nValues /= 2;
    // Mark the vertices belonging to this region.
    std::vector<bool> isInRegion(nVertices, false);
    const unsigned int nElements = m_elements.size();
    for (unsigned int j = 0; j < nElements; ++j) {
      if (m_elements[j].region != index) continue;
      for (int k = 0; k <= m_elements[j].type; ++k) {
        isInRegion[m_elements[j].vertex[k]] = true;
      }
    }
    unsigned int ivertex = 0;
    for (int j = 0; j < nValues; ++j) {
      // Read the next value.
      double val1, val2;
      if (field) {
        datafile >> val1 >> val2;
      } else {
        datafile >> val1;
      }
      // Find the next vertex belonging to the region.
      while (ivertex < nVertices) {
        if (isInRegion[ivertex]) break;
        ++ivertex;
      }
      // Check if there is a mismatch between the number of vertices
      // and the number of values.
      if (ivertex >= nVertices) {
        std::cerr << m_className << "::LoadWeightingField:\n"
                  << "    Dataset " << dataset
                  << " has more values than vertices in region " << name << "\n";
        ok = false;
        break;
      }
      if (field) {
        wf[ivertex] = {val1, val2};
      } else {
        wp[ivertex] = val1;
      }
      ++ivertex;
    }
  }
 
  if (!ok || (datafile.fail() && !datafile.eof())) {
    std::cerr << m_className << "::LoadWeightingField:\n"
              << "    Error reading file " << datafilename << "\n";
    datafile.close();
    return false;
  }
 
  datafile.close();
  return true;
}
 
 
 
bool ComponentTcad2d::LoadGrid(const std::string& gridfilename) {
  // Open the file containing the mesh description.
  std::ifstream gridfile;
  gridfile.open(gridfilename.c_str(), std::ios::in);
  if (!gridfile) {
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Could not open file " << gridfilename << ".\n";
    return false;
  }
 
  // Delete existing mesh information.
  Cleanup();
  // Count line numbers.
  unsigned int iLine = 0;
  // Get the number of regions.
  unsigned int nRegions = 0;
  while (!gridfile.fail()) {
    // Read one line and strip white space from the beginning of the line.
    std::string line;
    std::getline(gridfile, line);
    ++iLine;
    ltrim(line);
    // Find entry 'nb_regions'.
    if (line.substr(0, 10) != "nb_regions") continue;
    const std::string::size_type pEq = line.find('=');
    if (pEq == std::string::npos) {
      // No "=" sign found.
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Could not read number of regions.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(pEq + 1);
    std::istringstream data;
    data.str(line);
    data >> nRegions;
    break;
  }
  if (gridfile.eof()) {
    // Reached end of file.
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Could not find entry 'nb_regions' in file\n"
              << "    " << gridfilename << ".\n";
    Cleanup();
    gridfile.close();
    return false;
  } else if (gridfile.fail()) {
    // Error reading from the file.
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Error reading file " << gridfilename << " (line " << iLine
              << ").\n";
    Cleanup();
    gridfile.close();
    return false;
  }
  m_regions.resize(nRegions);
  for (unsigned int j = 0; j < nRegions; ++j) {
    m_regions[j].name = "";
    m_regions[j].drift = false;
    m_regions[j].medium = nullptr;
  }
 
  if (m_debug) {
    std::cout << m_className << "::LoadGrid:\n"
              << "    Found " << nRegions << " regions.\n";
  }
 
  // Get the region names.
  while (!gridfile.fail()) {
    std::string line;
    std::getline(gridfile, line);
    ++iLine;
    ltrim(line);
    // Find entry 'regions'.
    if (line.substr(0, 7) != "regions") continue;
    // Get region names (given in brackets).
    const std::string::size_type bra = line.find('[');
    const std::string::size_type ket = line.find(']');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      // No closed brackets [].
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Could not read region names.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    std::istringstream data;
    data.str(line);
    for (unsigned int j = 0; j < nRegions; ++j) {
      data >> m_regions[j].name;
      data.clear();
      // Assume by default that all regions are active.
      m_regions[j].drift = true;
      m_regions[j].medium = nullptr;
    }
    break;
  }
  if (gridfile.eof()) {
    // Reached end of file.
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Could not find entry 'regions' in file\n"
              << "    " << gridfilename << ".\n";
    Cleanup();
    gridfile.close();
    return false;
  } else if (gridfile.fail()) {
    // Error reading from the file.
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Error reading file " << gridfilename << " (line " << iLine
              << ").\n";
    Cleanup();
    gridfile.close();
    return false;
  }
 
  // Get the vertices.
  unsigned int nVertices = 0;
  while (!gridfile.fail()) {
    std::string line;
    std::getline(gridfile, line);
    ++iLine;
    ltrim(line);
    // Find section 'Vertices'.
    if (line.substr(0, 8) != "Vertices") continue;
    // Get number of vertices (given in brackets).
    const std::string::size_type bra = line.find('(');
    const std::string::size_type ket = line.find(')');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      // No closed brackets [].
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Could not read number of vertices.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    std::istringstream data;
    data.str(line);
    data >> nVertices;
    m_vertices.resize(nVertices);
    // Get the coordinates of this vertex.
    for (unsigned int j = 0; j < nVertices; ++j) {
      gridfile >> m_vertices[j].x >> m_vertices[j].y;
      // Change units from micron to cm.
      m_vertices[j].x *= 1.e-4;
      m_vertices[j].y *= 1.e-4;
      ++iLine;
    }
    break;
  }
  if (gridfile.eof()) {
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Could not find section 'Vertices' in file\n"
              << "    " << gridfilename << ".\n";
    Cleanup();
    gridfile.close();
    return false;
  } else if (gridfile.fail()) {
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Error reading file " << gridfilename << " (line " << iLine
              << ").\n";
    Cleanup();
    gridfile.close();
    return false;
  }
 
  // Get the "edges" (lines connecting two vertices).
  int nEdges = 0;
  // Temporary arrays for storing edge points.
  std::vector<int> edgeP1;
  std::vector<int> edgeP2;
  while (!gridfile.fail()) {
    std::string line;
    std::getline(gridfile, line);
    ++iLine;
    ltrim(line);
    // Find section 'Edges'.
    if (line.substr(0, 5) != "Edges") continue;
    // Get the number of edges (given in brackets).
    const std::string::size_type bra = line.find('(');
    const std::string::size_type ket = line.find(')');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      // No closed brackets ()
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Could not read number of edges.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    std::istringstream data;
    data.str(line);
    data >> nEdges;
    edgeP1.resize(nEdges);
    edgeP2.resize(nEdges);
    // Get the indices of the two endpoints.
    for (int j = 0; j < nEdges; ++j) {
      gridfile >> edgeP1[j] >> edgeP2[j];
      ++iLine;
    }
    break;
  }
  if (gridfile.eof()) {
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Could not find section 'Edges' in file\n"
              << "    " << gridfilename << ".\n";
    Cleanup();
    gridfile.close();
    return false;
  } else if (gridfile.fail()) {
    std::cerr << m_className << "::LoadGrid:\n"
              << "    Error reading file " << gridfilename << " (line " << iLine
              << ").\n";
    Cleanup();
    gridfile.close();
    return false;
  }
 
  for (int i = 0; i < nEdges; ++i) {
    // Make sure the indices of the edge endpoints are not out of range.
    if (edgeP1[i] < 0 || edgeP1[i] >= (int)nVertices || edgeP2[i] < 0 ||
        edgeP2[i] >= (int)nVertices) {
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Vertex index of edge " << i << " out of range.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    // Make sure the edge is non-degenerate.
    if (edgeP1[i] == edgeP2[i]) {
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Edge " << i << " is degenerate.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
  }
 
  // Get the elements.
  int edge0, edge1, edge2, edge3;
  int type;
  unsigned int nElements = 0;
  while (!gridfile.fail()) {
    std::string line;
    std::getline(gridfile, line);
    ++iLine;
    ltrim(line);
    // Find section 'Elements'.
    if (line.substr(0, 8) != "Elements") continue;
    // Get number of elements (given in brackets).
    const std::string::size_type bra = line.find('(');
    const std::string::size_type ket = line.find(')');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      // No closed brackets ().
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Could not read number of elements.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    std::istringstream data;
    data.str(line);
    data >> nElements;
    // Resize array of elements.
    m_elements.resize(nElements);
    // Get type and constituting edges of each element.
    for (unsigned int j = 0; j < nElements; ++j) {
      for (int k = nMaxVertices; k--;) m_elements[j].vertex[k] = -1;
      m_elements[j].neighbours.clear();
      ++iLine;
      gridfile >> type;
      switch (type) {
        case 1:
          // Line
          gridfile >> edge0 >> edge1;
          if (edge0 < 0) edge0 = -edge0 - 1;
          if (edge1 < 0) edge1 = -edge1 - 1;
          // Make sure the indices are not out of range.
          if (edge0 >= nEdges || edge1 >= nEdges) {
            std::cerr << m_className << "::LoadGrid:\n"
                      << "    Error reading file " << gridfilename << " (line "
                      << iLine << ").\n"
                      << "    Edge index out of range.\n";
            Cleanup();
            gridfile.close();
            return false;
          }
          // Get the vertices of this element.
          // Negative edge index means that the sequence of the two points
          // is supposed to be inverted.
          // The actual index is then given by "-index - 1".
          // Orientt the line such that the first point is on the left.
          if (m_vertices[edgeP1[edge0]].x > m_vertices[edgeP2[edge0]].x) {
            m_elements[j].vertex[0] = edgeP2[edge0];
            m_elements[j].vertex[1] = edgeP1[edge0];
          } else {
            m_elements[j].vertex[0] = edgeP1[edge0];
            m_elements[j].vertex[1] = edgeP2[edge0];
          }
          break;
        case 2:
          // Triangle
          gridfile >> edge0 >> edge1 >> edge2;
          // Make sure the indices are not out of range.
          if (edge0 < 0) edge0 = -edge0 - 1;
          if (edge1 < 0) edge1 = -edge1 - 1;
          if (edge2 < 0) edge2 = -edge2 - 1;
          if (edge0 >= nEdges || edge1 >= nEdges || edge2 >= nEdges) {
            std::cerr << m_className << "::LoadGrid:\n"
                      << "    Error reading file " << gridfilename << " (line "
                      << iLine << ").\n"
                      << "    Edge index out of range.\n";
            Cleanup();
            gridfile.close();
            return false;
          }
          m_elements[j].vertex[0] = edgeP1[edge0];
          m_elements[j].vertex[1] = edgeP2[edge0];
          if (edgeP1[edge1] != m_elements[j].vertex[0] &&
              edgeP1[edge1] != m_elements[j].vertex[1]) {
            m_elements[j].vertex[2] = edgeP1[edge1];
          } else {
            m_elements[j].vertex[2] = edgeP2[edge1];
          }
          // Rearrange vertices such that point 0 is on the left.
          while (m_vertices[m_elements[j].vertex[0]].x >
                     m_vertices[m_elements[j].vertex[1]].x ||
                 m_vertices[m_elements[j].vertex[0]].x >
                     m_vertices[m_elements[j].vertex[2]].x) {
            const int tmp = m_elements[j].vertex[0];
            m_elements[j].vertex[0] = m_elements[j].vertex[1];
            m_elements[j].vertex[1] = m_elements[j].vertex[2];
            m_elements[j].vertex[2] = tmp;
          }
          break;
        case 3:
          // Rectangle
          gridfile >> edge0 >> edge1 >> edge2 >> edge3;
          // Make sure the indices are not out of range.
          if (edge0 >= nEdges || -edge0 - 1 >= nEdges || edge1 >= nEdges ||
              -edge1 - 1 >= nEdges || edge2 >= nEdges || -edge2 - 1 >= nEdges ||
              edge3 >= nEdges || -edge3 - 1 >= nEdges) {
            std::cerr << m_className << "::LoadGrid:\n"
                      << "    Error reading file " << gridfilename << " (line "
                      << iLine << ").\n"
                      << "    Edge index out of range.\n";
            Cleanup();
            gridfile.close();
            return false;
          }
          if (edge0 >= 0)
            m_elements[j].vertex[0] = edgeP1[edge0];
          else
            m_elements[j].vertex[0] = edgeP2[-edge0 - 1];
          if (edge1 >= 0)
            m_elements[j].vertex[1] = edgeP1[edge1];
          else
            m_elements[j].vertex[1] = edgeP2[-edge1 - 1];
          if (edge2 >= 0)
            m_elements[j].vertex[2] = edgeP1[edge2];
          else
            m_elements[j].vertex[2] = edgeP2[-edge2 - 1];
          if (edge3 >= 0)
            m_elements[j].vertex[3] = edgeP1[edge3];
          else
            m_elements[j].vertex[3] = edgeP2[-edge3 - 1];
 
          // Rearrange vertices such that point 0 is on the left.
          while (m_vertices[m_elements[j].vertex[0]].x >
                     m_vertices[m_elements[j].vertex[1]].x ||
                 m_vertices[m_elements[j].vertex[0]].x >
                     m_vertices[m_elements[j].vertex[2]].x ||
                 m_vertices[m_elements[j].vertex[0]].x >
                     m_vertices[m_elements[j].vertex[3]].x) {
            const int tmp = m_elements[j].vertex[0];
            m_elements[j].vertex[0] = m_elements[j].vertex[1];
            m_elements[j].vertex[1] = m_elements[j].vertex[2];
            m_elements[j].vertex[2] = m_elements[j].vertex[3];
            m_elements[j].vertex[3] = tmp;
          }
          break;
        default:
          // Other element types are not permitted for 2d grids.
          std::cerr << m_className << "::LoadGrid:\n"
                    << "    Error reading file " << gridfilename << " (line "
                    << iLine << ").\n";
          std::cerr << "    Invalid element type (" << type
                    << ") for 2d mesh.\n";
          Cleanup();
          gridfile.close();
          return false;
          break;
      }
      m_elements[j].type = type;
      m_elements[j].region = -1;
    }
    break;
  }
  if (gridfile.eof()) {
    std::cerr << m_className << "::LoadGrid:\n";
    std::cerr << "    Could not find section 'Elements' in file\n";
    std::cerr << "    " << gridfilename << ".\n";
    Cleanup();
    gridfile.close();
    return false;
  } else if (gridfile.fail()) {
    std::cerr << m_className << "::LoadGrid:\n";
    std::cerr << "    Error reading file " << gridfilename << " (line " << iLine
              << ").\n";
    Cleanup();
    gridfile.close();
    return false;
  }
 
  // Assign regions to elements.
  std::string name;
  while (!gridfile.fail()) {
    std::string line;
    std::getline(gridfile, line);
    ltrim(line);
    // Find section 'Region'.
    if (line.substr(0, 6) != "Region") continue;
    // Get region name (given in brackets).
    std::string::size_type bra = line.find('(');
    std::string::size_type ket = line.find(')');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Could not read region name.\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    std::istringstream data;
    data.str(line);
    data >> name;
    data.clear();
    const int index = FindRegion(name);
    if (index == -1) {
      // Specified region name is not in the list.
      std::cerr << m_className << "::LoadGrid:\n"
                << "    Error reading file " << gridfilename << ".\n"
                << "    Unknown region " << name << ".\n";
      continue;
    }
    std::getline(gridfile, line);
    std::getline(gridfile, line);
    bra = line.find('(');
    ket = line.find(')');
    if (ket < bra || bra == std::string::npos || ket == std::string::npos) {
      // No closed brackets ().
      std::cerr << m_className << "::LoadGrid:\n";
      std::cerr << "    Error reading file " << gridfilename << ".\n";
      std::cerr << "    Could not read number of elements in region " << name
                << ".\n";
      Cleanup();
      gridfile.close();
      return false;
    }
    line = line.substr(bra + 1, ket - bra - 1);
    int nElementsRegion;
    int iElement;
    data.str(line);
    data >> nElementsRegion;
    data.clear();
    for (int j = 0; j < nElementsRegion; ++j) {
      gridfile >> iElement;
      m_elements[iElement].region = index;
    }
  }
 
  gridfile.close();
  if (gridfile.fail() && !gridfile.eof()) {
    std::cerr << m_className << "::LoadGrid:\n";
    std::cerr << "    Error reading file " << gridfilename << ".\n";
    Cleanup();
    return false;
  }
 
  return true;
}
 
void ComponentTcad2d::FindNeighbours() {
  const unsigned int nElements = m_elements.size();
  std::vector<std::vector<bool> > adjacent(nElements,
                                           std::vector<bool>(nElements, false));
 
  const double tol = 5.e-4;
  for (unsigned int i = 0; i < nElements; ++i) {
    const Element& ei = m_elements[i];
    for (unsigned int j = 0; j < nElements; ++j) {
      if (i == j || adjacent[i][j]) continue;
      const Element& ej = m_elements[j];
      if (ei.xmin > ej.xmax + tol || ei.xmax < ej.xmin - tol) continue;
      if (ei.ymin > ej.ymax + tol || ei.ymax < ej.ymin - tol) continue;
      for (unsigned int m = 0; m < nMaxVertices; ++m) {
        if (ei.vertex[m] < 0) break;
        for (unsigned int n = 0; n < nMaxVertices; ++n) {
          if (ei.vertex[n] < 0) break;
          if (ei.vertex[m] == ej.vertex[n]) {
            adjacent[i][j] = adjacent[j][i] = true;
            break;
          }
        }
        if (adjacent[i][j]) break;
      }
    }
  }
 
  for (unsigned int i = 0; i < nElements; ++i) {
    m_elements[i].neighbours.clear();
    for (unsigned int j = 0; j < nElements; ++j) {
      if (adjacent[i][j]) {
        m_elements[i].neighbours.push_back(j);
      }
    }
  }
}
 
void ComponentTcad2d::Cleanup() {
  // Vertices
  m_vertices.clear();
  // Elements
  m_elements.clear();
  // Regions
  m_regions.clear();
  // Weighting fields and potentials
  m_wf.clear();
  m_wp.clear();
}
 
unsigned int ComponentTcad2d::FindElement(
    const double x, const double y, 
    std::array<double, nMaxVertices>& w) const {
 
  w.fill(0.);
  if (m_lastElement >= 0) {
    // Check if the point is still located in the previously found element.
    const Element& last = m_elements[m_lastElement];
    if (x >= last.xmin && x <= last.xmax && y >= last.ymin && y <= last.ymax) {
      if (CheckElement(x, y, last, w)) return m_lastElement;
    }
    // The point is not in the previous element.
    // Check the adjacent elements.
    const unsigned int nNeighbours = last.neighbours.size();
    for (unsigned int i = 0; i < nNeighbours; ++i) {
      const Element& element = m_elements[last.neighbours[i]];
      if (x < element.xmin || x > element.xmax || y < element.ymin ||
          y > element.ymax)
        continue;
      if (!CheckElement(x, y, element, w)) continue;
      return last.neighbours[i];
    }
  }
 
  // The point is not in the previous element nor in the adjacent ones.
  // We have to loop over all elements.
  const unsigned int nElements = m_elements.size();
  for (unsigned int i = 0; i < nElements; ++i) {
    const Element& element = m_elements[i];
    if (x < element.xmin || x > element.xmax || y < element.ymin ||
        y > element.ymax)
      continue;
    if (CheckElement(x, y, element, w)) return i;
  }
  // Point is outside the mesh.
  if (m_debug) {
    std::cerr << m_className << "::FindElement:\n"
              << "    Point (" << x << ", " << y << ") is outside the mesh.\n";
  }
  return nElements;
}
 
bool ComponentTcad2d::CheckElement(const double x, const double y,
                                   const Element& element,
                                   std::array<double, nMaxVertices>& w) const {
  switch (element.type) {
    case 1:
      return CheckLine(x, y, element, w);
      break;
    case 2:
      return CheckTriangle(x, y, element, w);
      break;
    case 3:
      return CheckRectangle(x, y, element, w);
      break;
    default:
      std::cerr << m_className << "::CheckElement:\n"
                << "    Unknown element type. Program bug!\n";
      break;
  }
  return false;
}
 
bool ComponentTcad2d::CheckRectangle(const double x, const double y,
                                     const Element& element,
                                     std::array<double, nMaxVertices>& w) const {
  const Vertex& v0 = m_vertices[element.vertex[0]];
  const Vertex& v1 = m_vertices[element.vertex[1]];
  const Vertex& v3 = m_vertices[element.vertex[3]];
  if (y < v0.y || x > v3.x || y > v1.y) return false;
 
  // Map (x, y) to local variables (u, v) in [-1, 1].
  const double u = (x - 0.5 * (v0.x + v3.x)) / (v3.x - v0.x);
  const double v = (y - 0.5 * (v0.y + v1.y)) / (v1.y - v0.y);
  // Compute weighting factors for each corner.
  w[0] = (0.5 - u) * (0.5 - v);
  w[1] = (0.5 - u) * (0.5 + v);
  w[2] = (0.5 + u) * (0.5 + v);
  w[3] = (0.5 + u) * (0.5 - v);
  return true;
}
 
bool ComponentTcad2d::CheckTriangle(const double x, const double y,
                                    const Element& element,
                                    std::array<double, nMaxVertices>& w) const {
  const Vertex& v0 = m_vertices[element.vertex[0]];
  const Vertex& v1 = m_vertices[element.vertex[1]];
  const Vertex& v2 = m_vertices[element.vertex[2]];
  if (x > v1.x && x > v2.x) return false;
  if (y < v0.y && y < v1.y && y < v2.y) return false;
  if (y > v0.y && y > v1.y && y > v2.y) return false;
 
  // Map (x, y) onto local variables (b, c) such that
  // P = A + b * (B - A) + c * (C - A)
  // A point P is inside the triangle ABC if b, c > 0 and b + c < 1;
  // b, c are also weighting factors for points B, C
  const double v1x = v1.x - v0.x;
  const double v1y = v1.y - v0.y;
  const double v2x = v2.x - v0.x;
  const double v2y = v2.y - v0.y;
 
  w[1] = ((x - v0.x) * v2y - (y - v0.y) * v2x) / (v1x * v2y - v1y * v2x);
  if (w[1] < 0. || w[1] > 1.) return false;
  w[2] = ((v0.x - x) * v1y - (v0.y - y) * v1x) / (v1x * v2y - v1y * v2x);
  if (w[2] < 0. || w[1] + w[2] > 1.) return false;
 
  // Weighting factor for point A
  w[0] = 1. - w[1] - w[2];
 
  return true;
}
 
bool ComponentTcad2d::CheckLine(const double x, const double y,
                                const Element& element,
                                std::array<double, nMaxVertices>& w) const {
  const Vertex& v0 = m_vertices[element.vertex[0]];
  const Vertex& v1 = m_vertices[element.vertex[1]];
  if (x > v1.x) return false;
  if (y < v0.y && y < v1.y) return false;
  if (y > v0.y && y > v1.y) return false;
  const double tx = (x - v0.x) / (v1.x - v0.x);
  if (tx < 0. || tx > 1.) return false;
  const double ty = (y - v0.y) / (v1.y - v0.y);
  if (ty < 0. || ty > 1.) return false;
  if (tx == ty) {
    // Compute weighting factors for endpoints A, B
    w[0] = tx;
    w[1] = 1. - w[0];
    return true;
  }
  return false;
}
 
void ComponentTcad2d::Reset() {
  Cleanup();
  m_hasRangeZ = false;
  m_ready = false;
}
 
void ComponentTcad2d::UpdatePeriodicity() {
  if (!m_ready) {
    std::cerr << m_className << "::UpdatePeriodicity:\n"
              << "    Field map not available.\n";
    return;
  }
 
  // Check for conflicts.
  for (unsigned int i = 0; i < 3; ++i) {
    if (m_periodic[i] && m_mirrorPeriodic[i]) {
      std::cerr << m_className << "::UpdatePeriodicity:\n"
                << "    Both simple and mirror periodicity requested. Reset.\n";
      m_periodic[i] = m_mirrorPeriodic[i] = false;
    }
  }
 
  if (m_axiallyPeriodic[0] || m_axiallyPeriodic[1] || m_axiallyPeriodic[2]) {
    std::cerr << m_className << "::UpdatePeriodicity:\n"
              << "    Axial symmetry is not supported. Reset.\n";
    m_axiallyPeriodic.fill(false);
  }
 
  if (m_rotationSymmetric[0] || m_rotationSymmetric[1] ||
      m_rotationSymmetric[2]) {
    std::cerr << m_className << "::UpdatePeriodicity:\n"
              << "    Rotation symmetry is not supported. Reset.\n";
    m_rotationSymmetric.fill(false);
  }
}
 
int ComponentTcad2d::FindRegion(const std::string& name) const {
  const unsigned int nRegions = m_regions.size();
  for (unsigned int j = 0; j < nRegions; ++j) {
    if (name == m_regions[j].name) return j;
  }
  return -1;
}
 
void ComponentTcad2d::MapCoordinates(double& x, double& y, bool& xmirr,
                                     bool& ymirr) const {
  // In case of periodicity, reduce to the cell volume.
  xmirr = false;
  const double cellsx = m_xMaxBB - m_xMinBB;
  if (m_periodic[0]) {
    x = m_xMinBB + fmod(x - m_xMinBB, cellsx);
    if (x < m_xMinBB) x += cellsx;
  } else if (m_mirrorPeriodic[0]) {
    double xNew = m_xMinBB + fmod(x - m_xMinBB, cellsx);
    if (xNew < m_xMinBB) xNew += cellsx;
    const int nx = int(floor(0.5 + (xNew - x) / cellsx));
    if (nx != 2 * (nx / 2)) {
      xNew = m_xMinBB + m_xMaxBB - xNew;
      xmirr = true;
    }
    x = xNew;
  }
  ymirr = false;
  const double cellsy = m_yMaxBB - m_yMinBB;
  if (m_periodic[1]) {
    y = m_yMinBB + fmod(y - m_yMinBB, cellsy);
    if (y < m_yMinBB) y += cellsy;
  } else if (m_mirrorPeriodic[1]) {
    double yNew = m_yMinBB + fmod(y - m_yMinBB, cellsy);
    if (yNew < m_yMinBB) yNew += cellsy;
    const int ny = int(floor(0.5 + (yNew - y) / cellsy));
    if (ny != 2 * (ny / 2)) {
      yNew = m_yMinBB + m_yMaxBB - yNew;
      ymirr = true;
    }
    y = yNew;
  }
}
 
bool ComponentTcad2d::CheckTraps() const {
  for (const auto& trap : m_donors) {
    if (trap.xsece < 0. || trap.xsech < 0. || trap.conc < 0.) return false;
  }
 
  for (const auto& trap : m_acceptors) {
    if (trap.xsece < 0. || trap.xsech < 0. || trap.conc < 0.) return false;
  }
  return true;
}
}