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