Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_bernardi_raugel.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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/polynomials_p.h>
18 #include <deal.II/base/qprojector.h>
19 #include <deal.II/base/quadrature_lib.h>
20 #include <deal.II/base/std_cxx14/memory.h>
21 
22 #include <deal.II/dofs/dof_accessor.h>
23 
24 #include <deal.II/fe/fe.h>
25 #include <deal.II/fe/fe_bernardi_raugel.h>
26 #include <deal.II/fe/fe_tools.h>
27 #include <deal.II/fe/fe_values.h>
28 #include <deal.II/fe/mapping.h>
29 
30 #include <deal.II/grid/tria.h>
31 #include <deal.II/grid/tria_iterator.h>
32 
33 #include <iostream>
34 #include <sstream>
35 
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 template <int dim>
42  p,
43  FiniteElementData<dim>(get_dpo_vector(),
44  dim,
45  2,
46  FiniteElementData<dim>::Hdiv),
47  std::vector<bool>(PolynomialsBernardiRaugel<dim>::compute_n_pols(p),
48  true),
49  std::vector<ComponentMask>(PolynomialsBernardiRaugel<dim>::compute_n_pols(
50  p),
51  std::vector<bool>(dim, true)))
52 {
53  Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim));
54  Assert(p == 1, ExcMessage("Only BR1 elements are available"));
55 
56  // const unsigned int n_dofs = this->dofs_per_cell;
57 
58  this->mapping_type = mapping_none;
59  // These must be done first, since
60  // they change the evaluation of
61  // basis functions
62 
63  // Set up the generalized support
64  // points
66 }
67 
68 
69 
70 template <int dim>
71 std::string
73 {
74  std::ostringstream namebuf;
75  namebuf << "FE_BR<" << dim << ">(" << 1 << ")";
76 
77  return namebuf.str();
78 }
79 
80 
81 
82 template <int dim>
83 std::unique_ptr<FiniteElement<dim, dim>>
85 {
86  return std_cxx14::make_unique<FE_BernardiRaugel<dim>>(*this);
87 }
88 
89 
90 
91 template <int dim>
92 void
94  const std::vector<Vector<double>> &support_point_values,
95  std::vector<double> & nodal_values) const
96 {
97  Assert(support_point_values.size() == this->generalized_support_points.size(),
98  ExcDimensionMismatch(support_point_values.size(),
99  this->generalized_support_points.size()));
100  AssertDimension(support_point_values[0].size(), dim);
101  Assert(nodal_values.size() == this->dofs_per_cell,
102  ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
103 
104  std::vector<Tensor<1, dim>> normals;
105  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
106  {
107  Tensor<1, dim> normal;
108  normal[i / 2] = 1;
109  normals.push_back(normal);
110  }
111 
112  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
113  nodal_values[i] = support_point_values[i][i % dim];
114 
115  // Compute the support points values for the bubble functions
116  for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
117  i < dim * GeometryInfo<dim>::vertices_per_cell +
119  ++i)
120  {
121  nodal_values[i] = 0;
122  for (unsigned int j = 0; j < dim; ++j)
123  nodal_values[i] +=
124  support_point_values[i][j] *
125  normals[i - dim * GeometryInfo<dim>::vertices_per_cell][j];
126  }
127 }
128 
129 
130 
131 template <int dim>
132 std::vector<unsigned int>
134 {
135  // compute the number of unknowns per cell interior/face/edge
136  //
137  // there are <tt>dim</tt> degrees of freedom per vertex and there
138  // is 1 degree of freedom per edge in 2D (face in 3D)
139  std::vector<unsigned int> dpo(dim + 1, 0u);
140  dpo[0] = dim;
141  dpo[dim - 1] = 1u;
142 
143  return dpo;
144 }
145 
146 
147 
148 template <int dim>
149 void
151 {
152  // The support points for our shape functions are the vertices and
153  // the face midpoints, for a total of #vertices + #faces points
154  this->generalized_support_points.resize(this->dofs_per_cell);
155 
156  // We need dim copies of each vertex for the first dim*vertices_per_cell
157  // generalized support points
158  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
159  this->generalized_support_points[i] =
161 
162  // The remaining 2*dim points are the edge midpoints
163  for (unsigned int i = 0; i < dim; ++i)
164  {
165  for (unsigned int j = 0; j < 2; ++j)
166  {
167  Point<dim> p;
168  p[0] = 0.5;
169  p[1] = 0.5;
170  if (dim == 3)
171  p[2] = 0.5;
172  p[i] = j;
173 
174  const unsigned int k =
175  dim * GeometryInfo<dim>::vertices_per_cell + i * 2 + j;
176  this->generalized_support_points[k] = p;
177  }
178  }
179 }
180 
181 template class FE_BernardiRaugel<1>;
182 template class FE_BernardiRaugel<2>;
183 template class FE_BernardiRaugel<3>;
184 
185 DEAL_II_NAMESPACE_CLOSE
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
FE_BernardiRaugel(const unsigned int p=1)
virtual std::string get_name() const override
STL namespace.
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double >> &support_point_values, std::vector< double > &nodal_values) const override
static std::vector< unsigned int > get_dpo_vector()
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override