src/INTERP_KERNEL/VectorUtils.hxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   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 <numeric>
  26 #include <string>
  27 #include <cmath>
  28 #include <map>
  29
  30
  31 namespace INTERP_KERNEL
  32 {
  33   /// Precision used for tests of 3D part of INTERP_KERNEL
  34   const double VOL_PREC = 1.0e-6;
  35
  36   /// Default relative tolerance in epsilonEqualRelative
  37   const double DEFAULT_REL_TOL = 1.0e-6;
  38
  39   /// Default absolute tolerance in epsilonEqual and epsilonEqualRelative
  40   const double DEFAULT_ABS_TOL = 5.0e-12;
  41
  42   /**
  43    * @param a first point. Should point on a array of size at least equal to SPACEDIM.
  44    * @param b second point. Should point on a array of size at least equal to SPACEDIM.
  45    */
  46   template<int SPACEDIM>
  47   inline double getDistanceBtw2Pts(const double *a, const double *b)
  48   {
  49     double ret2=0.;
  50     for(int i=0;i<SPACEDIM;i++)
  51       ret2+=(a[i]-b[i])*(a[i]-b[i]);
  52     return sqrt(ret2);
  53   }
  54
  55   // -------------------------------------------------------------------
  56   // Math operations for vectors represented by double[3] - arrays
  57   // -------------------------------------------------------------------
  58
  59   /**
  60    * Copies a double[3] vector from src to dest
  61    *
  62    * @param src   source vector
  63    * @param dest  destination vector
  64    *
  65    */
  66   inline void copyVector3(const double* src, double* dest)
  67   {
  68     for(int i = 0 ; i < 3 ; ++i)
  69       dest[i] = src[i];
  70   }
  71
  72   /**
  73    * Creates a string representation of a double[3] vector
  74    *
  75    * @param  pt  a 3-vector
  76    * @return a string of the form [x, y, z]
  77    */
  78   inline const std::string vToStr(const double* pt)
  79   {
  80     std::stringstream ss(std::ios::out);
  81     ss << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
  82     return ss.str();
  83   }
  84
  85   /**
  86    * Adds a double[3] - vector to another one.
  87    *
  88    * @param v     vector v
  89    * @param res   vector in which to store the result res + v.
  90    */
  91   inline void add(const double* v, double* res)
  92   {
  93     res[0] += v[0];
  94     res[1] += v[1];
  95     res[2] += v[2];
  96   }
  97
  98   /**
  99    * Calculates the cross product of two double[3] - vectors.
 100    *
 101    * @param v1    vector v1
 102    * @param v2    vector v2
 103    * @param res   vector in which to store the result v1 x v2. It should not be one of v1 and v2.
 104    */
 105   inline void cross(const double* v1, const double* v2,double* res)
 106   {
 107     res[0] = v1[1]*v2[2] - v1[2]*v2[1];
 108     res[1] = v1[2]*v2[0] - v1[0]*v2[2];
 109     res[2] = v1[0]*v2[1] - v1[1]*v2[0];
 110   }
 111
 112   /**
 113    * Calculates the dot product of two double[3] - vectors
 114    *
 115    * @param v1   vector v1
 116    * @param v2   vector v2
 117    * @return   dot (scalar) product v1.v2
 118    */
 119   inline double dot(const double* v1, const double* v2)
 120   {
 121     return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
 122   }
 123
 124   /**
 125    * Calculates norm of a double[3] vector
 126    *
 127    * @param v  a vector v
 128    * @return euclidean norm of v
 129    */
 130   inline double norm(const double* v)
 131   {
 132     return sqrt(dot(v,v));
 133   }
 134
 135   /**
 136    * Compares doubles using an absolute tolerance
 137    * This is suitable mainly for comparisons with 0.0
 138    *
 139    * @param x         first value
 140    * @param y         second value
 141    * @param errTol    maximum allowed absolute difference that is to be treated as equality
 142    * @return  true if |x - y| < errTol, false otherwise
 143    */
 144   inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL)
 145   {
 146     return y < x ? x - y < errTol : y - x < errTol;
 147     //    return std::fabs(x - y) < errTol;
 148   }
 149
 150
 151   /**
 152    * Test whether two 3D vectors are colinear. The two vectors are expected to be of unit norm (not checked)
 153    * Implemented by checking that the norm of the cross product is null.
 154    */
 155   inline bool isColinear3D(const double *v1, const double *v2, const double eps = DEFAULT_ABS_TOL)
 156   {
 157     double cros[3];
 158     cross(v1, v2, cros);
 159     return epsilonEqual(dot(cros, cros), 0.0, eps);
 160   }
 161
 162
 163   /**
 164    * Compares doubles using a relative tolerance
 165    * This is suitable mainly for comparing larger values to each other. Before performing the relative test,
 166    * an absolute test is performed to guard from problems when comparing to 0.0
 167    *
 168    * @param x         first value
 169    * @param y         second value
 170    * @param relTol    maximum allowed relative difference that is to be treated as equality
 171    * @param absTol    maximum allowed absolute difference that is to be treated as equality
 172    * @return  true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise
 173    */
 174   inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
 175   {
 176     // necessary for comparing values close to zero
 177     // in order to avoid division by very small numbers
 178     if(std::fabs(x - y) < absTol)
 179       {
 180         return true;
 181       }
 182
 183     const double relError = std::fabs((x - y) / std::max(std::fabs(x), std::fabs(y)));
 184
 185     return relError < relTol;
 186   }
 187
 188   inline double sumOfAbsoluteValues(const double row[3])
 189   {
 190     double ret(0.);
 191     std::for_each(row,row+3,[&ret](double v) { ret += std::abs(v); });
 192     return ret;
 193   }
 194
 195   /*!
 196    * Returns the infinite norm of a 3x3 input matrix \a mat.
 197    * The max of absolute value of row sum.
 198    */
 199   inline double normInf(const double mat[9])
 200   {
 201     double ret(std::max(sumOfAbsoluteValues(mat),sumOfAbsoluteValues(mat+3)));
 202     return std::max(ret,sumOfAbsoluteValues(mat+6));
 203   }
 204
 205 }
 206
 207 #endif