src/INTERP_KERNEL/VectorUtils.hxx

   1 // Copyright (C) 2007-2022  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    * Caracteristic vector size (its biggest component, in absolute)
 164    */
 165   inline double caracteristicDimVector(const double *v)
 166   {
 167     double ret = 0;
 168     for (int i = 0; i < 3; i++)
 169       ret = std::max(ret, std::fabs(v[i]));
 170     return ret;
 171   }
 172
 173   /**
 174    * Compares doubles using a relative tolerance
 175    * This is suitable mainly for comparing larger values to each other. Before performing the relative test,
 176    * an absolute test is performed to guard from problems when comparing to 0.0
 177    *
 178    * @param x         first value
 179    * @param y         second value
 180    * @param relTol    maximum allowed relative difference that is to be treated as equality
 181    * @param absTol    maximum allowed absolute difference that is to be treated as equality
 182    * @return  true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise
 183    */
 184   inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
 185   {
 186     // necessary for comparing values close to zero
 187     // in order to avoid division by very small numbers
 188     if(std::fabs(x - y) < absTol)
 189       {
 190         return true;
 191       }
 192
 193     const double relError = std::fabs((x - y) / std::max(std::fabs(x), std::fabs(y)));
 194
 195     return relError < relTol;
 196   }
 197
 198   inline double sumOfAbsoluteValues(const double row[3])
 199   {
 200     double ret(0.);
 201     std::for_each(row,row+3,[&ret](double v) { ret += std::abs(v); });
 202     return ret;
 203   }
 204
 205   /*!
 206    * Returns the infinite norm of a 3x3 input matrix \a mat.
 207    * The max of absolute value of row sum.
 208    */
 209   inline double normInf(const double mat[9])
 210   {
 211     double ret(std::max(sumOfAbsoluteValues(mat),sumOfAbsoluteValues(mat+3)));
 212     return std::max(ret,sumOfAbsoluteValues(mat+6));
 213   }
 214
 215 }
 216
 217 #endif