Integration3D.cpp

#ifndef INTEGRATION3D_CC
#define INTEGRATION3D_CC

#include "Integration3D.h"
#include "recurCut.h"
#include "../../Numeric/Gauss.h"

#define ZERO_LS_TOL  1.e-7
#define EQUALITY_TOL 1.e-15

// cross product
inline void cross(const double *v1, const double *v2, double *v) {
  v[0] = v1[1] * v2[2] - v2[1] * v1[2];
  v[1] = v2[0] * v1[2] - v1[0] * v2[2];
  v[2] = v1[0] * v2[1] - v2[0] * v1[1];
}
// dot product
inline double dot(const double *v1, const double *v2) {
  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}
// norm
inline double norm(const double *v) {
  return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}
// substract
inline void vec (const double *x1, const double *x2, double *v) {
  v[0] = x2[0] - x1[0]; v[1] = x2[1] - x1[1]; v[2] = x2[2] - x1[2];
}
inline void vec (const DI_Point &p1, const DI_Point &p2, double *v) {
  v[0] = p2.x() - p1.x(); v[1] = p2.y() - p1.y(); v[2] = p2.z() - p1.z();
}
inline double sq2 (const double a) {return a * a;}

// determinant
inline double det3(double d11, double d12, double d13,
                   double d21, double d22, double d23,
                   double d31, double d32, double d33) {
  return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22);
}
inline double det4 (double d11, double d12, double d13, double d14,
                    double d21, double d22, double d23, double d24,
                    double d31, double d32, double d33, double d34,
                    double d41, double d42, double d43, double d44) {
  return d11 * det3(d22, d23, d24, d32, d33, d34, d42, d43, d44)
       - d21 * det3(d12, d13, d14, d32, d33, d34, d42, d43, d44)
       + d31 * det3(d12, d13, d14, d22, d23, d24, d42, d43, d44)
       - d41 * det3(d12, d13, d14, d22, d23, d24, d32, d33, d34);
}

// distance
inline double distance(const DI_Point &p1, const DI_Point &p2) {
  return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x())
            + (p1.y() - p2.y()) * (p1.y() - p2.y())
            + (p1.z() - p2.z()) * (p1.z() - p2.z()));
}
inline double distance(const DI_CuttingPoint &p1, const DI_CuttingPoint &p2) {
  return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x())
            + (p1.y() - p2.y()) * (p1.y() - p2.y())
            + (p1.z() - p2.z()) * (p1.z() - p2.z()));
}

// middle of 2 points
inline DI_Point middle (const DI_Point &p1, const DI_Point &p2) {
  return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
}
inline DI_Point middle (const DI_Point &p1, const DI_Point &p2, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
  return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2, e, RPNi);
}

// barycentre