Quadratic.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include <algorithm>
#include <boost/multi_array.hpp>
#include <complex>
#include <fstream>
#include <iomanip>
#include <list>
#include <map>
#include <set>
#include <sstream>
#include <stack>
#include <vector>
 
#include "MantidGeometry/Math/PolyBase.h"
#include "MantidGeometry/Math/mathSupport.h"
#include "MantidGeometry/Surfaces/BaseVisit.h"
#include "MantidGeometry/Surfaces/Quadratic.h"
#include "MantidGeometry/Surfaces/Surface.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/Tolerance.h"
#include "MantidKernel/V3D.h"
 
#include "MantidKernel/Strings.h"
 
namespace Mantid::Geometry {
namespace {
Kernel::Logger logger("Quadratic");
}
using Kernel::Tolerance;
using Kernel::V3D;
 
Quadratic::Quadratic()
    : Surface(), BaseEqn(10)
{}
 
double Quadratic::eqnValue(const Kernel::V3D &Pt) const
{
  double res(0.0);
  res += BaseEqn[0] * Pt[0] * Pt[0];
  res += BaseEqn[1] * Pt[1] * Pt[1];
  res += BaseEqn[2] * Pt[2] * Pt[2];
  res += BaseEqn[3] * Pt[0] * Pt[1];
  res += BaseEqn[4] * Pt[0] * Pt[2];
  res += BaseEqn[5] * Pt[1] * Pt[2];
  res += BaseEqn[6] * Pt[0];
  res += BaseEqn[7] * Pt[1];
  res += BaseEqn[8] * Pt[2];
  res += BaseEqn[9];
  return res;
}
 
int Quadratic::side(const Kernel::V3D &Pt) const
{
  double res = eqnValue(Pt);
  if (fabs(res) < Tolerance)
    return 0;
  return (res > 0) ? 1 : -1;
}
 
Kernel::V3D Quadratic::surfaceNormal(const Kernel::V3D &Pt) const
{
  Kernel::V3D N(2 * BaseEqn[0] * Pt[0] + BaseEqn[3] * Pt[1] + BaseEqn[4] * Pt[2] + BaseEqn[6],
                2 * BaseEqn[1] * Pt[1] + BaseEqn[3] * Pt[0] + BaseEqn[5] * Pt[2] + BaseEqn[7],
                2 * BaseEqn[2] * Pt[2] + BaseEqn[4] * Pt[0] + BaseEqn[5] * Pt[1] + BaseEqn[8]);
  N.normalize();
  return N;
}
 
void Quadratic::matrixForm(Kernel::Matrix<double> &A, Kernel::V3D &B, double &C) const
{
  A.setMem(3, 3); // set incase memory out
  for (int i = 0; i < 3; i++)
    A[i][i] = BaseEqn[i];
 
  A[0][1] = A[1][0] = BaseEqn[3] / 2.0;
  A[0][2] = A[2][0] = BaseEqn[4] / 2.0;
  A[1][2] = A[2][1] = BaseEqn[5] / 2.0;
 
  for (int i = 0; i < 3; i++)
    B[i] = BaseEqn[6 + i];
  C = BaseEqn[9];
}
 
double Quadratic::distance(const Kernel::V3D &Pt) const
{
  // Job 1 :: Create matrix and vector representation
  Kernel::Matrix<double> A(3, 3);
  Kernel::V3D B;
  double cc;
  matrixForm(A, B, cc);
 
  // Job 2 :: calculate the diagonal matrix
  Kernel::Matrix<double> D(3, 3);
  Kernel::Matrix<double> R(3, 3);
  if (!A.Diagonalise(R, D)) {
    logger.warning() << "Problem with matrix :: distance now guessed at\n";
    return distance(Pt);
  }
 
  Kernel::V3D alpha = R.Tprime() * Pt;
  Kernel::V3D beta = R.Tprime() * B;
 
  // Calculate fundermental equation:
  const double aa(alpha[0]);
  const double aa2(aa * aa);
  const double ab(alpha[1]);
  const double ab2(ab * ab);
  const double ac(alpha[2]);
  const double ac2(ac * ac);
 
  const double ba(beta[0]);
  const double ba2(ba * ba);
  const double bb(beta[1]);
  const double bb2(bb * bb);
  const double bc(beta[2]);
  const double bc2(bc * bc);
 
  const double da(D[0][0]);
  const double da2(da * da);
  const double db(D[1][1]);
  const double db2(db * db);
  const double dc(D[2][2]);
  const double dc2(dc * dc);
 
  mathLevel::PolyBase T(6);
 
  T[0] = aa * ba + ab * bb + ac * bc + cc + aa2 * da + ab2 * db + ac2 * dc;
 
  T[1] = -ba2 - bb2 - bc2 + 4 * ab * bb * da + 4 * ac * bc * da + 4 * cc * da + 4 * aa * ba * db + 4 * ac * bc * db +
         4 * cc * db + 4 * aa2 * da * db + 4 * ab2 * da * db + 4 * aa * ba * dc + 4 * ab * bb * dc + 4 * cc * dc +
         4 * aa2 * da * dc + 4 * ac2 * da * dc + 4 * ab2 * db * dc + 4 * ac2 * db * dc;
 
  T[2] = -ba2 * da - 4 * bb2 * da - 4 * bc2 * da + 4 * ab * bb * da2 + 4 * ac * bc * da2 + 4 * cc * da2 - 4 * ba2 * db -
         bb2 * db - 4 * bc2 * db + 16 * ac * bc * da * db + 16 * cc * da * db + 4 * ab2 * da2 * db + 4 * aa * ba * db2 +
         4 * ac * bc * db2 + 4 * cc * db2 + 4 * aa2 * da * db2 - 4 * ba2 * dc - 4 * bb2 * dc - bc2 * dc +
         16 * ab * bb * da * dc + 16 * cc * da * dc + 4 * ac2 * da2 * dc + 16 * aa * ba * db * dc + 16 * cc * db * dc +
         16 * aa2 * da * db * dc + 16 * ab2 * da * db * dc + 16 * ac2 * da * db * dc + 4 * ac2 * db2 * dc +
         4 * aa * ba * dc2 + 4 * ab * bb * dc2 + 4 * cc * dc2 + 4 * aa2 * da * dc2 + 4 * ab2 * db * dc2;
 
  T[3] = -4 * bb2 * da2 - 4 * bc2 * da2 - 4 * ba2 * da * db - 4 * bb2 * da * db - 16 * bc2 * da * db +
         16 * ac * bc * da2 * db + 16 * cc * da2 * db - 4 * ba2 * db2 - 4 * bc2 * db2 + 16 * ac * bc * da * db2 +
         16 * cc * da * db2 - 4 * ba2 * da * dc - 16 * bb2 * da * dc - 4 * bc2 * da * dc + 16 * ab * bb * da2 * dc +
         16 * cc * da2 * dc - 16 * ba2 * db * dc - 4 * bb2 * db * dc - 4 * bc2 * db * dc + 64 * cc * da * db * dc +
         16 * ab2 * da2 * db * dc + 16 * ac2 * da2 * db * dc + 16 * aa * ba * db2 * dc + 16 * cc * db2 * dc +
         16 * aa2 * da * db2 * dc + 16 * ac2 * da * db2 * dc - 4 * ba2 * dc2 - 4 * bb2 * dc2 + 16 * ab * bb * da * dc2 +
         16 * cc * da * dc2 + 16 * aa * ba * db * dc2 + 16 * cc * db * dc2 + 16 * aa2 * da * db * dc2 +
         16 * ab2 * da * db * dc2;
 
  T[4] = -4 * bb2 * da2 * db - 16 * bc2 * da2 * db - 4 * ba2 * da * db2 - 16 * bc2 * da * db2 +
         16 * ac * bc * da2 * db2 + 16 * cc * da2 * db2 - 16 * bb2 * da2 * dc - 4 * bc2 * da2 * dc -
         16 * ba2 * da * db * dc - 16 * bb2 * da * db * dc - 16 * bc2 * da * db * dc + 64 * cc * da2 * db * dc -
         16 * ba2 * db2 * dc - 4 * bc2 * db2 * dc + 64 * cc * da * db2 * dc + 16 * ac2 * da2 * db2 * dc -
         4 * ba2 * da * dc2 - 16 * bb2 * da * dc2 + 16 * ab * bb * da2 * dc2 + 16 * cc * da2 * dc2 -
         16 * ba2 * db * dc2 - 4 * bb2 * db * dc2 + 64 * cc * da * db * dc2 + 16 * ab2 * da2 * db * dc2 +
         16 * aa * ba * db2 * dc2 + 16 * cc * db2 * dc2 + 16 * aa2 * da * db2 * dc2;
 
  T[5] = -16 * bc2 * da2 * db2 - 16 * bb2 * da2 * db * dc - 16 * bc2 * da2 * db * dc - 16 * ba2 * da * db2 * dc -
         16 * bc2 * da * db2 * dc + 64 * cc * da2 * db2 * dc - 16 * bb2 * da2 * dc2 - 16 * ba2 * da * db * dc2 -
         16 * bb2 * da * db * dc2 + 64 * cc * da2 * db * dc2 - 16 * ba2 * db2 * dc2 + 64 * cc * da * db2 * dc2;
 
  T[6] = 16 * da * db * dc * (-bc2 * da * db - bb2 * da * dc - ba2 * db * dc + 4 * cc * da * db * dc);
 
  std::vector<double> TRange = T.realRoots(1e-10);
  if (TRange.empty())
    return -1.0;
 
  double Out = -1;
  Kernel::V3D xvec;
  std::vector<double>::const_iterator vc;
  for (vc = TRange.begin(); vc != TRange.end(); ++vc) {
    const double daI = 1.0 + 2 * (*vc) * da;
    const double dbI = 1.0 + 2 * (*vc) * db;
    const double dcI = 1.0 + 2 * (*vc) * dc;
    if ((daI * daI) > Tolerance || ((dbI * dbI) > Tolerance && (dcI * dcI) < Tolerance)) {
      Kernel::Matrix<double> DI(3, 3);
      DI[0][0] = 1.0 / daI;
      DI[1][1] = 1.0 / dbI;
      DI[2][2] = 1.0 / dcI;
      xvec = R * (DI * (alpha - beta * *vc)); // care here: to avoid 9*9 +9*3 in
                                              // favour of 9*3+9*3 ops.
      const double Dist = xvec.distance(Pt);
      if (Out < 0 || Out > Dist)
        Out = Dist;
    }
  }
  return Out;
}
 
bool Quadratic::onSurface(const Kernel::V3D &Pt) const
{
  const double res = eqnValue(Pt);
  return (std::abs(res) <= Tolerance);
}
 
void Quadratic::displace(const Kernel::V3D &Pt)
{
  BaseEqn[9] += Pt[0] * (Pt[0] * BaseEqn[0] - BaseEqn[6]) + Pt[1] * (Pt[1] * BaseEqn[1] - BaseEqn[7]) +
                Pt[2] * (Pt[2] * BaseEqn[2] - BaseEqn[8]) + BaseEqn[4] * Pt[0] * Pt[1] + BaseEqn[5] * Pt[0] * Pt[2] +
                BaseEqn[6] * Pt[1] * Pt[2];
  BaseEqn[6] += -2 * BaseEqn[0] * Pt[0] - BaseEqn[3] * Pt[1] - BaseEqn[4] * Pt[2];
  BaseEqn[7] += -2 * BaseEqn[1] * Pt[1] - BaseEqn[3] * Pt[0] - BaseEqn[5] * Pt[2];
  BaseEqn[8] += -2 * BaseEqn[2] * Pt[2] - BaseEqn[4] * Pt[0] - BaseEqn[5] * Pt[1];
}
 
void Quadratic::rotate(const Kernel::Matrix<double> &MX)
{
  Kernel::Matrix<double> MA = MX;
  MA.Invert();
  const double a(MA[0][0]), b(MA[0][1]), c(MA[0][2]);
  const double d(MA[1][0]), e(MA[1][1]), f(MA[1][2]);
  const double g(MA[2][0]), h(MA[2][1]), j(MA[2][2]);
  double B[9];
  B[0] = BaseEqn[0] * a * a + BaseEqn[1] * d * d + BaseEqn[2] * g * g + BaseEqn[3] * a * d + BaseEqn[4] * a * g +
         BaseEqn[5] * d * g;
 
  B[1] = BaseEqn[0] * b * b + BaseEqn[1] * e * e + BaseEqn[2] * h * h + BaseEqn[3] * b * e + BaseEqn[4] * b * h +
         BaseEqn[5] * e * h;
 
  B[2] = BaseEqn[0] * c * c + BaseEqn[1] * f * f + BaseEqn[2] * j * j + BaseEqn[3] * c * f + BaseEqn[4] * c * j +
         BaseEqn[5] * f * j;
 
  B[3] = 2.0 * (BaseEqn[0] * a * b + BaseEqn[1] * d * e + BaseEqn[2] * g * h) + BaseEqn[3] * (a * e + b * d) +
         BaseEqn[4] * (a * h + b * g) + BaseEqn[5] * (d * j + e * g);
 
  B[4] = 2.0 * (BaseEqn[0] * a * c + BaseEqn[1] * d * f + BaseEqn[2] * g * j) + BaseEqn[3] * (a * f + c * d) +
         BaseEqn[4] * (a * j + c * h) + BaseEqn[5] * (d * j + f * h);
 
  B[5] = 2.0 * (BaseEqn[0] * b * c + BaseEqn[1] * e * f + BaseEqn[2] * h * j) + BaseEqn[3] * (b * f + c * e) +
         BaseEqn[4] * (b * j + c * h) + BaseEqn[5] * (e * j + f * h);
 
  B[6] = BaseEqn[6] * a + BaseEqn[7] * d + BaseEqn[8] * g;
 
  B[7] = BaseEqn[6] * b + BaseEqn[7] * e + BaseEqn[8] * h;
 
  B[8] = BaseEqn[6] * c + BaseEqn[7] * f + BaseEqn[8] * j;
 
  for (int i = 0; i < 9; i++) // Item 9 left alone
    BaseEqn[i] = B[i];
}
 
void Quadratic::print() const
{
  Surface::print();
  for (int i = 0; i < 10; i++)
    logger.debug() << BaseEqn[i] << " ";
  logger.debug() << '\n';
}
 
void Quadratic::write(std::ostream &OX) const
{
  std::ostringstream cx;
  writeHeader(cx);
  cx.precision(Surface::Nprecision);
  cx << "gq ";
  for (int i = 0; i < 4; i++)
    cx << " " << BaseEqn[i] << " ";
  cx << " " << BaseEqn[5] << " ";
  cx << " " << BaseEqn[4] << " ";
  for (int i = 6; i < 10; i++)
    cx << " " << BaseEqn[i] << " ";
  Mantid::Kernel::Strings::writeMCNPX(cx.str(), OX);
}
 
std::unique_ptr<Quadratic> Quadratic::clone() const { return std::unique_ptr<Quadratic>(doClone()); }
 
} // namespace Mantid::Geometry