1 /*! 2 * \file src/NUMODIS/Utilities.cxx 3 * \brief 4 * \author Laurent Dupuy 5 * \date 9/06/2017 6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 7 * reserved. 8 * This project is publicly released under either the GNU GPL Licence 9 * or the CECILL-A licence. A copy of thoses licences are delivered 10 * with the sources of TFEL. CEA or EDF may also distribute this 11 * project under specific licensing conditions. 12 */ 13 14 #include <cerrno> 15 #include <numeric> 16 #include <ostream> 17 #include <algorithm> 18 19 #include "NUMODIS/Math/Utilities.hxx" 20 21 namespace numodis 22 { 23 24 namespace math 25 { 26 27 //=============================================================== 28 // dUnitVector 29 //--------------------------------------------------------------- 30 //! Returns the unit vector V corresponding to a vector U (3D) 31 //--------------------------------------------------------------- 32 // Note that we do not check whether U is of non-zero length. 33 // \param U input vector (3D) 34 // \param V input vector (3D) 35 //=============================================================== dUnitVector(const Vect3 & U,Vect3 & V)36 double dUnitVector(const Vect3& U, 37 Vect3& V) 38 { 39 40 double norm=U.Length(); 41 42 V=U/norm; 43 44 return norm; 45 46 } 47 48 //=============================================================== 49 // dUnitVector 50 //--------------------------------------------------------------- 51 //! Normalize vector U (3D) 52 //--------------------------------------------------------------- 53 // Note that we do not check whether U is of non-zero length. 54 // \param U input/output vector (3D) 55 //=============================================================== dUnitVector(Vect3 & U)56 double dUnitVector(Vect3& U) 57 { 58 59 double norm=U.Length(); 60 61 U/=norm; 62 63 return norm; 64 65 } 66 67 //=============================================================== 68 // iCollinear 69 //--------------------------------------------------------------- 70 //! Determine whether two vectors of int are collinear or not 71 //--------------------------------------------------------------- 72 /*! 73 74 \bug QUID DE VECTEURS NULS? QUI SE SERT DE CA ET POURQUOI ? 75 76 \param U first vector 77 \param V second vector 78 \return true if U and V are collinear, false otherwise 79 */ 80 //=============================================================== iCollinear(const std::vector<int> & U,const std::vector<int> & V)81 bool iCollinear(const std::vector<int>& U, 82 const std::vector<int>& V) 83 { 84 85 //---------------------------------------------------- 86 // dot products (compatible with 4 indices notations) 87 //---------------------------------------------------- 88 int UV=iScalProduct(U,V); 89 int UU=iScalProduct(U,U); 90 int VV=iScalProduct(V,V); 91 92 //------------------------- 93 // compare scalar products 94 //------------------------- 95 return (UU!=0 && VV!=0 ? (UU*VV==UV*UV) : (UU==VV)); 96 97 } 98 99 //=============================================================== 100 // iCollinear 101 //--------------------------------------------------------------- 102 //! Determine whether two std::vector<int> are collinear or not 103 //--------------------------------------------------------------- 104 /*! 105 \param U first vector 106 \param V second vector 107 \param pdirection +1/0/-1 108 \return true if U and V are collinear, false otherwise 109 */ 110 //=============================================================== iCollinear(const std::vector<int> & U,const std::vector<int> & V,int * pdirection)111 bool iCollinear(const std::vector<int>& U, 112 const std::vector<int>& V, 113 int* pdirection) 114 { 115 116 //---------------- 117 // initialization 118 //---------------- 119 *pdirection=0; 120 121 //-------------- 122 // compare size 123 //-------------- 124 if(U.size()!=V.size()) 125 return false; 126 127 //---------------------------------------------------- 128 // dot products (compatible with 4 indices notations) 129 //---------------------------------------------------- 130 int UV=iScalProduct(U,V); 131 int UU=iScalProduct(U,U); 132 int VV=iScalProduct(V,V); 133 134 //------------------------- 135 // compare scalar products 136 //------------------------- 137 if(UU*VV==UV*UV) 138 { 139 *pdirection=(UV<0 ? -1 : +1); 140 return true; 141 } 142 143 //--------------- 144 // default value 145 //--------------- 146 return false; 147 148 } 149 150 //=============================================================== 151 // iSortVector 152 //--------------------------------------------------------------- 153 //! Sort a vector of int 154 //--------------------------------------------------------------- 155 /*! \param U vector to be sorted then sorted */ 156 //=============================================================== iSortVector(std::vector<int> & U)157 void iSortVector(std::vector<int>& U) 158 { 159 160 std::vector<int> V(U.size()); 161 162 // copy as a vector 163 for(unsigned i=0; i<U.size(); i++) 164 V[i]=U[i]; 165 166 // sort as a vector 167 sort(V.begin(),V.end()); 168 169 // copy back 170 for(unsigned i=0; i<U.size(); i++) 171 U[i]=V[i]; 172 173 } 174 175 //=============================================================== 176 // iSortVector3FirstValue 177 //--------------------------------------------------------------- 178 //! Sort 3 first values of a vector of int 179 //--------------------------------------------------------------- 180 /*! \param U vector to be sorted then sorted */ 181 //=============================================================== iSortVector3FirstValue(std::vector<int> & U)182 void iSortVector3FirstValue(std::vector<int>& U) 183 { 184 if(U[0]<U[1]) { std::swap(U[0],U[1]); } 185 if(U[1]<U[2]) { std::swap(U[1],U[2]); } 186 else return; 187 if(U[0]<U[1]) { std::swap(U[0],U[1]); } 188 } 189 190 //=============================================================== 191 // GCD 192 //--------------------------------------------------------------- 193 //! Compute the greatest common divisor of an ensemble of integers 194 //--------------------------------------------------------------- 195 /*! 196 \param u input vector of int 197 \return greatest common divisor 198 */ 199 //=============================================================== GCD(const std::vector<int> & u)200 int GCD(const std::vector<int>& u) 201 { 202 // trivial situation 203 if(u.size()==1) 204 return std::abs(u[0]); 205 // default value 206 int gcd=std::abs(u[0]); 207 // find the overall gcd 208 for(unsigned i=1; i<u.size(); i++) 209 gcd=numodis::math::GCD(gcd,std::abs(u[i])); 210 // output value 211 return gcd; 212 } 213 214 //=============================================================== 215 // Epsilon 216 //--------------------------------------------------------------- 217 //! Permutation / Levi Civita symbol 218 //--------------------------------------------------------------- 219 /*! 220 221 This function should work for {0,1,2} and {1,2,3}. 222 223 \param i i 224 \param j j 225 \param k k 226 \return permutation(i,j,k) (-1, 0 or +1) 227 228 */ 229 //=============================================================== Epsilon(const unsigned i,const unsigned j,const unsigned k)230 int Epsilon(const unsigned i, 231 const unsigned j, 232 const unsigned k) 233 { 234 235 // look for identical terms 236 if(i==j || j==k || i==k) 237 return 0; 238 239 // look for odd permutations 240 if(i+2==j || j+2==k || k+2==i) 241 return -1; 242 243 // default value 244 return 1; 245 246 } 247 248 } // end namespace math 249 250 } // end of namespace numodis 251