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