mbsim 4.0.0 MBSim Kernel

gaussian_quadratur.h
1 /* Copyright (C) 2004-2015 MBSim Development Team
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16 *
17 * Contact: martin.o.foerg@googlemail.com
18 */
19
20 #ifndef GAUSSIAN_QUADRATUR_H_
21 #define GAUSSIAN_QUADRATUR_H_
22
23 #include <fmatvec/fmatvec.h>
24 #include <mbsim/functions/function.h>
25
26 #include <mbsim/mbsim_event.h>
27
28 namespace MBSim {
29
30 static double transformVariablePosition(double x, double l, double u) {
31 return (u - l) / 2. * x + (l + u) / 2.;
32 }
33
39 static double pW[15][2] = {{0., 2.}, {-0.57735026919, 1.}, {0.57735026919, 1}, {-0.774596669241, 5. / 9.}, {0., 8. / 9}, {0.774596669241, 5. / 9.}, {-0.861136311594053, 0.347854845137454}, {-0.339981043584856, 0.652145154862546}, {0.339981043584856, 0.652145154862546}, {0.861136311594053, 0.347854845137454}, {-0.906179845938664, 0.236926885056189}, {-0.538469310105683, 0.478628670499366}, {0, 0.568888888888889}, {0.538469310105683, 0.478628670499366}, {0.906179845938664, 0.236926885056189}};
40
41 template <class Ret>
42 class GaussLegendreQuadrature {
43 public:
44 GaussLegendreQuadrature(fmatvec::Function<Ret(double)> * function, int dimension, int subIntegrals = 1, int Nb = 3) :
45 function(function), dim(dimension), subIntegrals(subIntegrals), Nb(Nb) {
46 }
47 ~GaussLegendreQuadrature() {
48 }
49
56 const Ret integrate(double l, double u) {
57 double currLow = l;
58 double step = (u - l) / subIntegrals;
59
60 Ret result(dim, fmatvec::INIT, 0.);
61
62 for (size_t i = 0; i < subIntegrals; i++) {
63 result += integrateSingle(currLow, currLow + step);
64 currLow += step;
65 }
66
67 return result;
68 }
69
70 /* GETTER/SETTER*/
71 void setSubIntegrals(int subInts_) {
72 subIntegrals = subInts_;
73 }
74
78 fmatvec::Function<Ret(double)> * function;
79
83 int dim;
84
88 size_t subIntegrals;
89
93 size_t Nb;
94
95 protected:
96 const Ret integrateSingle(double l, double u) {
97 if (Nb > 5)
98 throw std::runtime_error("It is not possible to use more than 5 Gauss-Points atm!");
99
100 int posInpW = 0;
101
102 for (unsigned int i = 1; i < Nb; i++) {
103 posInpW += i;
104 }
105
106 Ret result(dim, fmatvec::INIT, 0.);
107 for (unsigned int i = 0; i < Nb; i++) {
108 int index = i + posInpW;
109 double p = transformVariablePosition(pW[index][0], l, u);
110 double w = pW[index][1];
111 result += w * (*function)(p);
112 }
113
114 return (u - l) / 2. * result;
115 }
116 };
117
118 namespace GaussTschebyschowQuadrature {
128 template <class Ret>
129 const Ret integrate(fmatvec::Function<Ret(double)> * function, double l, double u, size_t Nb = 3) {
130 Ret result((*function)(l).size(), fmatvec::INIT, 0.);
131 for (size_t i = 0; i < Nb; i++) {
132 double x = cos((2. * i - 1) / (2. * Nb) * M_PI);
133 double p = transformVariablePosition(x, l, u);
134 result += sqrt(1 - x * x) * (*function)(p);
135 }
136
137 return (u - l) / 2. * M_PI / Nb * result;
138 }
139 }
140 }
141
142 #endif /* GAUSSIAN_QUADRATUR_H_ */
MBSim::GaussLegendreQuadrature::dim
int dim
dimension of the result
Definition: gaussian_quadratur.h:83
MBSim::GaussLegendreQuadrature
Definition: gaussian_quadratur.h:42
MBSim::GaussLegendreQuadrature::Nb
size_t Nb
number of gauss points that should be used
Definition: gaussian_quadratur.h:93
MBSim::GaussLegendreQuadrature::integrate
const Ret integrate(double l, double u)
function to integrate a MBSim-Function numerically using the gaussian quadrature
Definition: gaussian_quadratur.h:56
MBSim::GaussLegendreQuadrature::subIntegrals
size_t subIntegrals
number of integrals the original one should be split to
Definition: gaussian_quadratur.h:88
MBSim
namespace MBSim
Definition: bilateral_constraint.cc:30