Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
polynomial.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
18#include <deal.II/base/point.h>
23
24#include <algorithm>
25#include <cmath>
26#include <limits>
27#include <memory>
28
30
31
32
33// have a lock that guarantees that at most one thread is changing and
34// accessing the @p{coefficients} arrays of classes implementing
35// polynomials with tables. make this lock local to this file.
36//
37// having only one lock for all of these classes is probably not going
38// to be a problem since we only need it on very rare occasions. if
39// someone finds this is a bottleneck, feel free to replace it by a
40// more fine-grained solution
41namespace
42{
43 std::mutex coefficients_lock;
44}
45
46
47
48namespace Polynomials
49{
50 // -------------------- class Polynomial ---------------- //
51
52
53 template <typename number>
54 Polynomial<number>::Polynomial(const std::vector<number> &a)
55 : coefficients(a)
56 , in_lagrange_product_form(false)
57 , lagrange_weight(1.)
58 {}
59
60
61
62 template <typename number>
63 Polynomial<number>::Polynomial(const unsigned int n)
64 : coefficients(n + 1, 0.)
65 , in_lagrange_product_form(false)
66 , lagrange_weight(1.)
67 {}
68
69
70
71 template <typename number>
72 Polynomial<number>::Polynomial(const std::vector<Point<1>> &supp,
73 const unsigned int center)
74 : in_lagrange_product_form(true)
75 {
76 Assert(supp.size() > 0, ExcEmptyObject());
77 AssertIndexRange(center, supp.size());
78
79 lagrange_support_points.reserve(supp.size() - 1);
80 number tmp_lagrange_weight = 1.;
81 for (unsigned int i = 0; i < supp.size(); ++i)
82 if (i != center)
83 {
84 lagrange_support_points.push_back(supp[i](0));
85 tmp_lagrange_weight *= supp[center](0) - supp[i](0);
86 }
87
88 // check for underflow and overflow
89 Assert(std::fabs(tmp_lagrange_weight) > std::numeric_limits<number>::min(),
90 ExcMessage("Underflow in computation of Lagrange denominator."));
91 Assert(std::fabs(tmp_lagrange_weight) < std::numeric_limits<number>::max(),
92 ExcMessage("Overflow in computation of Lagrange denominator."));
93
94 lagrange_weight = static_cast<number>(1.) / tmp_lagrange_weight;
95 }
96
97
98
99 template <typename number>
100 void
101 Polynomial<number>::value(const number x, std::vector<number> &values) const
102 {
103 Assert(values.size() > 0, ExcZero());
104
105 value(x, values.size() - 1, values.data());
106 }
107
108
109
110 template <typename number>
111 void
113 {
114 // should only be called when the product form is active
115 Assert(in_lagrange_product_form == true, ExcInternalError());
116 Assert(coefficients.size() == 0, ExcInternalError());
117
118 // compute coefficients by expanding the product (x-x_i) term by term
119 coefficients.resize(lagrange_support_points.size() + 1);
120 if (lagrange_support_points.size() == 0)
121 coefficients[0] = 1.;
122 else
123 {
124 coefficients[0] = -lagrange_support_points[0];
125 coefficients[1] = 1.;
126 for (unsigned int i = 1; i < lagrange_support_points.size(); ++i)
127 {
128 coefficients[i + 1] = 1.;
129 for (unsigned int j = i; j > 0; --j)
130 coefficients[j] = (-lagrange_support_points[i] * coefficients[j] +
131 coefficients[j - 1]);
132 coefficients[0] *= -lagrange_support_points[i];
133 }
134 }
135 for (unsigned int i = 0; i < lagrange_support_points.size() + 1; ++i)
136 coefficients[i] *= lagrange_weight;
137
138 // delete the product form data
139 std::vector<number> new_points;
140 lagrange_support_points.swap(new_points);
141 in_lagrange_product_form = false;
142 lagrange_weight = 1.;
143 }
144
145
146
147 template <typename number>
148 void
149 Polynomial<number>::scale(std::vector<number> &coefficients,
150 const number factor)
151 {
152 number f = 1.;
153 for (typename std::vector<number>::iterator c = coefficients.begin();
154 c != coefficients.end();
155 ++c)
156 {
157 *c *= f;
158 f *= factor;
159 }
160 }
161
162
163
164 template <typename number>
165 void
166 Polynomial<number>::scale(const number factor)
167 {
168 // to scale (x-x_0)*(x-x_1)*...*(x-x_n), scale
169 // support points by 1./factor and the weight
170 // likewise
171 if (in_lagrange_product_form == true)
172 {
173 number inv_fact = number(1.) / factor;
174 number accumulated_fact = 1.;
175 for (unsigned int i = 0; i < lagrange_support_points.size(); ++i)
176 {
177 lagrange_support_points[i] *= inv_fact;
178 accumulated_fact *= factor;
179 }
180 lagrange_weight *= accumulated_fact;
182 // otherwise, use the function above
183 else
184 scale(coefficients, factor);
185 }
186
187
188
189 template <typename number>
190 void
191 Polynomial<number>::multiply(std::vector<number> &coefficients,
192 const number factor)
193 {
194 for (typename std::vector<number>::iterator c = coefficients.begin();
195 c != coefficients.end();
196 ++c)
197 *c *= factor;
198 }
199
201
202 template <typename number>
205 {
206 if (in_lagrange_product_form == true)
207 lagrange_weight *= s;
208 else
209 {
210 for (typename std::vector<number>::iterator c = coefficients.begin();
211 c != coefficients.end();
212 ++c)
213 *c *= s;
214 }
215 return *this;
216 }
217
218
220 template <typename number>
223 {
224 // if we are in Lagrange form, just append the
225 // new points
226 if (in_lagrange_product_form == true && p.in_lagrange_product_form == true)
227 {
228 lagrange_weight *= p.lagrange_weight;
229 lagrange_support_points.insert(lagrange_support_points.end(),
230 p.lagrange_support_points.begin(),
232 }
233
234 // cannot retain product form, recompute...
235 else if (in_lagrange_product_form == true)
236 transform_into_standard_form();
238 // need to transform p into standard form as
239 // well if necessary. copy the polynomial to
240 // do this
241 std::unique_ptr<Polynomial<number>> q_data;
242 const Polynomial<number> * q = nullptr;
244 {
245 q_data = std::make_unique<Polynomial<number>>(p);
247 q = q_data.get();
248 }
249 else
250 q = &p;
251
252 // Degree of the product
253 unsigned int new_degree = this->degree() + q->degree();
254
255 std::vector<number> new_coefficients(new_degree + 1, 0.);
256
257 for (unsigned int i = 0; i < q->coefficients.size(); ++i)
258 for (unsigned int j = 0; j < this->coefficients.size(); ++j)
259 new_coefficients[i + j] += this->coefficients[j] * q->coefficients[i];
260 this->coefficients = new_coefficients;
261
262 return *this;
263 }
265
266
267 template <typename number>
270 {
271 // Lagrange product form cannot reasonably be
272 // retained after polynomial addition. we
273 // could in theory check if either this
274 // polynomial or the other is a zero
275 // polynomial and retain it, but we actually
276 // currently (r23974) assume that the addition
277 // of a zero polynomial changes the state and
278 // tests equivalence.
279 if (in_lagrange_product_form == true)
280 transform_into_standard_form();
281
282 // need to transform p into standard form as
283 // well if necessary. copy the polynomial to
284 // do this
285 std::unique_ptr<Polynomial<number>> q_data;
286 const Polynomial<number> * q = nullptr;
287 if (p.in_lagrange_product_form == true)
288 {
289 q_data = std::make_unique<Polynomial<number>>(p);
291 q = q_data.get();
293 else
294 q = &p;
295
296 // if necessary expand the number
297 // of coefficients we store
298 if (q->coefficients.size() > coefficients.size())
299 coefficients.resize(q->coefficients.size(), 0.);
300
301 for (unsigned int i = 0; i < q->coefficients.size(); ++i)
302 coefficients[i] += q->coefficients[i];
303
304 return *this;
305 }
306
307
308
309 template <typename number>
310 Polynomial<number> &
312 {
313 // Lagrange product form cannot reasonably be
314 // retained after polynomial addition
315 if (in_lagrange_product_form == true)
316 transform_into_standard_form();
317
318 // need to transform p into standard form as
319 // well if necessary. copy the polynomial to
320 // do this
321 std::unique_ptr<Polynomial<number>> q_data;
322 const Polynomial<number> * q = nullptr;
323 if (p.in_lagrange_product_form == true)
324 {
325 q_data = std::make_unique<Polynomial<number>>(p);
327 q = q_data.get();
328 }
329 else
330 q = &p;
331
332 // if necessary expand the number
333 // of coefficients we store
334 if (q->coefficients.size() > coefficients.size())
335 coefficients.resize(q->coefficients.size(), 0.);
336
337 for (unsigned int i = 0; i < q->coefficients.size(); ++i)
338 coefficients[i] -= q->coefficients[i];
339
340 return *this;
341 }
342
343
344
345 template <typename number>
346 bool
348 {
349 // need to distinguish a few cases based on
350 // whether we are in product form or not. two
351 // polynomials can still be the same when they
352 // are on different forms, but the expansion
353 // is the same
354 if (in_lagrange_product_form == true && p.in_lagrange_product_form == true)
355 return ((lagrange_weight == p.lagrange_weight) &&
356 (lagrange_support_points == p.lagrange_support_points));
357 else if (in_lagrange_product_form == true)
358 {
359 Polynomial<number> q = *this;
361 return (q.coefficients == p.coefficients);
362 }
363 else if (p.in_lagrange_product_form == true)
364 {
365 Polynomial<number> q = p;
367 return (q.coefficients == coefficients);
368 }
369 else
370 return (p.coefficients == coefficients);
371 }
372
373
374
375 template <typename number>
376 template <typename number2>
377 void
378 Polynomial<number>::shift(std::vector<number> &coefficients,
379 const number2 offset)
380 {
381 // too many coefficients cause overflow in
382 // the binomial coefficient used below
383 Assert(coefficients.size() < 31, ExcNotImplemented());
384
385 // Copy coefficients to a vector of
386 // accuracy given by the argument
387 std::vector<number2> new_coefficients(coefficients.begin(),
388 coefficients.end());
389
390 // Traverse all coefficients from
391 // c_1. c_0 will be modified by
392 // higher degrees, only.
393 for (unsigned int d = 1; d < new_coefficients.size(); ++d)
394 {
395 const unsigned int n = d;
396 // Binomial coefficients are
397 // needed for the
398 // computation. The rightmost
399 // value is unity.
400 unsigned int binomial_coefficient = 1;
401
402 // Powers of the offset will be
403 // needed and computed
404 // successively.
405 number2 offset_power = offset;
406
407 // Compute (x+offset)^d
408 // and modify all values c_k
409 // with k<d.
410 // The coefficient in front of
411 // x^d is not modified in this step.
412 for (unsigned int k = 0; k < d; ++k)
413 {
414 // Recursion from Bronstein
415 // Make sure no remainders
416 // occur in integer
417 // division.
418 binomial_coefficient = (binomial_coefficient * (n - k)) / (k + 1);
419
420 new_coefficients[d - k - 1] +=
421 new_coefficients[d] * binomial_coefficient * offset_power;
422 offset_power *= offset;
423 }
424 // The binomial coefficient
425 // should have gone through a
426 // whole row of Pascal's
427 // triangle.
428 Assert(binomial_coefficient == 1, ExcInternalError());
429 }
430
431 // copy new elements to old vector
432 coefficients.assign(new_coefficients.begin(), new_coefficients.end());
433 }
434
435
436
437 template <typename number>
438 template <typename number2>
439 void
440 Polynomial<number>::shift(const number2 offset)
441 {
442 // shift is simple for a polynomial in product
443 // form, (x-x_0)*(x-x_1)*...*(x-x_n). just add
444 // offset to all shifts
445 if (in_lagrange_product_form == true)
446 {
447 for (unsigned int i = 0; i < lagrange_support_points.size(); ++i)
448 lagrange_support_points[i] -= offset;
449 }
450 else
451 // do the shift in any case
452 shift(coefficients, offset);
453 }
454
455
456
457 template <typename number>
460 {
461 // no simple form possible for Lagrange
462 // polynomial on product form
463 if (degree() == 0)
464 return Monomial<number>(0, 0.);
465
466 std::unique_ptr<Polynomial<number>> q_data;
467 const Polynomial<number> * q = nullptr;
468 if (in_lagrange_product_form == true)
469 {
470 q_data = std::make_unique<Polynomial<number>>(*this);
472 q = q_data.get();
473 }
474 else
475 q = this;
476
477 std::vector<number> newcoefficients(q->coefficients.size() - 1);
478 for (unsigned int i = 1; i < q->coefficients.size(); ++i)
479 newcoefficients[i - 1] = number(i) * q->coefficients[i];
480
481 return Polynomial<number>(newcoefficients);
482 }
483
484
485
486 template <typename number>
489 {
490 // no simple form possible for Lagrange
491 // polynomial on product form
492 std::unique_ptr<Polynomial<number>> q_data;
493 const Polynomial<number> * q = nullptr;
494 if (in_lagrange_product_form == true)
495 {
496 q_data = std::make_unique<Polynomial<number>>(*this);
498 q = q_data.get();
499 }
500 else
501 q = this;
502
503 std::vector<number> newcoefficients(q->coefficients.size() + 1);
504 newcoefficients[0] = 0.;
505 for (unsigned int i = 0; i < q->coefficients.size(); ++i)
506 newcoefficients[i + 1] = q->coefficients[i] / number(i + 1.);
507
508 return Polynomial<number>(newcoefficients);
509 }
510
511
512
513 template <typename number>
514 void
515 Polynomial<number>::print(std::ostream &out) const
516 {
517 if (in_lagrange_product_form == true)
518 {
519 out << lagrange_weight;
520 for (unsigned int i = 0; i < lagrange_support_points.size(); ++i)
521 out << " (x-" << lagrange_support_points[i] << ")";
522 out << std::endl;
523 }
524 else
525 for (int i = degree(); i >= 0; --i)
526 {
527 out << coefficients[i] << " x^" << i << std::endl;
528 }
529 }
530
531
532 template <typename number>
533 std::size_t
535 {
536 return (MemoryConsumption::memory_consumption(coefficients) +
537 MemoryConsumption::memory_consumption(in_lagrange_product_form) +
538 MemoryConsumption::memory_consumption(lagrange_support_points) +
540 }
541
542
543
544 // ------------------ class Monomial -------------------------- //
545
546 template <typename number>
547 std::vector<number>
548 Monomial<number>::make_vector(unsigned int n, double coefficient)
549 {
550 std::vector<number> result(n + 1, 0.);
551 result[n] = coefficient;
552 return result;
553 }
554
555
556
557 template <typename number>
558 Monomial<number>::Monomial(unsigned int n, double coefficient)
559 : Polynomial<number>(make_vector(n, coefficient))
560 {}
561
562
563
564 template <typename number>
565 std::vector<Polynomial<number>>
567 {
568 std::vector<Polynomial<number>> v;
569 v.reserve(degree + 1);
570 for (unsigned int i = 0; i <= degree; ++i)
571 v.push_back(Monomial<number>(i));
572 return v;
573 }
574
575
576
577 // ------------------ class LagrangeEquidistant --------------- //
578
579 namespace internal
580 {
581 namespace LagrangeEquidistantImplementation
582 {
583 std::vector<Point<1>>
585 {
586 std::vector<Point<1>> points(n + 1);
587 const double one_over_n = 1. / n;
588 for (unsigned int k = 0; k <= n; ++k)
589 points[k](0) = static_cast<double>(k) * one_over_n;
590 return points;
591 }
592 } // namespace LagrangeEquidistantImplementation
593 } // namespace internal
594
595
596
598 const unsigned int support_point)
599 : Polynomial<double>(internal::LagrangeEquidistantImplementation::
600 generate_equidistant_unit_points(n),
601 support_point)
602 {
603 Assert(coefficients.size() == 0, ExcInternalError());
604
605 // For polynomial order up to 3, we have precomputed weights. Use these
606 // weights instead of the product form
607 if (n <= 3)
608 {
609 this->in_lagrange_product_form = false;
610 this->lagrange_weight = 1.;
611 std::vector<double> new_support_points;
612 this->lagrange_support_points.swap(new_support_points);
613 this->coefficients.resize(n + 1);
614 compute_coefficients(n, support_point, this->coefficients);
615 }
616 }
617
618
619
620 void
622 const unsigned int support_point,
623 std::vector<double> &a)
624 {
625 AssertIndexRange(support_point, n + 1);
626
627 unsigned int n_functions = n + 1;
628 AssertIndexRange(support_point, n_functions);
629 double const *x = nullptr;
630
631 switch (n)
632 {
633 case 1:
634 {
635 static const double x1[4] = {1.0, -1.0, 0.0, 1.0};
636 x = &x1[0];
637 break;
638 }
639 case 2:
640 {
641 static const double x2[9] = {
642 1.0, -3.0, 2.0, 0.0, 4.0, -4.0, 0.0, -1.0, 2.0};
643 x = &x2[0];
644 break;
645 }
646 case 3:
647 {
648 static const double x3[16] = {1.0,
649 -11.0 / 2.0,
650 9.0,
651 -9.0 / 2.0,
652 0.0,
653 9.0,
654 -45.0 / 2.0,
655 27.0 / 2.0,
656 0.0,
657 -9.0 / 2.0,
658 18.0,
659 -27.0 / 2.0,
660 0.0,
661 1.0,
662 -9.0 / 2.0,
663 9.0 / 2.0};
664 x = &x3[0];
665 break;
666 }
667 default:
668 Assert(false, ExcInternalError())
669 }
670
671 Assert(x != nullptr, ExcInternalError());
672 for (unsigned int i = 0; i < n_functions; ++i)
673 a[i] = x[support_point * n_functions + i];
674 }
675
676
677
678 std::vector<Polynomial<double>>
680 {
681 if (degree == 0)
682 // create constant polynomial
683 return std::vector<Polynomial<double>>(
684 1, Polynomial<double>(std::vector<double>(1, 1.)));
685 else
686 {
687 // create array of Lagrange
688 // polynomials
689 std::vector<Polynomial<double>> v;
690 for (unsigned int i = 0; i <= degree; ++i)
691 v.push_back(LagrangeEquidistant(degree, i));
692 return v;
693 }
694 }
695
696
697
698 //----------------------------------------------------------------------//
699
700
701 std::vector<Polynomial<double>>
702 generate_complete_Lagrange_basis(const std::vector<Point<1>> &points)
703 {
704 std::vector<Polynomial<double>> p;
705 p.reserve(points.size());
706
707 for (unsigned int i = 0; i < points.size(); ++i)
708 p.emplace_back(points, i);
709 return p;
710 }
711
712
713
714 // ------------------ class Legendre --------------- //
715
716
717
718 Legendre::Legendre(const unsigned int k)
719 : Polynomial<double>(0)
720 {
721 this->coefficients.clear();
722 this->in_lagrange_product_form = true;
723 this->lagrange_support_points.resize(k);
724
725 // the roots of a Legendre polynomial are exactly the points in the
726 // Gauss-Legendre quadrature formula
727 if (k > 0)
728 {
729 QGauss<1> gauss(k);
730 for (unsigned int i = 0; i < k; ++i)
731 this->lagrange_support_points[i] = gauss.get_points()[i][0];
732 }
733
734 // compute the abscissa in zero of the product of monomials. The exact
735 // value should be sqrt(2*k+1), so set the weight to that value.
736 double prod = 1.;
737 for (unsigned int i = 0; i < k; ++i)
738 prod *= this->lagrange_support_points[i];
739 this->lagrange_weight = std::sqrt(double(2 * k + 1)) / prod;
740 }
741
742
743
744 std::vector<Polynomial<double>>
745 Legendre::generate_complete_basis(const unsigned int degree)
746 {
747 std::vector<Polynomial<double>> v;
748 v.reserve(degree + 1);
749 for (unsigned int i = 0; i <= degree; ++i)
750 v.push_back(Legendre(i));
751 return v;
752 }
753
754
755
756 // ------------------ class Lobatto -------------------- //
757
758
759 Lobatto::Lobatto(const unsigned int p)
760 : Polynomial<double>(compute_coefficients(p))
761 {}
762
763 std::vector<double>
764 Lobatto::compute_coefficients(const unsigned int p)
765 {
766 switch (p)
767 {
768 case 0:
769 {
770 std::vector<double> coefficients(2);
771
772 coefficients[0] = 1.0;
773 coefficients[1] = -1.0;
774 return coefficients;
775 }
776
777 case 1:
778 {
779 std::vector<double> coefficients(2);
780
781 coefficients[0] = 0.0;
782 coefficients[1] = 1.0;
783 return coefficients;
784 }
785
786 case 2:
787 {
788 std::vector<double> coefficients(3);
789
790 coefficients[0] = 0.0;
791 coefficients[1] = -1.0 * std::sqrt(3.);
792 coefficients[2] = std::sqrt(3.);
793 return coefficients;
794 }
795
796 default:
797 {
798 std::vector<double> coefficients(p + 1);
799 std::vector<double> legendre_coefficients_tmp1(p);
800 std::vector<double> legendre_coefficients_tmp2(p - 1);
801
802 coefficients[0] = -1.0 * std::sqrt(3.);
803 coefficients[1] = 2.0 * std::sqrt(3.);
804 legendre_coefficients_tmp1[0] = 1.0;
805
806 for (unsigned int i = 2; i < p; ++i)
807 {
808 for (unsigned int j = 0; j < i - 1; ++j)
809 legendre_coefficients_tmp2[j] = legendre_coefficients_tmp1[j];
810
811 for (unsigned int j = 0; j < i; ++j)
812 legendre_coefficients_tmp1[j] = coefficients[j];
813
814 coefficients[0] =
815 std::sqrt(2 * i + 1.) *
816 ((1.0 - 2 * i) * legendre_coefficients_tmp1[0] /
817 std::sqrt(2 * i - 1.) +
818 (1.0 - i) * legendre_coefficients_tmp2[0] /
819 std::sqrt(2 * i - 3.)) /
820 i;
821
822 for (unsigned int j = 1; j < i - 1; ++j)
823 coefficients[j] =
824 std::sqrt(2 * i + 1.) *
825 (std::sqrt(2 * i - 1.) *
826 (2.0 * legendre_coefficients_tmp1[j - 1] -
827 legendre_coefficients_tmp1[j]) +
828 (1.0 - i) * legendre_coefficients_tmp2[j] /
829 std::sqrt(2 * i - 3.)) /
830 i;
831
832 coefficients[i - 1] = std::sqrt(4 * i * i - 1.) *
833 (2.0 * legendre_coefficients_tmp1[i - 2] -
834 legendre_coefficients_tmp1[i - 1]) /
835 i;
836 coefficients[i] = 2.0 * std::sqrt(4 * i * i - 1.) *
837 legendre_coefficients_tmp1[i - 1] / i;
838 }
839
840 for (int i = p; i > 0; --i)
841 coefficients[i] = coefficients[i - 1] / i;
842
843 coefficients[0] = 0.0;
844 return coefficients;
845 }
846 }
847 }
848
849 std::vector<Polynomial<double>>
851 {
852 std::vector<Polynomial<double>> basis(p + 1);
853
854 for (unsigned int i = 0; i <= p; ++i)
855 basis[i] = Lobatto(i);
856
857 return basis;
858 }
859
860
861
862 // ------------------ class Hierarchical --------------- //
863
864
865 // Reserve space for polynomials up to degree 19. Should be sufficient
866 // for the start.
867 std::vector<std::unique_ptr<const std::vector<double>>>
869
870
871
872 Hierarchical::Hierarchical(const unsigned int k)
873 : Polynomial<double>(get_coefficients(k))
874 {}
875
876
877
878 void
880 {
881 unsigned int k = k_;
882
883 // first make sure that no other
884 // thread intercepts the operation
885 // of this function
886 // for this, acquire the lock
887 // until we quit this function
888 std::lock_guard<std::mutex> lock(coefficients_lock);
889
890 // The first 2 coefficients
891 // are hard-coded
892 if (k == 0)
893 k = 1;
894 // check: does the information
895 // already exist?
896 if ((recursive_coefficients.size() < k + 1) ||
897 (recursive_coefficients[k].get() == nullptr))
898 // no, then generate the
899 // respective coefficients
900 {
901 // make sure that there is enough
902 // space in the array for the
903 // coefficients, so we have to resize
904 // it to size k+1
905
906 // but it's more complicated than
907 // that: we call this function
908 // recursively, so if we simply
909 // resize it to k+1 here, then
910 // compute the coefficients for
911 // degree k-1 by calling this
912 // function recursively, then it will
913 // reset the size to k -- not enough
914 // for what we want to do below. the
915 // solution therefore is to only
916 // resize the size if we are going to
917 // *increase* it
918 if (recursive_coefficients.size() < k + 1)
919 recursive_coefficients.resize(k + 1);
920
921 if (k <= 1)
922 {
923 // create coefficients
924 // vectors for k=0 and k=1
925 //
926 // allocate the respective
927 // amount of memory and
928 // later assign it to the
929 // coefficients array to
930 // make it const
931 std::vector<double> c0(2);
932 c0[0] = 1.;
933 c0[1] = -1.;
934
935 std::vector<double> c1(2);
936 c1[0] = 0.;
937 c1[1] = 1.;
938
939 // now make these arrays
940 // const
942 std::make_unique<const std::vector<double>>(std::move(c0));
944 std::make_unique<const std::vector<double>>(std::move(c1));
945 }
946 else if (k == 2)
947 {
948 coefficients_lock.unlock();
950 coefficients_lock.lock();
951
952 std::vector<double> c2(3);
953
954 const double a = 1.; // 1./8.;
955
956 c2[0] = 0. * a;
957 c2[1] = -4. * a;
958 c2[2] = 4. * a;
959
961 std::make_unique<const std::vector<double>>(std::move(c2));
962 }
963 else
964 {
965 // for larger numbers,
966 // compute the coefficients
967 // recursively. to do so,
968 // we have to release the
969 // lock temporarily to
970 // allow the called
971 // function to acquire it
972 // itself
973 coefficients_lock.unlock();
975 coefficients_lock.lock();
976
977 std::vector<double> ck(k + 1);
978
979 const double a = 1.; // 1./(2.*k);
980
981 ck[0] = -a * (*recursive_coefficients[k - 1])[0];
982
983 for (unsigned int i = 1; i <= k - 1; ++i)
984 ck[i] = a * (2. * (*recursive_coefficients[k - 1])[i - 1] -
985 (*recursive_coefficients[k - 1])[i]);
986
987 ck[k] = a * 2. * (*recursive_coefficients[k - 1])[k - 1];
988 // for even degrees, we need
989 // to add a multiple of
990 // basis fcn phi_2
991 if ((k % 2) == 0)
992 {
993 double b = 1.; // 8.;
994 // for (unsigned int i=1; i<=k; ++i)
995 // b /= 2.*i;
996
997 ck[1] += b * (*recursive_coefficients[2])[1];
998 ck[2] += b * (*recursive_coefficients[2])[2];
999 }
1000 // finally assign the newly
1001 // created vector to the
1002 // const pointer in the
1003 // coefficients array
1005 std::make_unique<const std::vector<double>>(std::move(ck));
1006 }
1007 }
1008 }
1009
1010
1011
1012 const std::vector<double> &
1014 {
1015 // first make sure the coefficients
1016 // get computed if so necessary
1018
1019 // then get a pointer to the array
1020 // of coefficients. do that in a MT
1021 // safe way
1022 std::lock_guard<std::mutex> lock(coefficients_lock);
1023 return *recursive_coefficients[k];
1024 }
1025
1026
1027
1028 std::vector<Polynomial<double>>
1029 Hierarchical::generate_complete_basis(const unsigned int degree)
1030 {
1031 if (degree == 0)
1032 // create constant
1033 // polynomial. note that we
1034 // can't use the other branch
1035 // of the if-statement, since
1036 // calling the constructor of
1037 // this class with argument
1038 // zero does _not_ create the
1039 // constant polynomial, but
1040 // rather 1-x
1041 return std::vector<Polynomial<double>>(
1042 1, Polynomial<double>(std::vector<double>(1, 1.)));
1043 else
1044 {
1045 std::vector<Polynomial<double>> v;
1046 v.reserve(degree + 1);
1047 for (unsigned int i = 0; i <= degree; ++i)
1048 v.push_back(Hierarchical(i));
1049 return v;
1050 }
1051 }
1052
1053
1054
1055 // ------------------ HermiteInterpolation --------------- //
1056
1058 : Polynomial<double>(0)
1059 {
1060 this->coefficients.clear();
1061 this->in_lagrange_product_form = true;
1062
1063 this->lagrange_support_points.resize(3);
1064 if (p == 0)
1065 {
1066 this->lagrange_support_points[0] = -0.5;
1067 this->lagrange_support_points[1] = 1.;
1068 this->lagrange_support_points[2] = 1.;
1069 this->lagrange_weight = 2.;
1070 }
1071 else if (p == 1)
1072 {
1073 this->lagrange_support_points[0] = 0.;
1074 this->lagrange_support_points[1] = 0.;
1075 this->lagrange_support_points[2] = 1.5;
1076 this->lagrange_weight = -2.;
1077 }
1078 else if (p == 2)
1079 {
1080 this->lagrange_support_points[0] = 0.;
1081 this->lagrange_support_points[1] = 1.;
1082 this->lagrange_support_points[2] = 1.;
1083 }
1084 else if (p == 3)
1085 {
1086 this->lagrange_support_points[0] = 0.;
1087 this->lagrange_support_points[1] = 0.;
1088 this->lagrange_support_points[2] = 1.;
1089 }
1090 else
1091 {
1092 this->lagrange_support_points.resize(4);
1093 this->lagrange_support_points[0] = 0.;
1094 this->lagrange_support_points[1] = 0.;
1095 this->lagrange_support_points[2] = 1.;
1096 this->lagrange_support_points[3] = 1.;
1097 this->lagrange_weight = 16.;
1098
1099 if (p > 4)
1100 {
1101 Legendre legendre(p - 4);
1102 (*this) *= legendre;
1103 }
1104 }
1105 }
1106
1107
1108 std::vector<Polynomial<double>>
1110 {
1111 Assert(n >= 3,
1112 ExcNotImplemented("Hermite interpolation makes no sense for "
1113 "degrees less than three"));
1114 std::vector<Polynomial<double>> basis(n + 1);
1115
1116 for (unsigned int i = 0; i <= n; ++i)
1117 basis[i] = HermiteInterpolation(i);
1118
1119 return basis;
1120 }
1121
1122
1123 // ------------------ HermiteLikeInterpolation --------------- //
1124 namespace
1125 {
1126 // Finds the zero position x_star such that the mass matrix entry (0,1)
1127 // with the Hermite polynomials evaluates to zero. The function has
1128 // originally been derived by a secant method for the integral entry
1129 // l_0(x) * l_1(x) but we only need to do one iteration because the zero
1130 // x_star is linear in the integral value.
1131 double
1132 find_support_point_x_star(const std::vector<double> &jacobi_roots)
1133 {
1134 // Initial guess for the support point position values: The zero turns
1135 // out to be between zero and the first root of the Jacobi polynomial,
1136 // but the algorithm is agnostic about that, so simply choose two points
1137 // that are sufficiently far apart.
1138 double guess_left = 0;
1139 double guess_right = 0.5;
1140 const unsigned int degree = jacobi_roots.size() + 3;
1141
1142 // Compute two integrals of the product of l_0(x) * l_1(x)
1143 // l_0(x) =
1144 // (x-y)*(x-jacobi_roots(0))*...*(x-jacobi_roos(degree-4))*(x-1)*(x-1)
1145 // l_1(x) =
1146 // (x-0)*(x-jacobi_roots(0))*...*(x-jacobi_roots(degree-4))*(x-1)*(x-1)
1147 // where y is either guess_left or guess_right for the two integrals.
1148 // Note that the polynomials are not yet normalized here, which is not
1149 // necessary because we are only looking for the x_star where the matrix
1150 // entry is zero, for which the constants do not matter.
1151 QGauss<1> gauss(degree + 1);
1152 double integral_left = 0, integral_right = 0;
1153 for (unsigned int q = 0; q < gauss.size(); ++q)
1154 {
1155 const double x = gauss.point(q)[0];
1156 double poly_val_common = x;
1157 for (unsigned int j = 0; j < degree - 3; ++j)
1158 poly_val_common *= Utilities::fixed_power<2>(x - jacobi_roots[j]);
1159 poly_val_common *= Utilities::fixed_power<4>(x - 1.);
1160 integral_left +=
1161 gauss.weight(q) * (poly_val_common * (x - guess_left));
1162 integral_right +=
1163 gauss.weight(q) * (poly_val_common * (x - guess_right));
1164 }
1165
1166 // compute guess by secant method. Due to linearity in the root x_star,
1167 // this is the correct position after this single step
1168 return guess_right - (guess_right - guess_left) /
1169 (integral_right - integral_left) * integral_right;
1170 }
1171 } // namespace
1172
1173
1174
1176 const unsigned int index)
1177 : Polynomial<double>(0)
1178 {
1179 AssertIndexRange(index, degree + 1);
1180
1181 this->coefficients.clear();
1182 this->in_lagrange_product_form = true;
1183
1184 this->lagrange_support_points.resize(degree);
1185
1186 if (degree == 0)
1187 this->lagrange_weight = 1.;
1188 else if (degree == 1)
1189 {
1190 if (index == 0)
1191 {
1192 this->lagrange_support_points[0] = 1.;
1193 this->lagrange_weight = -1.;
1194 }
1195 else
1196 {
1197 this->lagrange_support_points[0] = 0.;
1198 this->lagrange_weight = 1.;
1199 }
1200 }
1201 else if (degree == 2)
1202 {
1203 if (index == 0)
1204 {
1205 this->lagrange_support_points[0] = 1.;
1206 this->lagrange_support_points[1] = 1.;
1207 this->lagrange_weight = 1.;
1208 }
1209 else if (index == 1)
1210 {
1211 this->lagrange_support_points[0] = 0;
1212 this->lagrange_support_points[1] = 1;
1213 this->lagrange_weight = -2.;
1214 }
1215 else
1216 {
1217 this->lagrange_support_points[0] = 0.;
1218 this->lagrange_support_points[1] = 0.;
1219 this->lagrange_weight = 1.;
1220 }
1221 }
1222 else if (degree == 3)
1223 {
1224 // 4 Polynomials with degree 3
1225 // entries (1,0) and (3,2) of the mass matrix will be equal to 0
1226 //
1227 // | x 0 x x |
1228 // | 0 x x x |
1229 // M = | x x x 0 |
1230 // | x x 0 x |
1231 //
1232 if (index == 0)
1233 {
1234 this->lagrange_support_points[0] = 2. / 7.;
1235 this->lagrange_support_points[1] = 1.;
1236 this->lagrange_support_points[2] = 1.;
1237 this->lagrange_weight = -3.5;
1238 }
1239 else if (index == 1)
1240 {
1241 this->lagrange_support_points[0] = 0.;
1242 this->lagrange_support_points[1] = 1.;
1243 this->lagrange_support_points[2] = 1.;
1244
1245 // this magic value 5.5 is obtained when evaluating the general
1246 // formula below for the degree=3 case
1247 this->lagrange_weight = 5.5;
1248 }
1249 else if (index == 2)
1250 {
1251 this->lagrange_support_points[0] = 0.;
1252 this->lagrange_support_points[1] = 0.;
1253 this->lagrange_support_points[2] = 1.;
1254 this->lagrange_weight = -5.5;
1255 }
1256 else if (index == 3)
1257 {
1258 this->lagrange_support_points[0] = 0.;
1259 this->lagrange_support_points[1] = 0.;
1260 this->lagrange_support_points[2] = 5. / 7.;
1261 this->lagrange_weight = 3.5;
1262 }
1263 }
1264 else
1265 {
1266 // Higher order Polynomials degree>=4: the entries (1,0) and
1267 // (degree,degree-1) of the mass matrix will be equal to 0
1268 //
1269 // | x 0 x x x x x |
1270 // | 0 x x x . . . x x x |
1271 // | x x x 0 0 x x |
1272 // | x x 0 x 0 x x |
1273 // | . . . |
1274 // M = | . . . |
1275 // | . . . |
1276 // | x x 0 0 x x x |
1277 // | x x x x . . . x x 0 |
1278 // | x x x x x 0 x |
1279 //
1280 // We find the inner points as the zeros of the Jacobi polynomials
1281 // with alpha = beta = 4 which is the polynomial with the kernel
1282 // (1-x)^4 (1+x)^4. Since polynomials (1-x)^2 (1+x)^2 are contained
1283 // in every interior polynomial (bubble function), their product
1284 // leads us to the orthogonality condition of the Jacobi(4,4)
1285 // polynomials.
1286
1287 std::vector<double> jacobi_roots =
1288 jacobi_polynomial_roots<double>(degree - 3, 4, 4);
1289 AssertDimension(jacobi_roots.size(), degree - 3);
1290
1291 // iteration from variable support point N with secant method
1292 // initial values
1293
1294 this->lagrange_support_points.resize(degree);
1295 if (index == 0)
1296 {
1297 const double auxiliary_zero =
1298 find_support_point_x_star(jacobi_roots);
1299 this->lagrange_support_points[0] = auxiliary_zero;
1300 for (unsigned int m = 0; m < degree - 3; ++m)
1301 this->lagrange_support_points[m + 1] = jacobi_roots[m];
1302 this->lagrange_support_points[degree - 2] = 1.;
1303 this->lagrange_support_points[degree - 1] = 1.;
1304
1305 // ensure that the polynomial evaluates to one at x=0
1306 this->lagrange_weight = 1. / this->value(0);
1307 }
1308 else if (index == 1)
1309 {
1310 this->lagrange_support_points[0] = 0.;
1311 for (unsigned int m = 0; m < degree - 3; ++m)
1312 this->lagrange_support_points[m + 1] = jacobi_roots[m];
1313 this->lagrange_support_points[degree - 2] = 1.;
1314 this->lagrange_support_points[degree - 1] = 1.;
1315
1316 // Select the weight to make the derivative of the sum of P_0 and
1317 // P_1 in zero to be 0. The derivative in x=0 is simply given by
1318 // p~(0)/auxiliary_zero+p~'(0) + a*p~(0), where p~(x) is the
1319 // Lagrange polynomial in all points except the first one which is
1320 // the same for P_0 and P_1, and a is the weight we seek here. If
1321 // we solve this for a, we obtain the desired property. Since the
1322 // basis is nodal for all interior points, this property ensures
1323 // that the sum of all polynomials with weight 1 is one.
1324 std::vector<Point<1>> points(degree);
1325 double ratio = 1.;
1326 for (unsigned int i = 0; i < degree; ++i)
1327 {
1328 points[i][0] = this->lagrange_support_points[i];
1329 if (i > 0)
1330 ratio *= -this->lagrange_support_points[i];
1331 }
1332 Polynomial<double> helper(points, 0);
1333 std::vector<double> value_and_grad(2);
1334 helper.value(0., value_and_grad);
1335 Assert(std::abs(value_and_grad[0]) > 1e-10,
1336 ExcInternalError("There should not be a zero at x=0."));
1337
1338 const double auxiliary_zero =
1339 find_support_point_x_star(jacobi_roots);
1340 this->lagrange_weight =
1341 (1. / auxiliary_zero - value_and_grad[1] / value_and_grad[0]) /
1342 ratio;
1343 }
1344 else if (index >= 2 && index < degree - 1)
1345 {
1346 this->lagrange_support_points[0] = 0.;
1347 this->lagrange_support_points[1] = 0.;
1348 for (unsigned int m = 0, c = 2; m < degree - 3; ++m)
1349 if (m + 2 != index)
1350 this->lagrange_support_points[c++] = jacobi_roots[m];
1351 this->lagrange_support_points[degree - 2] = 1.;
1352 this->lagrange_support_points[degree - 1] = 1.;
1353
1354 // ensure that the polynomial evaluates to one at the respective
1355 // nodal point
1356 this->lagrange_weight = 1. / this->value(jacobi_roots[index - 2]);
1357 }
1358 else if (index == degree - 1)
1359 {
1360 this->lagrange_support_points[0] = 0.;
1361 this->lagrange_support_points[1] = 0.;
1362 for (unsigned int m = 0; m < degree - 3; ++m)
1363 this->lagrange_support_points[m + 2] = jacobi_roots[m];
1364 this->lagrange_support_points[degree - 1] = 1.;
1365
1366 std::vector<Point<1>> points(degree);
1367 double ratio = 1.;
1368 for (unsigned int i = 0; i < degree; ++i)
1369 {
1370 points[i][0] = this->lagrange_support_points[i];
1371 if (i < degree - 1)
1372 ratio *= 1. - this->lagrange_support_points[i];
1373 }
1374 Polynomial<double> helper(points, degree - 1);
1375 std::vector<double> value_and_grad(2);
1376 helper.value(1., value_and_grad);
1377 Assert(std::abs(value_and_grad[0]) > 1e-10,
1378 ExcInternalError("There should not be a zero at x=1."));
1379
1380 const double auxiliary_zero =
1381 find_support_point_x_star(jacobi_roots);
1382 this->lagrange_weight =
1383 (-1. / auxiliary_zero - value_and_grad[1] / value_and_grad[0]) /
1384 ratio;
1385 }
1386 else if (index == degree)
1387 {
1388 const double auxiliary_zero =
1389 find_support_point_x_star(jacobi_roots);
1390 this->lagrange_support_points[0] = 0.;
1391 this->lagrange_support_points[1] = 0.;
1392 for (unsigned int m = 0; m < degree - 3; ++m)
1393 this->lagrange_support_points[m + 2] = jacobi_roots[m];
1394 this->lagrange_support_points[degree - 1] = 1. - auxiliary_zero;
1395
1396 // ensure that the polynomial evaluates to one at x=1
1397 this->lagrange_weight = 1. / this->value(1.);
1398 }
1399 }
1400 }
1401
1402
1403
1404 std::vector<Polynomial<double>>
1406 {
1407 std::vector<Polynomial<double>> basis(degree + 1);
1408
1409 for (unsigned int i = 0; i <= degree; ++i)
1410 basis[i] = HermiteLikeInterpolation(degree, i);
1411
1412 return basis;
1413 }
1414
1415} // namespace Polynomials
1416
1417// ------------------ explicit instantiations --------------- //
1418
1419#ifndef DOXYGEN
1420namespace Polynomials
1421{
1422 template class Polynomial<float>;
1423 template class Polynomial<double>;
1424 template class Polynomial<long double>;
1425
1426 template void
1427 Polynomial<float>::shift(const float offset);
1428 template void
1429 Polynomial<float>::shift(const double offset);
1430 template void
1431 Polynomial<double>::shift(const double offset);
1432 template void
1433 Polynomial<long double>::shift(const long double offset);
1434 template void
1435 Polynomial<float>::shift(const long double offset);
1436 template void
1437 Polynomial<double>::shift(const long double offset);
1438
1439 template class Monomial<float>;
1440 template class Monomial<double>;
1441 template class Monomial<long double>;
1442} // namespace Polynomials
1443#endif // DOXYGEN
1444
Definition point.h:112
HermiteInterpolation(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
HermiteLikeInterpolation(const unsigned int degree, const unsigned int index)
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Definition polynomial.h:574
Hierarchical(const unsigned int p)
static void compute_coefficients(const unsigned int p)
static const std::vector< double > & get_coefficients(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
LagrangeEquidistant(const unsigned int n, const unsigned int support_point)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Legendre(const unsigned int p)
std::vector< double > compute_coefficients(const unsigned int p)
Lobatto(const unsigned int p=0)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Monomial(const unsigned int n, const double coefficient=1.)
static std::vector< number > make_vector(unsigned int n, const double coefficient)
number value(const number x) const
Definition polynomial.h:817
bool operator==(const Polynomial< number > &p) const
std::vector< number > coefficients
Definition polynomial.h:302
Polynomial< number > primitive() const
Polynomial< number > & operator+=(const Polynomial< number > &p)
Polynomial< number > derivative() const
void scale(const number factor)
Polynomial< number > & operator-=(const Polynomial< number > &p)
std::vector< number > lagrange_support_points
Definition polynomial.h:314
void shift(const number2 offset)
void print(std::ostream &out) const
static void multiply(std::vector< number > &coefficients, const number factor)
Polynomial< number > & operator*=(const double s)
virtual std::size_t memory_consumption() const
unsigned int degree() const
Definition polynomial.h:800
const std::vector< Point< dim > > & get_points() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > center
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
std::vector< Point< 1 > > generate_equidistant_unit_points(const unsigned int n)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)