src/INTERP_KERNEL/VectorUtils.hxx

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