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