Reference documentation for deal.II version GIT 6823c776d5 2023-12-04 20:20: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\}}\)
fe_bernardi_raugel.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2023 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 
20 
22 
23 #include <deal.II/fe/fe.h>
25 #include <deal.II/fe/fe_tools.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/fe/mapping.h>
28 
29 #include <deal.II/grid/tria.h>
31 
32 #include <iostream>
33 #include <memory>
34 #include <sstream>
35 
36 
38 
39 // TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
40 // adjust_line_dof_index_for_line_orientation_table fields, and write tests
41 // similar to bits/face_orientation_and_fe_q_*
42 
43 template <int dim>
45  : FE_PolyTensor<dim>(
48  dim,
49  2,
50  FiniteElementData<dim>::Hdiv),
51  std::vector<bool>(PolynomialsBernardiRaugel<dim>::n_polynomials(p), true),
52  std::vector<ComponentMask>(PolynomialsBernardiRaugel<dim>::n_polynomials(
53  p),
54  ComponentMask(std::vector<bool>(dim, true))))
55 {
56  Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim));
57  Assert(p == 1, ExcMessage("Only BR1 elements are available"));
58 
59  // const unsigned int n_dofs = this->n_dofs_per_cell();
60 
61  this->mapping_kind = {mapping_none};
62  // These must be done first, since
63  // they change the evaluation of
64  // basis functions
65 
66  // Set up the generalized support
67  // points
69 
70  // We need to initialize the dof permutation table and the one for the sign
71  // change.
73 }
74 
75 
76 
77 template <int dim>
78 std::string
80 {
81  std::ostringstream namebuf;
82  namebuf << "FE_BR<" << dim << ">(" << 1 << ")";
83 
84  return namebuf.str();
85 }
86 
87 
88 template <int dim>
89 void
91 {
92  // for 1d and 2d, do nothing
93  if (dim < 3)
94  return;
95 
96  // TODO: Implement this for this class
97  return;
98 }
99 
100 
101 template <int dim>
102 std::unique_ptr<FiniteElement<dim, dim>>
104 {
105  return std::make_unique<FE_BernardiRaugel<dim>>(*this);
106 }
107 
108 
109 
110 template <int dim>
111 void
113  const std::vector<Vector<double>> &support_point_values,
114  std::vector<double> &nodal_values) const
115 {
116  Assert(support_point_values.size() == this->generalized_support_points.size(),
117  ExcDimensionMismatch(support_point_values.size(),
118  this->generalized_support_points.size()));
119  AssertDimension(support_point_values[0].size(), dim);
120  Assert(nodal_values.size() == this->n_dofs_per_cell(),
121  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
122 
123  std::vector<Tensor<1, dim>> normals;
124  for (const unsigned int i : GeometryInfo<dim>::face_indices())
125  {
126  Tensor<1, dim> normal;
127  normal[i / 2] = 1;
128  normals.push_back(normal);
129  }
130 
131  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
132  nodal_values[i] = support_point_values[i][i % dim];
133 
134  // Compute the support points values for the bubble functions
135  for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
136  i < dim * GeometryInfo<dim>::vertices_per_cell +
138  ++i)
139  {
140  nodal_values[i] = 0;
141  for (unsigned int j = 0; j < dim; ++j)
142  nodal_values[i] +=
143  support_point_values[i][j] *
144  normals[i - dim * GeometryInfo<dim>::vertices_per_cell][j];
145  }
146 }
147 
148 
149 
150 template <int dim>
151 std::vector<unsigned int>
153 {
154  // compute the number of unknowns per cell interior/face/edge
155  //
156  // there are <tt>dim</tt> degrees of freedom per vertex and there
157  // is 1 degree of freedom per edge in 2d (face in 3d)
158  std::vector<unsigned int> dpo(dim + 1, 0u);
159  dpo[0] = dim;
160  dpo[dim - 1] = 1u;
161 
162  return dpo;
163 }
164 
165 
166 
167 template <int dim>
168 void
170 {
171  // The support points for our shape functions are the vertices and
172  // the face midpoints, for a total of #vertices + #faces points
173  this->generalized_support_points.resize(this->n_dofs_per_cell());
174 
175  // We need dim copies of each vertex for the first dim*vertices_per_cell
176  // generalized support points
177  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
178  this->generalized_support_points[i] =
180 
181  // The remaining 2*dim points are the edge midpoints
182  for (unsigned int i = 0; i < dim; ++i)
183  {
184  for (unsigned int j = 0; j < 2; ++j)
185  {
186  Point<dim> p;
187  p[0] = 0.5;
188  p[1] = 0.5;
189  if (dim == 3)
190  p[2] = 0.5;
191  p[i] = j;
192 
193  const unsigned int k =
194  dim * GeometryInfo<dim>::vertices_per_cell + i * 2 + j;
195  this->generalized_support_points[k] = p;
196  }
197  }
198 }
199 
200 template class FE_BernardiRaugel<2>;
201 template class FE_BernardiRaugel<3>;
202 
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()
void initialize_quad_dof_index_permutation_and_sign_change()
FE_BernardiRaugel(const unsigned int p=1)
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual std::string get_name() const override
std::vector< MappingKind > mapping_kind
Definition: point.h:112
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
static ::ExceptionBase & ExcMessage(std::string arg1)
@ mapping_none
Definition: mapping.h:83
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)