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