60 for (
unsigned i = 0; i < d; ++i)
65 for (
unsigned i = 0; i < d; ++i)
76 for (
unsigned i = 0; i < d; ++i) {
77 target(i) = target(i) + scale*incr(i);
81 inline void multadd (
double& target,
const double scale,
const double incr)
88 target = target + scale*incr;
97 {0.166666666666667, 0.166666666666667, 0.666666666666667, 0.166666666666667},
98 {0.166666666666667, 0.666666666666667, 0.166666666666667, 0.166666666666667},
99 {0.666666666666667, 0.166666666666667, 0.166666666666667, 0.166666666666667},
100 {0.0, 0.0, 0.0, 0.0},
101 {0.0, 0.0, 0.0, 0.0},
102 {0.0, 0.0, 0.0, 0.0},
103 {0.0, 0.0, 0.0, 0.0},
104 {0.0, 0.0, 0.0, 0.0},
105 {0.0, 0.0, 0.0, 0.0},
106 {0.0, 0.0, 0.0, 0.0},
107 {0.0, 0.0, 0.0, 0.0},
108 {0.0, 0.0, 0.0, 0.0},
109 {0.0, 0.0, 0.0, 0.0},
110 {0.0, 0.0, 0.0, 0.0},
111 {0.0, 0.0, 0.0, 0.0},
117 {0.445948490915965, 0.445948490915965, 0.108103018168070, 0.111690794839005},
118 {0.445948490915965, 0.108103018168070, 0.445948490915965, 0.111690794839005},
119 {0.108103018168070, 0.445948490915965, 0.445948490915965, 0.111690794839005},
120 {0.091576213509771, 0.091576213509771, 0.816847572980458, 0.054975871827661},
121 {0.091576213509771, 0.816847572980458, 0.091576213509771, 0.054975871827661},
122 {0.816847572980458, 0.091576213509771, 0.091576213509771, 0.054975871827661},
123 {0.0, 0.0, 0.0, 0.0},
124 {0.0, 0.0, 0.0, 0.0},
125 {0.0, 0.0, 0.0, 0.0},
126 {0.0, 0.0, 0.0, 0.0},
127 {0.0, 0.0, 0.0, 0.0},
128 {0.0, 0.0, 0.0, 0.0},
129 {0.0, 0.0, 0.0, 0.0},
130 {0.0, 0.0, 0.0, 0.0},
131 {0.0, 0.0, 0.0, 0.0},
137 {0.333333333333333, 0.333333333333333, 0.333333333333333, 0.1125},
138 {0.470142064105115, 0.470142064105115, 0.059715871789770, 0.066197076394253},
139 {0.470142064105115, 0.059715871789770, 0.470142064105115, 0.066197076394253},
140 {0.059715871789770, 0.470142064105115, 0.470142064105115, 0.066197076394253},
141 {0.101286507323456, 0.101286507323456, 0.797426985353088, 0.062969590272414},
142 {0.101286507323456, 0.797426985353088, 0.101286507323456, 0.062969590272414},
143 {0.797426985353088, 0.101286507323456, 0.101286507323456, 0.062969590272414},
144 {0.0, 0.0, 0.0, 0.0},
145 {0.0, 0.0, 0.0, 0.0},
146 {0.0, 0.0, 0.0, 0.0},
147 {0.0, 0.0, 0.0, 0.0},
148 {0.0, 0.0, 0.0, 0.0},
149 {0.0, 0.0, 0.0, 0.0},
150 {0.0, 0.0, 0.0, 0.0},
151 {0.0, 0.0, 0.0, 0.0},
158 {0.333333333333333, 0.333333333333333, 0.333333333333333, 0.072157803838893},
159 {0.081414823414554, 0.459292588292722, 0.459292588292722, 0.047545817133642},
160 {0.459292588292722, 0.081414823414554, 0.459292588292722, 0.047545817133642},
161 {0.459292588292722, 0.459292588292722, 0.081414823414554, 0.047545817133642},
162 {0.898905543365937, 0.050547228317031, 0.050547228317031, 0.016229248811599},
163 {0.050547228317031, 0.898905543365937, 0.050547228317031, 0.016229248811599},
164 {0.050547228317031, 0.050547228317031, 0.898905543365937, 0.016229248811599},
165 {0.658861384496479, 0.170569307751760, 0.170569307751761, 0.051608685267359},
166 {0.170569307751760, 0.658861384496479, 0.170569307751761, 0.051608685267359},
167 {0.170569307751760, 0.170569307751761, 0.658861384496479, 0.051608685267359},
168 {0.008394777409957, 0.728492392955404, 0.263112829634639, 0.013615157087217},
169 {0.728492392955404, 0.008394777409957, 0.263112829634639, 0.013615157087217},
170 {0.728492392955404, 0.263112829634639, 0.008394777409957, 0.013615157087217},
171 {0.008394777409957, 0.263112829634639, 0.728492392955404, 0.013615157087217},
172 {0.263112829634639, 0.008394777409957, 0.728492392955404, 0.013615157087217},
173 {0.263112829634639, 0.728492392955404, 0.008394777409957, 0.013615157087217}
178 static const unsigned nbPts[4] = {3, 6, 7, 16};
180 template <
typename T,
typename I>
196 std::cout <<
"Unavailable Gauss order: min is 1, max is 3" << n << std::endl;
197 order = (n < 1) ? 1 : 3;
203 const Vect3 points[3] = { Trg.
s1(), Trg.
s2(), Trg.
s3() };
204 return triangle_integration(fc, points);
212 Vect3 crossprod = (points[1] - points[0])^(points[2] - points[0]);
213 double S = crossprod.
norm();
215 for (
unsigned i = 0; i <
nbPts[order];++i) {
216 Vect3 v(0.0, 0.0, 0.0);
217 for (
unsigned j = 0; j < 3; ++j) {
226 template <
typename T,
typename I>
237 inline double norm(
const double a) {
return fabs(a); }
241 const Vect3 points[3] = { Trg.
s1(), Trg.
s2(), Trg.
s3() };
242 T I0 = base::triangle_integration(fc, points);
243 return adaptive_integration(fc, points, I0, 0);
252 Vect3 newpoint0(0.0, 0.0, 0.0);
253 multadd(newpoint0, 0.5, points[0]);
254 multadd(newpoint0, 0.5, points[1]);
255 Vect3 newpoint1(0.0, 0.0, 0.0);
256 multadd(newpoint1, 0.5, points[1]);
257 multadd(newpoint1, 0.5, points[2]);
258 Vect3 newpoint2(0.0, 0.0, 0.0);
259 multadd(newpoint2, 0.5, points[2]);
260 multadd(newpoint2, 0.5, points[0]);
261 Vect3 points1[3] = {points[0], newpoint0, newpoint2};
262 Vect3 points2[3] = {points[1], newpoint1, newpoint0};
263 Vect3 points3[3] = {points[2], newpoint2, newpoint1};
264 Vect3 points4[3] = {newpoint0, newpoint1, newpoint2};
265 T I1 = base::triangle_integration(fc, points1);
266 T I2 = base::triangle_integration(fc, points2);
267 T I3 = base::triangle_integration(fc, points3);
268 T I4 = base::triangle_integration(fc, points4);
270 if ( norm(I0-sum) > tolerance*norm(I0) ) {
273 I1 = adaptive_integration(fc, points1, I1, n);
274 I2 = adaptive_integration(fc, points2, I2, n);
275 I3 = adaptive_integration(fc, points3, I3, n);
276 I4 = adaptive_integration(fc, points4, I4, n);
static const double cordBars[4][16][4]
static const unsigned nbPts[4]
void multadd(const double &d, const Vect3 &v)
AdaptiveIntegrator(double tol)
double norm(const Vect3 &a)
Vect3 operator()(int i) const
void setOrder(const unsigned n)
const Vertex & s2() const
virtual T integrate(const I &fc, const Triangle &Trg)
virtual T integrate(const I &fc, const Triangle &Trg)
double norm(const double a)
Vect3array< d > operator*(double x) const
Vect3 & operator()(int i)
void multadd(Vect3array< d > &target, const double scale, const Vect3array< d > &incr)
T adaptive_integration(const I &fc, const Vect3 *points, T I0, unsigned n)
T triangle_integration(const I &fc, const Vect3 points[3])
const Vertex & s1() const
const Vertex & s3() const