src/INTERP_KERNEL/VectorUtils.hxx

   1 //  Copyright (C) 2007-2008  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.
   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 #ifndef __VECTORUTILS_HXX__
  20 #define __VECTORUTILS_HXX__
  21
  22 #include <sstream>
  23 #include <numeric>
  24 #include <string>
  25 #include <cmath>
  26 #include <map>
  27
  28 /// Precision used for tests of 3D part of INTERP_KERNEL
  29 #define VOL_PREC 1.0e-6
  30
  31 /// Default relative tolerance in epsilonEqualRelative
  32 #define DEFAULT_REL_TOL 1.0e-6
  33
  34 /// Default absolute tolerance in epsilonEqual and epsilonEqualRelative
  35 #define DEFAULT_ABS_TOL 5.0e-12
  36
  37 namespace INTERP_KERNEL
  38 {
  39   /**
  40    * @param a first point. Should point on a array of size at least equal to SPACEDIM.
  41    * @param b second point. Should point on a array of size at least equal to SPACEDIM.
  42    */
  43   template<int SPACEDIM>
  44   inline double getDistanceBtw2Pts(const double *a, const double *b)
  45   {
  46     double ret2=0.;
  47     for(int i=0;i<SPACEDIM;i++)
  48       ret2+=(a[i]-b[i])*(a[i]-b[i]);
  49     return sqrt(ret2);
  50   }
  51
  52   // -------------------------------------------------------------------
  53   // Math operations for vectors represented by double[3] - arrays
  54   // -------------------------------------------------------------------
  55
  56   /**
  57    * Copies a double[3] vector from src to dest
  58    *
  59    * @param src   source vector
  60    * @param dest  destination vector
  61    *
  62    */
  63   inline void copyVector3(const double* src, double* dest)
  64   {
  65     for(int i = 0 ; i < 3 ; ++i)
  66       dest[i] = src[i];
  67   }
  68
  69   /**
  70    * Creates a string representation of a double[3] vector
  71    *
  72    * @param  pt  a 3-vector
  73    * @return a string of the form [x, y, z]
  74    */
  75   inline const std::string vToStr(const double* pt)
  76   {
  77     std::stringstream ss(std::ios::out);
  78     ss << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
  79     return ss.str();
  80   }
  81
  82   /**
  83    * Calculates the cross product of two double[3] - vectors.
  84    *
  85    * @param v1    vector v1
  86    * @param v2    vector v2
  87    * @param res   vector in which to store the result v1 x v2. It should not be one of v1 and v2.
  88    */
  89   inline void cross(const double* v1, const double* v2,double* res)
  90   {
  91     res[0] = v1[1]*v2[2] - v1[2]*v2[1];
  92     res[1] = v1[2]*v2[0] - v1[0]*v2[2];
  93     res[2] = v1[0]*v2[1] - v1[1]*v2[0];
  94   }
  95
  96   /**
  97    * Calculates the dot product of two double[3] - vectors
  98    *
  99    * @param v1   vector v1
 100    * @param v2   vector v2
 101    * @return   dot (scalar) product v1.v2
 102    */
 103   inline double dot(const double* v1, const double* v2)
 104   {
 105     return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
 106   }
 107
 108   /**
 109    * Calculates norm of a double[3] vector
 110    *
 111    * @param v  a vector v
 112    * @return euclidean norm of v
 113    */
 114   inline double norm(const double* v)
 115   {
 116     return sqrt(dot(v,v));
 117   }
 118
 119   /**
 120    * Compares doubles using an absolute tolerance
 121    * This is suitable mainly for comparisons with 0.0
 122    *
 123    * @param x         first value
 124    * @param y         second value
 125    * @param errTol    maximum allowed absolute difference that is to be treated as equality
 126    * @return  true if |x - y| < errTol, false otherwise
 127    */
 128   inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL)
 129   {
 130     return y < x ? x - y < errTol : y - x < errTol;
 131     //    return std::fabs(x - y) < errTol;
 132   }
 133
 134   /**
 135    * Compares doubles using a relative tolerance
 136    * This is suitable mainly for comparing larger values to each other. Before performing the relative test,
 137    * an absolute test is performed to guard from problems when comparing to 0.0
 138    *
 139    * @param x         first value
 140    * @param y         second value
 141    * @param relTol    maximum allowed relative difference that is to be treated as equality
 142    * @param absTol    maximum allowed absolute difference that is to be treated as equality
 143    * @return  true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise
 144    */
 145   inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
 146   {
 147     // necessary for comparing values close to zero
 148     // in order to avoid division by very small numbers
 149     if(std::fabs(x - y) < absTol)
 150       {
 151         return true;
 152       }
 153
 154     const double relError = std::fabs((x - y) / std::max(std::fabs(x), std::fabs(y)));
 155
 156     return relError < relTol;
 157   }
 158
 159 }
 160
 161 #endif