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