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