Reference documentation for deal.II version 9.4.1
\(\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// Copyright (C) 2004 - 2021 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>
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
43template <int dim>
45 : FE_PolyTensor<dim>(
47 FiniteElementData<dim>(get_dpo_vector(),
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 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
77template <int dim>
78std::string
80{
81 std::ostringstream namebuf;
82 namebuf << "FE_BR<" << dim << ">(" << 1 << ")";
83
84 return namebuf.str();
85}
86
87
88template <int dim>
89void
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
101template <int dim>
102std::unique_ptr<FiniteElement<dim, dim>>
104{
105 return std::make_unique<FE_BernardiRaugel<dim>>(*this);
106}
107
108
109
110template <int dim>
111void
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 (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
150template <int dim>
151std::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
167template <int dim>
168void
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
200template class FE_BernardiRaugel<2>;
201template class FE_BernardiRaugel<3>;
202
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
Definition: tensor.h:503
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ mapping_none
Definition: mapping.h:77
STL namespace.