1 // ============================================================================= 2 // PROJECT CHRONO - http://projectchrono.org 3 // 4 // Copyright (c) 2014 projectchrono.org 5 // All rights reserved. 6 // 7 // Use of this source code is governed by a BSD-style license that can be found 8 // in the LICENSE file at the top level of the distribution and at 9 // http://projectchrono.org/license-chrono.txt. 10 // 11 // ============================================================================= 12 // Authors: Alessandro Tasora 13 // ============================================================================= 14 15 #ifndef CHQUADRATURE 16 #define CHQUADRATURE 17 18 #include <vector> 19 20 #include "chrono/core/ChApiCE.h" 21 #include "chrono/core/ChMatrix.h" 22 23 namespace chrono { 24 25 /// Class to store polynomial roots and weights for the Gauss-Legendre 26 /// quadrature. It is automatically managed by ChQuadrature. 27 class ChApi ChQuadratureTables { 28 public: 29 ChQuadratureTables(int order_from = 1, int order_to = 10); 30 31 std::vector<std::vector<double> > Weight; 32 std::vector<std::vector<double> > Lroots; 33 34 void PrintTables(); 35 36 private: 37 void glege_roots(ChMatrixDynamic<>& lcoef, int N, int ntable); 38 }; 39 40 41 42 /// Class to store polynomial roots and weights for quadrature 43 /// over a triangle. Assumes 2 natural 'area' coordinates u v, 44 /// the 3rd assumed 1-u-v. 45 /// D. A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle, Int. J. Num. 46 /// Meth. Engng, 21(1985):1129-1148. 47 class ChApi ChQuadratureTablesTriangle { 48 public: 49 ChQuadratureTablesTriangle(); 50 51 std::vector<std::vector<double> > Weight; 52 std::vector<std::vector<double> > LrootsU; 53 std::vector<std::vector<double> > LrootsV; 54 }; 55 56 /// Class to store polynomial roots and weights for quadrature 57 /// over a tetrahedron. Assumes 3 natural 'volume' coordinates u v w, 58 /// the 4th assumed 1-u-v-w 59 class ChApi ChQuadratureTablesTetrahedron { 60 public: 61 ChQuadratureTablesTetrahedron(); 62 63 std::vector<std::vector<double> > Weight; 64 std::vector<std::vector<double> > LrootsU; 65 std::vector<std::vector<double> > LrootsV; 66 std::vector<std::vector<double> > LrootsW; 67 }; 68 69 70 /// Base class for 1D integrand T=f(x) to be used in ChQuadrature. 71 /// Since the class is templated, the computed valued can be either 72 /// a simple 'double' or a more complex object, like ChMatrixNM<..>. 73 /// You must inherit your custom class from this base class, and 74 /// implement your own Evaluate() method, for example: 75 /// 76 /// class MySine : public ChIntegrable1D<double> 77 /// { 78 /// public: 79 /// void Evaluate (double& result, const double x) { result = sin(x); } 80 /// }; 81 /// 82 template <class T = double> 83 class ChIntegrable1D { 84 public: ~ChIntegrable1D()85 virtual ~ChIntegrable1D() {} 86 87 /// Evaluate the function at point x , that is 88 /// result T = f(x) 89 virtual void Evaluate(T& result, const double x) = 0; 90 }; 91 92 /// As ChIntegrable1D, but for 2D integrand T=f(x,y) 93 /// to be used in ChQuadrature. 94 template <class T = double> 95 class ChIntegrable2D { 96 public: ~ChIntegrable2D()97 virtual ~ChIntegrable2D() {} 98 99 /// Evaluate the function at point x,y , that is 100 /// result T = f(x,y) 101 virtual void Evaluate(T& result, const double x, const double y) = 0; 102 }; 103 104 /// As ChIntegrable1D, but for 3D integrand T=f(x,y,z) 105 /// to be used in ChQuadrature. 106 template <class T = double> 107 class ChIntegrable3D { 108 public: ~ChIntegrable3D()109 virtual ~ChIntegrable3D() {} 110 111 /// Evaluate the function at point x,y,z , that is 112 /// result T = f(x,y,z) 113 virtual void Evaluate(T& result, const double x, const double y, const double z) = 0; 114 }; 115 116 /// Class to perform Gauss-Legendre quadrature, in 1D, 2D, 3D. 117 /// It integrates a function on a nD domain using the Gauss quadrature; this is mostly 118 /// useful in the case that the integrand is polynomial, because the result is exact 119 /// if the order of quadrature is greater or equal to the degree of the polynomial. 120 /// Integrate() functions are static; no need to allocate an instance of this class. 121 122 class ChApi ChQuadrature { 123 public: 124 /// Integrate the integrand T = f(x) over the 1D interval [xA, xB], 125 /// with desired order of quadrature. Best if integrand is polynomial. 126 /// For order in 1..10, precomputed polynomial coefficients are used for max speed. 127 template <class T> Integrate1D(T & result,ChIntegrable1D<T> & integrand,const double a,const double b,const int order)128 static void Integrate1D(T& result, ///< result is returned here 129 ChIntegrable1D<T>& integrand, ///< this is the integrand 130 const double a, ///< min limit for x domain 131 const double b, ///< min limit for x domain 132 const int order) ///< order of integration 133 { 134 ChQuadratureTables* mtables = 0; 135 std::vector<double>* lroots; 136 std::vector<double>* weight; 137 bool static_tables; 138 139 if ((unsigned int)order <= GetStaticTables()->Lroots.size()) { 140 mtables = GetStaticTables(); 141 lroots = &mtables->Lroots[order - 1]; 142 weight = &mtables->Weight[order - 1]; 143 static_tables = true; 144 } else { 145 mtables = new ChQuadratureTables(order, order); 146 mtables->PrintTables(); 147 lroots = &mtables->Lroots[0]; 148 weight = &mtables->Weight[0]; 149 static_tables = false; 150 } 151 152 double c1 = (b - a) / 2; 153 double c2 = (b + a) / 2; 154 155 result *= 0; // as result = 0, but works also for matrices. 156 T val; // temporary value for loop 157 158 for (unsigned int i = 0; i < lroots->size(); i++) { 159 integrand.Evaluate(val, (c1 * lroots->at(i) + c2)); 160 val *= weight->at(i); 161 result += val; 162 } 163 result *= c1; // result = c1 * sum; 164 165 if (!static_tables) 166 delete mtables; 167 } 168 169 /// Integrate the integrand T = f(x,y) over the 2D interval [xA, xB][yA, yB], 170 /// with desired order of quadrature. Best if integrand is polynomial. 171 /// For order in 1..10, precomputed polynomial coefficients are used for max speed. 172 template <class T> Integrate2D(T & result,ChIntegrable2D<T> & integrand,const double Xa,const double Xb,const double Ya,const double Yb,const int order)173 static void Integrate2D(T& result, ///< result is returned here 174 ChIntegrable2D<T>& integrand, ///< this is the integrand 175 const double Xa, ///< min limit for x domain 176 const double Xb, ///< min limit for x domain 177 const double Ya, ///< min limit for y domain 178 const double Yb, ///< min limit for y domain 179 const int order) ///< order of integration 180 { 181 ChQuadratureTables* mtables = 0; 182 std::vector<double>* lroots; 183 std::vector<double>* weight; 184 bool static_tables; 185 186 if ((unsigned int)order <= GetStaticTables()->Lroots.size()) { 187 mtables = GetStaticTables(); 188 lroots = &mtables->Lroots[order - 1]; 189 weight = &mtables->Weight[order - 1]; 190 static_tables = true; 191 } else { 192 mtables = new ChQuadratureTables(order, order); 193 mtables->PrintTables(); 194 lroots = &mtables->Lroots[0]; 195 weight = &mtables->Weight[0]; 196 static_tables = false; 197 } 198 199 double Xc1 = (Xb - Xa) / 2; 200 double Xc2 = (Xb + Xa) / 2; 201 double Yc1 = (Yb - Ya) / 2; 202 double Yc2 = (Yb + Ya) / 2; 203 204 result *= 0; // as result = 0, but works also for matrices. 205 T val; // temporary value for loop 206 207 for (unsigned int ix = 0; ix < lroots->size(); ix++) 208 for (unsigned int iy = 0; iy < lroots->size(); iy++) { 209 integrand.Evaluate(val, (Xc1 * lroots->at(ix) + Xc2), (Yc1 * lroots->at(iy) + Yc2)); 210 val *= (weight->at(ix) * weight->at(iy)); 211 result += val; 212 } 213 result *= (Xc1 * Yc1); 214 215 if (!static_tables) 216 delete mtables; 217 } 218 219 /// Integrate the integrand T = f(x,y,z) over the 3D interval [xA, xB][yA, yB][zA, zB], 220 /// with desired order of quadrature. Best if integrand is polynomial. 221 /// For order in 1..10, precomputed polynomial coefficients are used for max speed. 222 template <class T> Integrate3D(T & result,ChIntegrable3D<T> & integrand,const double Xa,const double Xb,const double Ya,const double Yb,const double Za,const double Zb,const int order)223 static void Integrate3D(T& result, ///< result is returned here 224 ChIntegrable3D<T>& integrand, ///< this is the integrand 225 const double Xa, ///< min limit for x domain 226 const double Xb, ///< min limit for x domain 227 const double Ya, ///< min limit for y domain 228 const double Yb, ///< min limit for y domain 229 const double Za, ///< min limit for z domain 230 const double Zb, ///< min limit for z domain 231 const int order) ///< order of integration 232 { 233 ChQuadratureTables* mtables = 0; 234 std::vector<double>* lroots; 235 std::vector<double>* weight; 236 bool static_tables; 237 238 if ((unsigned int)order <= GetStaticTables()->Lroots.size()) { 239 mtables = GetStaticTables(); 240 lroots = &mtables->Lroots[order - 1]; 241 weight = &mtables->Weight[order - 1]; 242 static_tables = true; 243 } else { 244 mtables = new ChQuadratureTables(order, order); 245 mtables->PrintTables(); 246 lroots = &mtables->Lroots[0]; 247 weight = &mtables->Weight[0]; 248 static_tables = false; 249 } 250 251 double Xc1 = (Xb - Xa) / 2; 252 double Xc2 = (Xb + Xa) / 2; 253 double Yc1 = (Yb - Ya) / 2; 254 double Yc2 = (Yb + Ya) / 2; 255 double Zc1 = (Zb - Za) / 2; 256 double Zc2 = (Zb + Za) / 2; 257 258 result *= 0; // as result = 0, but works also for matrices. 259 T val; // temporary value for loop 260 261 for (unsigned int ix = 0; ix < lroots->size(); ix++) 262 for (unsigned int iy = 0; iy < lroots->size(); iy++) 263 for (unsigned int iz = 0; iz < lroots->size(); iz++) { 264 integrand.Evaluate(val, (Xc1 * lroots->at(ix) + Xc2), (Yc1 * lroots->at(iy) + Yc2), 265 (Zc1 * lroots->at(iz) + Zc2)); 266 val *= (weight->at(ix) * weight->at(iy) * weight->at(iz)); 267 result += val; 268 } 269 result *= (Xc1 * Yc1 * Zc1); 270 271 if (!static_tables) 272 delete mtables; 273 } 274 275 276 /// Special case of 2D integration: integrate the integrand T = f(u,v) over a triangle, 277 /// with desired order of quadrature. Best if integrand is polynomial. 278 /// Two triangle coordinates are assumed to be 'area' coordinates u,v in [0...1]. The third is assumed 1-u-v. 279 /// For order in 1..5. Use precomputed polynomial coefficients for max speed. 280 template <class T> Integrate2Dtriangle(T & result,ChIntegrable2D<T> & integrand,const int order)281 static void Integrate2Dtriangle(T& result, ///< result is returned here 282 ChIntegrable2D<T>& integrand, ///< this is the integrand 283 const int order) ///< order of integration 284 { 285 if ((unsigned int)order > GetStaticTablesTriangle()->Weight.size()) 286 throw ChException("Too high order of quadrature for triangle. Use lower order."); 287 288 ChQuadratureTablesTriangle* mtables = GetStaticTablesTriangle(); 289 std::vector<double>* lrootsU = &mtables->LrootsU[order - 1]; 290 std::vector<double>* lrootsV = &mtables->LrootsV[order - 1]; 291 std::vector<double>* weight = &mtables->Weight[order - 1]; 292 293 result *= 0; // as result = 0, but works also for matrices. 294 T val; // temporary value for loop 295 296 for (unsigned int i = 0; i < weight->size(); i++) { 297 integrand.Evaluate(val, lrootsU->at(i), lrootsV->at(i)); 298 val *= weight->at(i) * 0.5; // the 1/2 coefficient is not in the table 299 result += val; 300 } 301 } 302 303 304 /// Special case of 3D integration: integrate the integrand T = f(u,v,w) over a tetrahedron, 305 /// with desired order of quadrature. Best if integrand is polynomial. 306 /// Three tetrahedron coordinates are assumed to be 'volume' coordinates u,v,w in [0...1]. The 4th is assumed 1-u-v-w. 307 /// For order in 1..5. Use precomputed polynomial coefficients for max speed. 308 template <class T> Integrate3Dtetrahedron(T & result,ChIntegrable3D<T> & integrand,const int order)309 static void Integrate3Dtetrahedron(T& result, ///< result is returned here 310 ChIntegrable3D<T>& integrand, ///< this is the integrand 311 const int order) ///< order of integration 312 { 313 if ((unsigned int)order > GetStaticTablesTetrahedron()->Weight.size()) 314 throw ChException("Too high order of quadrature for tetrahedron. Use lower order."); 315 316 ChQuadratureTablesTetrahedron* mtables = GetStaticTablesTetrahedron(); 317 std::vector<double>* lrootsU = &mtables->LrootsU[order - 1]; 318 std::vector<double>* lrootsV = &mtables->LrootsV[order - 1]; 319 std::vector<double>* lrootsW = &mtables->LrootsW[order - 1]; 320 std::vector<double>* weight = &mtables->Weight[order - 1]; 321 322 result *= 0; // as result = 0, but works also for matrices. 323 T val; // temporary value for loop 324 325 for (unsigned int i = 0; i < weight->size(); i++) { 326 integrand.Evaluate(val, lrootsU->at(i), lrootsV->at(i), lrootsW->at(i)); 327 val *= weight->at(i) *(1./6.); // the 1/6 coefficient is not in the table 328 result += val; 329 } 330 } 331 332 333 /// Access a statically-allocated set of tables, from 0 to a 10th order, 334 /// with precomputed tables. 335 static ChQuadratureTables* GetStaticTables(); 336 337 /// Access a statically-allocated set of tables for tetrahedron quadrature, 338 /// with 5 precomputed tables. 339 static ChQuadratureTablesTriangle* GetStaticTablesTriangle(); 340 341 /// Access a statically-allocated set of tables for tetrahedron quadrature, 342 /// with 5 precomputed tables. Use Dunavant theory. 343 static ChQuadratureTablesTetrahedron* GetStaticTablesTetrahedron(); 344 }; 345 346 } // end namespace chrono 347 348 #endif 349