Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
polynomials_rannacher_turek.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/geometry_info.h>
18 #include <deal.II/base/polynomials_rannacher_turek.h>
19 
20 DEAL_II_NAMESPACE_OPEN
21 
22 
23 template <int dim>
25 {
26  Assert(dim == 2, ExcNotImplemented());
27 }
28 
29 
30 
31 template <int dim>
32 double
34  const Point<dim> & p) const
35 {
36  Assert(dim == 2, ExcNotImplemented());
37  if (i == 0)
38  {
39  return (0.75 - 2.5 * p(0) + 1.5 * p(1) +
40  1.5 * (p(0) * p(0) - p(1) * p(1)));
41  }
42  else if (i == 1)
43  {
44  return (-0.25 - 0.5 * p(0) + 1.5 * p(1) +
45  1.5 * (p(0) * p(0) - p(1) * p(1)));
46  }
47  else if (i == 2)
48  {
49  return (0.75 + 1.5 * p(0) - 2.5 * p(1) -
50  1.5 * (p(0) * p(0) - p(1) * p(1)));
51  }
52  else if (i == 3)
53  {
54  return (-0.25 + 1.5 * p(0) - 0.5 * p(1) -
55  1.5 * (p(0) * p(0) - p(1) * p(1)));
56  }
57 
58  Assert(false, ExcNotImplemented());
59  return 0;
60 }
61 
62 
63 
64 template <int dim>
67  const Point<dim> & p) const
68 {
69  Assert(dim == 2, ExcNotImplemented());
70  Tensor<1, dim> grad;
71  if (i == 0)
72  {
73  grad[0] = -2.5 + 3 * p(0);
74  grad[1] = 1.5 - 3 * p(1);
75  }
76  else if (i == 1)
77  {
78  grad[0] = -0.5 + 3.0 * p(0);
79  grad[1] = 1.5 - 3.0 * p(1);
80  }
81  else if (i == 2)
82  {
83  grad[0] = 1.5 - 3.0 * p(0);
84  grad[1] = -2.5 + 3.0 * p(1);
85  }
86  else if (i == 3)
87  {
88  grad[0] = 1.5 - 3.0 * p(0);
89  grad[1] = -0.5 + 3.0 * p(1);
90  }
91  else
92  {
93  Assert(false, ExcNotImplemented());
94  }
95 
96  return grad;
97 }
98 
99 
100 
101 template <int dim>
104  const unsigned int i,
105  const Point<dim> & /*p*/) const
106 {
107  Assert(dim == 2, ExcNotImplemented());
108  Tensor<2, dim> grad_grad;
109  if (i == 0)
110  {
111  grad_grad[0][0] = 3;
112  grad_grad[0][1] = 0;
113  grad_grad[1][0] = 0;
114  grad_grad[1][1] = -3;
115  }
116  else if (i == 1)
117  {
118  grad_grad[0][0] = 3;
119  grad_grad[0][1] = 0;
120  grad_grad[1][0] = 0;
121  grad_grad[1][1] = -3;
122  }
123  else if (i == 2)
124  {
125  grad_grad[0][0] = -3;
126  grad_grad[0][1] = 0;
127  grad_grad[1][0] = 0;
128  grad_grad[1][1] = 3;
129  }
130  else if (i == 3)
131  {
132  grad_grad[0][0] = -3;
133  grad_grad[0][1] = 0;
134  grad_grad[1][0] = 0;
135  grad_grad[1][1] = 3;
136  }
137  return grad_grad;
138 }
139 
140 
141 
142 template <int dim>
143 void
145  const Point<dim> & unit_point,
146  std::vector<double> & values,
147  std::vector<Tensor<1, dim>> &grads,
148  std::vector<Tensor<2, dim>> &grad_grads,
149  std::vector<Tensor<3, dim>> &third_derivatives,
150  std::vector<Tensor<4, dim>> &fourth_derivatives) const
151 {
152  const unsigned int n_pols = ::GeometryInfo<dim>::faces_per_cell;
153  Assert(values.size() == n_pols || values.size() == 0,
154  ExcDimensionMismatch(values.size(), n_pols));
155  Assert(grads.size() == n_pols || grads.size() == 0,
156  ExcDimensionMismatch(grads.size(), n_pols));
157  Assert(grad_grads.size() == n_pols || grad_grads.size() == 0,
158  ExcDimensionMismatch(grad_grads.size(), n_pols));
159  Assert(third_derivatives.size() == n_pols || third_derivatives.size() == 0,
160  ExcDimensionMismatch(third_derivatives.size(), n_pols));
161  Assert(fourth_derivatives.size() == n_pols || fourth_derivatives.size() == 0,
162  ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
163 
164  for (unsigned int i = 0; i < n_pols; ++i)
165  {
166  if (values.size() != 0)
167  {
168  values[i] = compute_value(i, unit_point);
169  }
170  if (grads.size() != 0)
171  {
172  grads[i] = compute_grad(i, unit_point);
173  }
174  if (grad_grads.size() != 0)
175  {
176  grad_grads[i] = compute_grad_grad(i, unit_point);
177  }
178  if (third_derivatives.size() != 0)
179  {
180  third_derivatives[i] = compute_derivative<3>(i, unit_point);
181  }
182  if (fourth_derivatives.size() != 0)
183  {
184  fourth_derivatives[i] = compute_derivative<4>(i, unit_point);
185  }
186  }
187 }
188 
189 
190 // explicit instantiations
191 #include "polynomials_rannacher_turek.inst"
192 
193 DEAL_II_NAMESPACE_CLOSE
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcNotImplemented()
void compute(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim >> &grads, std::vector< Tensor< 2, dim >> &grad_grads, std::vector< Tensor< 3, dim >> &third_derivatives, std::vector< Tensor< 4, dim >> &fourth_derivatives) const