src/INTERP_KERNEL/VectorUtils.hxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #ifndef __VECTORUTILS_HXX__
  21 #define __VECTORUTILS_HXX__
  22
  23 #include <algorithm>
  24 #include <sstream>
  25 #include <iomanip>
  26 #include <numeric>
  27 #include <string>
  28 #include <cmath>
  29 #include <map>
  30
  31
  32 namespace INTERP_KERNEL
  33 {
  34   /// Precision used for tests of 3D part of INTERP_KERNEL
  35   const double VOL_PREC = 1.0e-6;
  36
  37   /// Default relative tolerance in epsilonEqualRelative
  38   const double DEFAULT_REL_TOL = 1.0e-6;
  39
  40   /// Default absolute tolerance in epsilonEqual and epsilonEqualRelative
  41   const double DEFAULT_ABS_TOL = 5.0e-12;
  42
  43   /**
  44    * @param a first point. Should point on a array of size at least equal to SPACEDIM.
  45    * @param b second point. Should point on a array of size at least equal to SPACEDIM.
  46    */
  47   template<int SPACEDIM>
  48   inline double getDistanceBtw2Pts(const double *a, const double *b)
  49   {
  50     double ret2=0.;
  51     for(int i=0;i<SPACEDIM;i++)
  52       ret2+=(a[i]-b[i])*(a[i]-b[i]);
  53     return sqrt(ret2);
  54   }
  55
  56   // -------------------------------------------------------------------
  57   // Math operations for vectors represented by double[3] - arrays
  58   // -------------------------------------------------------------------
  59
  60   /**
  61    * Copies a double[3] vector from src to dest
  62    *
  63    * @param src   source vector
  64    * @param dest  destination vector
  65    *
  66    */
  67   inline void copyVector3(const double* src, double* dest)
  68   {
  69     for(int i = 0 ; i < 3 ; ++i)
  70       dest[i] = src[i];
  71   }
  72
  73   /**
  74    * Creates a string representation of a double[3] vector
  75    *
  76    * @param  pt  a 3-vector
  77    * @return a string of the form [x, y, z]
  78    */
  79   inline const std::string vToStr(const double* pt)
  80   {
  81     std::stringstream ss(std::ios::out);
  82     ss << std::setprecision(16) << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
  83     return ss.str();
  84   }
  85
  86   /**
  87    * Adds a double[3] - vector to another one.
  88    *
  89    * @param v     vector v
  90    * @param res   vector in which to store the result res + v.
  91    */
  92   inline void add(const double* v, double* res)
  93   {
  94     res[0] += v[0];
  95     res[1] += v[1];
  96     res[2] += v[2];
  97   }
  98
  99   /**
 100    * Calculates the cross product of two double[3] - vectors.
 101    *
 102    * @param v1    vector v1
 103    * @param v2    vector v2
 104    * @param res   vector in which to store the result v1 x v2. It should not be one of v1 and v2.
 105    */
 106   inline void cross(const double* v1, const double* v2,double* res)
 107   {
 108     res[0] = v1[1]*v2[2] - v1[2]*v2[1];
 109     res[1] = v1[2]*v2[0] - v1[0]*v2[2];
 110     res[2] = v1[0]*v2[1] - v1[1]*v2[0];
 111   }
 112
 113   /**
 114    * Calculates the dot product of two double[3] - vectors
 115    *
 116    * @param v1   vector v1
 117    * @param v2   vector v2
 118    * @return   dot (scalar) product v1.v2
 119    */
 120   inline double dot(const double* v1, const double* v2)
 121   {
 122     return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
 123   }
 124
 125   /**
 126    * Calculates norm of a double[3] vector
 127    *
 128    * @param v  a vector v
 129    * @return euclidean norm of v
 130    */
 131   inline double norm(const double* v)
 132   {
 133     return sqrt(dot(v,v));
 134   }
 135
 136   /**
 137    * Compares doubles using an absolute tolerance
 138    * This is suitable mainly for comparisons with 0.0
 139    *
 140    * @param x         first value
 141    * @param y         second value
 142    * @param errTol    maximum allowed absolute difference that is to be treated as equality
 143    * @return  true if |x - y| < errTol, false otherwise
 144    */
 145   inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL)
 146   {
 147     return y < x ? x - y < errTol : y - x < errTol;
 148     //    return std::fabs(x - y) < errTol;
 149   }
 150
 151
 152   /**
 153    * Test whether two 3D vectors are colinear. The two vectors are expected to be of unit norm (not checked)
 154    * Implemented by checking that the norm of the cross product is null.
 155    */
 156   inline bool isColinear3D(const double *v1, const double *v2, const double eps = DEFAULT_ABS_TOL)
 157   {
 158     double cros[3];
 159     cross(v1, v2, cros);
 160     return epsilonEqual(dot(cros, cros), 0.0, eps);
 161   }
 162
 163   /**
 164    * Test whether three 3D points are colinear.
 165    * Implemented by calling overload isColinear3D(v1,v2, eps)
 166    */
 167   inline bool isColinear3DPts(const double *p1, const double *p2, const double *p3, const double eps = DEFAULT_ABS_TOL)
 168   {
 169     // check if these points are colinear
 170     double vec1[3];
 171     // p2-p1
 172     std::transform(p2, p2 + 3, p1, vec1, std::minus<double>());
 173     double vec2[3];
 174     // p3-p2
 175     std::transform(p3, p3 + 3, p2, vec2, std::minus<double>());
 176     return isColinear3D(vec1, vec2, eps);
 177   }
 178
 179   /**
 180    * Caracteristic vector size (its biggest component, in absolute)
 181    */
 182   inline double caracteristicDimVector(const double *v)
 183   {
 184     double ret = 0;
 185     for (int i = 0; i < 3; i++)
 186       ret = std::max(ret, std::fabs(v[i]));
 187     return ret;
 188   }
 189
 190   /**
 191    * Compares doubles using a relative tolerance
 192    * This is suitable mainly for comparing larger values to each other. Before performing the relative test,
 193    * an absolute test is performed to guard from problems when comparing to 0.0
 194    *
 195    * @param x         first value
 196    * @param y         second value
 197    * @param relTol    maximum allowed relative difference that is to be treated as equality
 198    * @param absTol    maximum allowed absolute difference that is to be treated as equality
 199    * @return  true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise
 200    */
 201   inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
 202   {
 203     // necessary for comparing values close to zero
 204     // in order to avoid division by very small numbers
 205     if(std::fabs(x - y) < absTol)
 206       {
 207         return true;
 208       }
 209
 210     const double relError = std::fabs((x - y) / std::max(std::fabs(x), std::fabs(y)));
 211
 212     return relError < relTol;
 213   }
 214
 215   inline double sumOfAbsoluteValues(const double row[3])
 216   {
 217     double ret(0.);
 218     std::for_each(row,row+3,[&ret](double v) { ret += std::abs(v); });
 219     return ret;
 220   }
 221
 222   /*!
 223    * Returns the infinite norm of a 3x3 input matrix \a mat.
 224    * The max of absolute value of row sum.
 225    */
 226   inline double normInf(const double mat[9])
 227   {
 228     double ret(std::max(sumOfAbsoluteValues(mat),sumOfAbsoluteValues(mat+3)));
 229     return std::max(ret,sumOfAbsoluteValues(mat+6));
 230   }
 231
 232 }
 233
 234 #endif