Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
fe_bernardi_raugel.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
19
21
22#include <deal.II/fe/fe.h>
24#include <deal.II/fe/fe_tools.h>
26#include <deal.II/fe/mapping.h>
27
28#include <deal.II/grid/tria.h>
30
31#include <iostream>
32#include <memory>
33#include <sstream>
34
35
37
38// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
39// adjust_line_dof_index_for_line_orientation_table fields, and write tests
40// similar to bits/face_orientation_and_fe_q_*
41
42template <int dim>
44 : FE_PolyTensor<dim>(
46 FiniteElementData<dim>(get_dpo_vector(),
47 dim,
48 2,
49 FiniteElementData<dim>::Hdiv),
50 std::vector<bool>(PolynomialsBernardiRaugel<dim>::n_polynomials(p), true),
51 std::vector<ComponentMask>(PolynomialsBernardiRaugel<dim>::n_polynomials(
52 p),
53 ComponentMask(std::vector<bool>(dim, true))))
54{
55 Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim));
56 Assert(p == 1, ExcMessage("Only BR1 elements are available"));
57
58 // const unsigned int n_dofs = this->n_dofs_per_cell();
59
60 this->mapping_kind = {mapping_none};
61 // These must be done first, since
62 // they change the evaluation of
63 // basis functions
64
65 // Set up the generalized support
66 // points
68
69 // We need to initialize the dof permutation table and the one for the sign
70 // change.
72}
73
74
75
76template <int dim>
77std::string
79{
80 std::ostringstream namebuf;
81 namebuf << "FE_BR<" << dim << ">(" << 1 << ")";
82
83 return namebuf.str();
84}
85
86
87template <int dim>
88void
90{
91 // for 1d and 2d, do nothing
92 if (dim < 3)
93 return;
94
95 // TODO: Implement this for this class
96 return;
97}
98
99
100template <int dim>
101std::unique_ptr<FiniteElement<dim, dim>>
103{
104 return std::make_unique<FE_BernardiRaugel<dim>>(*this);
105}
106
107
108
109template <int dim>
110void
112 const std::vector<Vector<double>> &support_point_values,
113 std::vector<double> &nodal_values) const
114{
115 Assert(support_point_values.size() == this->generalized_support_points.size(),
116 ExcDimensionMismatch(support_point_values.size(),
117 this->generalized_support_points.size()));
118 AssertDimension(support_point_values[0].size(), dim);
119 Assert(nodal_values.size() == this->n_dofs_per_cell(),
120 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
121
122 std::vector<Tensor<1, dim>> normals;
123 for (const unsigned int i : GeometryInfo<dim>::face_indices())
124 {
125 Tensor<1, dim> normal;
126 normal[i / 2] = 1;
127 normals.push_back(normal);
128 }
129
130 for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
131 nodal_values[i] = support_point_values[i][i % dim];
132
133 // Compute the support points values for the bubble functions
134 for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
135 i < dim * GeometryInfo<dim>::vertices_per_cell +
137 ++i)
138 {
139 nodal_values[i] = 0;
140 for (unsigned int j = 0; j < dim; ++j)
141 nodal_values[i] +=
142 support_point_values[i][j] *
143 normals[i - dim * GeometryInfo<dim>::vertices_per_cell][j];
144 }
145}
146
147
148
149template <int dim>
150std::vector<unsigned int>
152{
153 // compute the number of unknowns per cell interior/face/edge
154 //
155 // there are <tt>dim</tt> degrees of freedom per vertex and there
156 // is 1 degree of freedom per edge in 2d (face in 3d)
157 std::vector<unsigned int> dpo(dim + 1, 0u);
158 dpo[0] = dim;
159 dpo[dim - 1] = 1u;
160
161 return dpo;
162}
163
164
165
166template <int dim>
167void
169{
170 // The support points for our shape functions are the vertices and
171 // the face midpoints, for a total of #vertices + #faces points
172 this->generalized_support_points.resize(this->n_dofs_per_cell());
173
174 // We need dim copies of each vertex for the first dim*vertices_per_cell
175 // generalized support points
176 for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
177 this->generalized_support_points[i] =
179
180 // The remaining 2*dim points are the edge midpoints
181 for (unsigned int i = 0; i < dim; ++i)
182 {
183 for (unsigned int j = 0; j < 2; ++j)
184 {
185 Point<dim> p;
186 p[0] = 0.5;
187 p[1] = 0.5;
188 if (dim == 3)
189 p[2] = 0.5;
190 p[i] = j;
191
192 const unsigned int k =
193 dim * GeometryInfo<dim>::vertices_per_cell + i * 2 + j;
194 this->generalized_support_points[k] = p;
195 }
196 }
197}
198
199template class FE_BernardiRaugel<2>;
200template class FE_BernardiRaugel<3>;
201
static std::vector< unsigned int > get_dpo_vector()
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
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:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ mapping_none
Definition mapping.h:84
STL namespace.
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()