GaussQuadratureFormula.hxx 12.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
///
////// \file
///
///
/// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Mon, 29 Sep 2014 10:25:42 +0200
/// Copyright (c) Inria. All rights reserved.
///
/// \ingroup FiniteElementGroup
/// \addtogroup FiniteElementGroup
/// \{
11

12 13
#ifndef MOREFEM_x_FINITE_ELEMENT_x_QUADRATURE_RULES_x_GAUSS_QUADRATURE_FORMULA_HXX_
# define MOREFEM_x_FINITE_ELEMENT_x_QUADRATURE_RULES_x_GAUSS_QUADRATURE_FORMULA_HXX_
14 15


16 17
#include <iostream> // \todo DEV TMP #907

18
namespace MoReFEM
19
{
20 21


22 23
    namespace QuadratureNS
    {
24 25


26
        template<GaussFormula QuadratureFormulaT>
27 28 29 30 31 32 33 34 35 36
        void ComputeGaussFormulas(const unsigned int Ngauss_points,
                                  std::vector<double>& points,
                                  std::vector<double>& weights)
        {
            // For details about that function, take a look at
            //
            //     P. N. swarztrauber, Computing the points and weights for
            //     Gauss-Legendre quadrature, SIAM J. Sci. Comput.,
            //     24(2002) pp. 945-954.
            //
37

38
            assert(Ngauss_points > 0u && "Added in MoReFEM; case 0 caused crashes...");
39

40
            const double eps = std::numeric_limits<double>::epsilon();
41

42 43
            const double pi = static_cast<double>(M_PI);
            const double half_pi = static_cast<double>(M_PI_2);
44

45
            // For Gauss-Lobatto, we compute internal nodes.
46
            assert(!(QuadratureFormulaT == GaussFormula::Gauss_Lobatto && Ngauss_points < 2u));
47

48
            const unsigned int Nquadrature_point
49
                = (QuadratureFormulaT == GaussFormula::Gauss_Lobatto ? Ngauss_points - 1u : Ngauss_points);
50

51 52 53 54
            // points(i) will be equal to cos(theta(i))
            std::vector<double> theta(Ngauss_points);
            points.resize(Ngauss_points);
            weights.resize(Ngauss_points);
55

56 57 58 59 60
            // coefficient of Fourier transform of Legendre polynom and derivative
            // we need the second derivative for Gauss-Lobatto points
            std::vector<double> coef(Nquadrature_point, 0.);
            std::vector<double> deriv_coef(Nquadrature_point, 0.);
            std::vector<double> second_deriv_coef(Nquadrature_point, 0.);
61

62 63 64 65
            // exact expressions are used for order 1, 2, 3
            switch (Nquadrature_point)
            {
                case 1:
66
                    Internal::GaussQuadratureNS::ComputeExactFormula<QuadratureFormulaT, 1u>::Perform(points, weights);
67 68
                    return;
                case 2:
69
                    Internal::GaussQuadratureNS::ComputeExactFormula<QuadratureFormulaT, 2u>::Perform(points, weights);
70 71
                    return;
                case 3:
72
                    Internal::GaussQuadratureNS::ComputeExactFormula<QuadratureFormulaT, 3u>::Perform(points, weights);
73 74 75 76
                    return;
                default:
                    break;
            }
77

78 79 80 81
            // for order greater than 3, numerical computation
            const unsigned int Nparity_points = Nquadrature_point % 2u;
            const unsigned int Nhalf_points = Ngauss_points / 2u;
            const unsigned int Nhalf = (Nquadrature_point + 1u) / 2u;
82

83 84
            double zero = 0.0;
            double cz = zero;
85

86
            Internal::GaussQuadratureNS::ComputeFourierCoef(cz, coef, deriv_coef, second_deriv_coef);
87

88
            double dtheta = zero;
89
            if (QuadratureFormulaT == GaussFormula::Gauss_Lobatto && (Nparity_points==1))
90
                dtheta = half_pi / (Nhalf - 1);
91
            else
92
                dtheta = half_pi / Nhalf;
93

94 95
            const double dthalf = 0.5 * dtheta;
            const double cmax = 0.2 * dtheta;
96

97 98 99 100
            // the zeros of Legendre polynom
            double zprev = zero, zlast = zero, zhold = zero;
            double pb = zero, dpb = zero, ddpb = zero, dcor = zero, sgnd = zero;
            unsigned int nix;
101

102 103 104
            //
            //     estimate first point next to theta = pi/2
            //
105
            if (Nparity_points != 0)
106 107 108 109
            {
                // Nquadrature_point = 2 Nhalf-1
                // if odd the first zero is at the middle pi/2
                // and the following pi/2-pi/(2*Nhalf)
110
                if (QuadratureFormulaT == GaussFormula::Gauss_Lobatto)
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
                {
                    zero = half_pi - dthalf;
                    zprev = half_pi + dthalf;
                    nix = Nhalf - 1;
                }
                else
                {
                    zero = half_pi - dtheta;
                    zprev = half_pi;
                    nix = Nhalf - 1; // index of the point
                }
            }
            else
            {
                // if even, no zero on the middle, the first zero is on pi/2-pi/(4*Nhalf)
126
                if (QuadratureFormulaT == GaussFormula::Gauss_Lobatto)
127 128 129
                {
                    zero = half_pi - dtheta;
                    zprev = half_pi;
130
                    nix = Nhalf - 1; // index of the point
131 132 133
                }
                else
                {
134 135
                    zero = half_pi - dthalf;
                    zprev = half_pi + dthalf;
136 137 138
                    nix = Nhalf;
                }
            }
139

140
            bool each_point_not_computed = true;
141

142 143
            while (each_point_not_computed)
            {
144
                int Niteration = 0;
145
                bool test_convergence_newton = true;
146
                double residu(1.0);
147

148
                while (test_convergence_newton && Niteration < 100)
149
                {
150
                    ++Niteration;
151
                    zlast = zero;
152

153 154 155
                    //
                    //     newton iterations
                    //
156
                    Internal::GaussQuadratureNS::ComputeLegendrePolAndDerivative(Nquadrature_point, zero, cz, coef,deriv_coef,
157
                                                             second_deriv_coef,pb, dpb, ddpb);
158

159
                    if (QuadratureFormulaT == GaussFormula::Gauss_Lobatto)
160
                        dcor = dpb / ddpb;
161
                    else
162
                        dcor = pb / dpb;
163

164
                    sgnd = 1.0;
165
                    if (!NumericNS::IsZero(dcor))
166
                        sgnd = dcor/std::fabs(dcor);
167

168
                    // we don't move the point further than 0.2*delta_theta
169 170 171
                    double tmp = std::fabs(dcor);
                    dcor = sgnd * std::min(tmp, cmax);
                    zero = zero - dcor;
172
                    double residu_prec = residu;
173
                    residu = std::fabs(zero - zlast);
174
                    // we check if the stopping criteria are reached
175 176 177
                    if ((std::fabs(zero - zlast) <= eps * std::fabs(zero))
                        || ((Niteration > 5) && (residu_prec < residu)))
                        test_convergence_newton = false;
178
                }
179

180 181
                assert(nix - 1 < theta.size());
                theta[nix - 1] = zero;
182 183 184 185 186
                zhold = zero;
                //      weights(nix) = (Nquadrature_point+Nquadrature_point+1)/(dpb*dpb)
                //
                //     yakimiw's formula permits using old pb and dpb
                //
187
                if (QuadratureFormulaT == GaussFormula::Gauss_Lobatto)
188
                {
189
                    double tmp = NumericNS::Square(pb - dcor * dpb + 0.5 * dcor * dcor * ddpb);
190
                    assert(!NumericNS::IsZero(tmp));
191 192
                    assert(nix - 1 < weights.size());
                    weights[nix - 1] = 1.0 / tmp;
193 194 195
                }
                else
                {
196
                    assert(!NumericNS::IsZero(std::sin(zlast)));
197
                    double tmp = NumericNS::Square(dpb + pb * std::cos(zlast) / std::sin(zlast));
198
                    assert(!NumericNS::IsZero(tmp));
199 200
                    assert(nix - 1 < weights.size());
                    weights[nix - 1] = static_cast<double>(2u * Nquadrature_point + 1u) / tmp;
201
                }
202

203
                --nix;
204

205 206
                if (nix == 0)
                    each_point_not_computed = false;
207 208
                else if (nix <= (Nhalf - 1))
                    zero += zero - zprev;
209

210 211
                zprev = zhold;
            }
212

213 214 215
            //
            //     extend points and weights via symmetries
            //
216
            if ((QuadratureFormulaT == GaussFormula::Gauss) && (Nparity_points != 0))
217
            {
218
                assert(Nhalf - 1 < theta.size());
219
                theta[Nhalf - 1] = half_pi;
220 221 222
                Internal::GaussQuadratureNS::ComputeLegendrePolAndDerivative(Nquadrature_point,
                                                                             half_pi, cz, coef,deriv_coef,
                                                                             second_deriv_coef, pb, dpb, ddpb);
223

224
                assert(!NumericNS::IsZero(dpb));
225
                assert(Nhalf - 1 < weights.size());
226
                weights[Nhalf - 1] = static_cast<double>(2u * Nquadrature_point + 1u) / NumericNS::Square(dpb);
227
            }
228

229
            if ((QuadratureFormulaT == GaussFormula::Gauss_Lobatto) && (Nparity_points == 0))
230
            {
231
                assert(Nhalf - 1 < theta.size());
232
                theta[Nhalf - 1] = half_pi;
233
                Internal::GaussQuadratureNS::ComputeLegendrePolAndDerivative(Nquadrature_point, half_pi, cz, coef,deriv_coef,
234
                                                         second_deriv_coef,pb, dpb, ddpb);
235
                assert(!NumericNS::IsZero(pb));
236
                assert(Nhalf - 1 < weights.size());
237
                weights[Nhalf - 1] = 1. / NumericNS::Square(pb);
238
            }
239

240
            // DISP(Nhalf); DISP(theta);DISP(weights);
241
            if (QuadratureFormulaT == GaussFormula::Gauss_Lobatto)
242
            {
243 244
                assert(theta.size() >= Nhalf);
                assert(weights.size() >= Nhalf);
245

246
                for (auto i = static_cast<int>(Nhalf) - 1; i >= 0; --i)
247
                {
248
                    theta[static_cast<std::size_t>(i + 1)] = theta[static_cast<std::size_t>(i)];
249
                    weights[static_cast<std::size_t>(i + 1)] = weights[static_cast<std::size_t>(i)];
250
                }
251 252


253
                theta[0] = 0.0;
254 255 256
                Internal::GaussQuadratureNS::ComputeLegendrePolAndDerivative(Nquadrature_point,
                                                                             theta[0], cz, coef,deriv_coef,
                                                                             second_deriv_coef, pb, dpb, ddpb);
257

258
                // DISP(cz); DISP(pb);
259
                weights[0] = 1.0 / (NumericNS::Square(pb));
260
            }
261

262 263 264
            // DISP(theta);
            for (unsigned int i = 0; i < Nhalf_points; i++)
            {
265 266
                assert(Ngauss_points-i - 1 < weights.size());
                assert(Ngauss_points-i - 1 < theta.size());
267 268
                weights[Ngauss_points-i - 1] = weights[i];
                theta[Ngauss_points-i - 1] = pi - theta[i];
269
            }
270

271
            // DISP(theta); DISP(weights);
272
            double sum = 0.;
273

274 275
            assert(Ngauss_points <= weights.size());
            assert(Ngauss_points <= theta.size());
276 277
            for (unsigned int i = 0; i < Ngauss_points; i++)
                sum += weights[i];
278

279
            assert(!NumericNS::IsZero(sum));
280

281 282
            for (unsigned int i = 0; i < Ngauss_points; i++)
            {
283 284
                weights[i] = 2. *  weights[i] / sum;
                points[i] = std::cos(theta[i]);
285
            }
286

287
            //reverse the arrays
288
            const auto Nweight = weights.size();
289

290
            for (auto i = 0ul; i < Nweight / 2 ; i++)
291 292
            {
                double tmp = weights[i];
293 294
                assert(weights.size() >= i + 1);
                std::size_t index = weights.size() - i - 1;
295

296 297
                weights[i] = weights[index];
                weights[index] = tmp;
298
            }
299

300
            const auto Npoints = points.size();
301

302
            for (auto i = 0ul; i < Npoints / 2; i++)
303 304
            {
                double tmp = points[i];
305
                assert(points.size() >= i + 1);
306
                std::size_t index = Nweight - i - 1;
307

308 309
                points[i] = points[index];
                points[index] = tmp;
310
            }
311

312
        }
313 314


315
    } // namespace QuadratureNS
316 317


318
} // namespace MoReFEM
319 320


321 322 323
/// @} // addtogroup FiniteElementGroup


324
#endif // MOREFEM_x_FINITE_ELEMENT_x_QUADRATURE_RULES_x_GAUSS_QUADRATURE_FORMULA_HXX_