1 // Copyright (C) 2007-2024 CEA, EDF
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #ifndef __VECTORUTILS_HXX__
21 #define __VECTORUTILS_HXX__
32 namespace INTERP_KERNEL
34 /// Precision used for tests of 3D part of INTERP_KERNEL
35 const double VOL_PREC = 1.0e-6;
37 /// Default relative tolerance in epsilonEqualRelative
38 const double DEFAULT_REL_TOL = 1.0e-6;
40 /// Default absolute tolerance in epsilonEqual and epsilonEqualRelative
41 const double DEFAULT_ABS_TOL = 5.0e-12;
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.
47 template<int SPACEDIM>
48 inline double getDistanceBtw2Pts(const double *a, const double *b)
51 for(int i=0;i<SPACEDIM;i++)
52 ret2+=(a[i]-b[i])*(a[i]-b[i]);
56 // -------------------------------------------------------------------
57 // Math operations for vectors represented by double[3] - arrays
58 // -------------------------------------------------------------------
61 * Copies a double[3] vector from src to dest
63 * @param src source vector
64 * @param dest destination vector
67 inline void copyVector3(const double* src, double* dest)
69 for(int i = 0 ; i < 3 ; ++i)
74 * Creates a string representation of a double[3] vector
76 * @param pt a 3-vector
77 * @return a string of the form [x, y, z]
79 inline const std::string vToStr(const double* pt)
81 std::stringstream ss(std::ios::out);
82 ss << std::setprecision(16) << "[" << pt[0] << ", " << pt[1] << ", " << pt[2] << "]";
87 * Adds a double[3] - vector to another one.
90 * @param res vector in which to store the result res + v.
92 inline void add(const double* v, double* res)
100 * Calculates the cross product of two double[3] - vectors.
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.
106 inline void cross(const double* v1, const double* v2,double* res)
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];
114 * Calculates the dot product of two double[3] - vectors
116 * @param v1 vector v1
117 * @param v2 vector v2
118 * @return dot (scalar) product v1.v2
120 inline double dot(const double* v1, const double* v2)
122 return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
126 * Calculates norm of a double[3] vector
128 * @param v a vector v
129 * @return euclidean norm of v
131 inline double norm(const double* v)
133 return sqrt(dot(v,v));
137 * Compares doubles using an absolute tolerance
138 * This is suitable mainly for comparisons with 0.0
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
145 inline bool epsilonEqual(const double x, const double y, const double errTol = DEFAULT_ABS_TOL)
147 return y < x ? x - y < errTol : y - x < errTol;
148 // return std::fabs(x - y) < errTol;
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.
156 inline bool isColinear3D(const double *v1, const double *v2, const double eps = DEFAULT_ABS_TOL)
160 return epsilonEqual(dot(cros, cros), 0.0, eps);
164 * Test whether three 3D points are colinear.
165 * Implemented by calling overload isColinear3D(v1,v2, eps)
167 inline bool isColinear3DPts(const double *p1, const double *p2, const double *p3, const double eps = DEFAULT_ABS_TOL)
169 // check if these points are colinear
172 std::transform(p2, p2 + 3, p1, vec1, std::minus<double>());
175 std::transform(p3, p3 + 3, p2, vec2, std::minus<double>());
176 return isColinear3D(vec1, vec2, eps);
180 * Caracteristic vector size (its biggest component, in absolute)
182 inline double caracteristicDimVector(const double *v)
185 for (int i = 0; i < 3; i++)
186 ret = std::max(ret, std::fabs(v[i]));
191 * Compares doubles using a relative tolerance
192 * This is suitable mainly for comparing larger values to each other. Before performing the relative test,
193 * an absolute test is performed to guard from problems when comparing to 0.0
195 * @param x first value
196 * @param y second value
197 * @param relTol maximum allowed relative difference that is to be treated as equality
198 * @param absTol maximum allowed absolute difference that is to be treated as equality
199 * @return true if |x - y| <= absTol or |x - y|/max(|x|,|y|) <= relTol, false otherwise
201 inline bool epsilonEqualRelative(const double x, const double y, const double relTol = DEFAULT_REL_TOL, const double absTol = DEFAULT_ABS_TOL)
203 // necessary for comparing values close to zero
204 // in order to avoid division by very small numbers
205 if(std::fabs(x - y) < absTol)
210 const double relError = std::fabs((x - y) / std::max(std::fabs(x), std::fabs(y)));
212 return relError < relTol;
215 inline double sumOfAbsoluteValues(const double row[3])
218 std::for_each(row,row+3,[&ret](double v) { ret += std::abs(v); });
223 * Returns the infinite norm of a 3x3 input matrix \a mat.
224 * The max of absolute value of row sum.
226 inline double normInf(const double mat[9])
228 double ret(std::max(sumOfAbsoluteValues(mat),sumOfAbsoluteValues(mat+3)));
229 return std::max(ret,sumOfAbsoluteValues(mat+6));