Reference documentation for deal.II version 8.5.1
polynomials_bernstein.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2016 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/polynomials_bernstein.h>
17 
18 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
19 #include <boost/math/special_functions/binomial.hpp>
20 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
21 
22 #include <vector>
23 
24 DEAL_II_NAMESPACE_OPEN
25 
26 namespace
27 {
28  template <typename number>
29  std::vector<number>
30  get_bernstein_coefficients (
31  const unsigned int k, const unsigned int n)
32  {
33  Assert(n>0, ExcMessage("Bernstein polynomial needs to be of degree > 0."));
34  AssertIndexRange(k, n+1);
35  std::vector<number> coeff(n + 1, number(0.0));
36  for (unsigned int i = k; i < n + 1; ++i)
37  {
38  coeff[i] = ((i - k) % 2 == 0 ? 1 : -1)
39  * boost::math::binomial_coefficient<number>(n, i)
40  * boost::math::binomial_coefficient<number>(i, k);
41  }
42  return coeff;
43  }
44 }
45 
46 template <typename number>
48  const unsigned int index, const unsigned int degree)
49  :
50  Polynomials::Polynomial<number>(
51  get_bernstein_coefficients<number>(index, degree))
52 {
53 }
54 
55 
56 #include "polynomials_bernstein.inst"
57 
58 DEAL_II_NAMESPACE_CLOSE
#define AssertIndexRange(index, range)
Definition: exceptions.h:1170
PolynomialsBernstein(const unsigned int index, const unsigned int degree)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:313