Reference documentation for deal.II version 9.4.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\}}\)
fe_p1nc.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2015 - 2022 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/fe/fe_p1nc.h>
18
19#include <memory>
20
22
23
25 : FiniteElement<2, 2>(
27 std::vector<bool>(1, false),
28 std::vector<ComponentMask>(4, ComponentMask(1, true)))
29{
30 // face support points: 2 end vertices
31 unit_face_support_points[0].resize(2);
32 unit_face_support_points[0][0][0] = 0.0;
33 unit_face_support_points[0][1][0] = 1.0;
34
35 // initialize constraints matrix
37}
38
39
40
41std::string
43{
44 return "FE_P1NC";
45}
46
47
48
51{
53
54 if (flags & update_values)
56 if (flags & update_gradients)
57 out |= update_gradients;
58 if (flags & update_normal_vectors)
60 if (flags & update_hessians)
61 out |= update_hessians;
62
63 return out;
64}
65
66
67
68std::unique_ptr<FiniteElement<2, 2>>
70{
71 return std::make_unique<FE_P1NC>(*this);
72}
73
74
75
76std::vector<unsigned int>
78{
79 std::vector<unsigned int> dpo(3);
80 dpo[0] = 1; // dofs per object: vertex
81 dpo[1] = 0; // line
82 dpo[2] = 0; // quad
83 return dpo;
84}
85
86
87
91{
92 // edge midpoints
93 const Point<2> mpt[4] = {(cell->vertex(0) + cell->vertex(2)) / 2,
94 (cell->vertex(1) + cell->vertex(3)) / 2,
95 (cell->vertex(0) + cell->vertex(1)) / 2,
96 (cell->vertex(2) + cell->vertex(3)) / 2};
97
98 // center point
99 const Point<2> cpt =
100 (cell->vertex(0) + cell->vertex(1) + cell->vertex(2) + cell->vertex(3)) / 4;
101
102 const double det = (mpt[0](0) - mpt[1](0)) * (mpt[2](1) - mpt[3](1)) -
103 (mpt[2](0) - mpt[3](0)) * (mpt[0](1) - mpt[1](1));
104
106 coeffs[0][0] =
107 ((mpt[2](1) - mpt[3](1)) * (0.5) - (mpt[0](1) - mpt[1](1)) * (0.5)) / det;
108 coeffs[1][0] =
109 ((mpt[2](1) - mpt[3](1)) * (-0.5) - (mpt[0](1) - mpt[1](1)) * (0.5)) / det;
110 coeffs[2][0] =
111 ((mpt[2](1) - mpt[3](1)) * (0.5) - (mpt[0](1) - mpt[1](1)) * (-0.5)) / det;
112 coeffs[3][0] =
113 ((mpt[2](1) - mpt[3](1)) * (-0.5) - (mpt[0](1) - mpt[1](1)) * (-0.5)) / det;
114
115 coeffs[0][1] =
116 (-(mpt[2](0) - mpt[3](0)) * (0.5) + (mpt[0](0) - mpt[1](0)) * (0.5)) / det;
117 coeffs[1][1] =
118 (-(mpt[2](0) - mpt[3](0)) * (-0.5) + (mpt[0](0) - mpt[1](0)) * (0.5)) / det;
119 coeffs[2][1] =
120 (-(mpt[2](0) - mpt[3](0)) * (0.5) + (mpt[0](0) - mpt[1](0)) * (-0.5)) / det;
121 coeffs[3][1] =
122 (-(mpt[2](0) - mpt[3](0)) * (-0.5) + (mpt[0](0) - mpt[1](0)) * (-0.5)) /
123 det;
124
125 coeffs[0][2] = 0.25 - cpt(0) * coeffs[0][0] - cpt(1) * coeffs[0][1];
126 coeffs[1][2] = 0.25 - cpt(0) * coeffs[1][0] - cpt(1) * coeffs[1][1];
127 coeffs[2][2] = 0.25 - cpt(0) * coeffs[2][0] - cpt(1) * coeffs[2][1];
128 coeffs[3][2] = 0.25 - cpt(0) * coeffs[3][0] - cpt(1) * coeffs[3][1];
129
130 return coeffs;
131}
132
133
134
135std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
137 const UpdateFlags update_flags,
138 const Mapping<2, 2> &,
139 const Quadrature<2> &quadrature,
141 &output_data) const
142{
143 auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
144
145 data_ptr->update_each = requires_update_flags(update_flags);
146
147 const unsigned int n_q_points = quadrature.size();
148 output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
149
150 // this is a linear element, so its second derivatives are zero
151 if (data_ptr->update_each & update_hessians)
152 output_data.shape_hessians.fill(Tensor<2, 2>());
153
154 return data_ptr;
155}
156
157
158
159std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
161 const UpdateFlags update_flags,
162 const Mapping<2, 2> &,
163 const hp::QCollection<1> &quadrature,
165 &output_data) const
166{
167 AssertDimension(quadrature.size(), 1);
168
169 auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
170
171 data_ptr->update_each = requires_update_flags(update_flags);
172
173 const unsigned int n_q_points = quadrature[0].size();
174 output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
175
176 // this is a linear element, so its second derivatives are zero
177 if (data_ptr->update_each & update_hessians)
178 output_data.shape_hessians.fill(Tensor<2, 2>());
179
180 return data_ptr;
181}
182
183
184
185std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
187 const UpdateFlags update_flags,
188 const Mapping<2, 2> &,
189 const Quadrature<1> &quadrature,
191 &output_data) const
192{
193 auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
194
195 data_ptr->update_each = requires_update_flags(update_flags);
196
197 const unsigned int n_q_points = quadrature.size();
198 output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
199
200 // this is a linear element, so its second derivatives are zero
201 if (data_ptr->update_each & update_hessians)
202 output_data.shape_hessians.fill(Tensor<2, 2>());
203
204 return data_ptr;
205}
206
207
208
209void
213 const Quadrature<2> &,
214 const Mapping<2, 2> &,
217 & mapping_data,
218 const FiniteElement<2, 2>::InternalDataBase &fe_internal,
220 const
221{
222 const UpdateFlags flags(fe_internal.update_each);
223
224 const unsigned int n_q_points = mapping_data.quadrature_points.size();
225
226 // linear shape functions
228
229 // compute on the cell
230 if (flags & update_values)
231 for (unsigned int i = 0; i < n_q_points; ++i)
232 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
233 output_data.shape_values[k][i] =
234 (coeffs[k][0] * mapping_data.quadrature_points[i](0) +
235 coeffs[k][1] * mapping_data.quadrature_points[i](1) + coeffs[k][2]);
236
237 if (flags & update_gradients)
238 for (unsigned int i = 0; i < n_q_points; ++i)
239 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
240 output_data.shape_gradients[k][i] =
241 Point<2>(coeffs[k][0], coeffs[k][1]);
242}
243
244
245
246void
249 const unsigned int face_no,
250 const hp::QCollection<1> & quadrature,
251 const Mapping<2, 2> & mapping,
253 const ::internal::FEValuesImplementation::MappingRelatedData<2, 2> &,
254 const InternalDataBase &fe_internal,
256 &output_data) const
257{
258 AssertDimension(quadrature.size(), 1);
259
260 const UpdateFlags flags(fe_internal.update_each);
261
262 // linear shape functions
264
265 // compute on the face
266 const Quadrature<2> quadrature_on_face =
268 quadrature[0],
269 face_no);
270
271 if (flags & update_values)
272 for (unsigned int i = 0; i < quadrature_on_face.size(); ++i)
273 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
274 {
275 const Point<2> quadrature_point =
276 mapping.transform_unit_to_real_cell(cell,
277 quadrature_on_face.point(i));
278
279 output_data.shape_values[k][i] =
280 (coeffs[k][0] * quadrature_point(0) +
281 coeffs[k][1] * quadrature_point(1) + coeffs[k][2]);
282 }
283
284 if (flags & update_gradients)
285 for (unsigned int i = 0; i < quadrature_on_face.size(); ++i)
286 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
287 output_data.shape_gradients[k][i] =
288 Point<2>(coeffs[k][0], coeffs[k][1]);
289}
290
291
292
293void
296 const unsigned int face_no,
297 const unsigned int sub_no,
298 const Quadrature<1> & quadrature,
299 const Mapping<2, 2> & mapping,
301 const ::internal::FEValuesImplementation::MappingRelatedData<2, 2> &,
302 const InternalDataBase &fe_internal,
304 &output_data) const
305{
306 const UpdateFlags flags(fe_internal.update_each);
307
308 // linear shape functions
310
311 // compute on the subface
312 const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(
313 this->reference_cell(), quadrature, face_no, sub_no);
314
315 if (flags & update_values)
316 for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i)
317 {
318 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
319 {
320 const Point<2> quadrature_point =
322 cell, quadrature_on_subface.point(i));
323
324 output_data.shape_values[k][i] =
325 (coeffs[k][0] * quadrature_point(0) +
326 coeffs[k][1] * quadrature_point(1) + coeffs[k][2]);
327 }
328 }
329
330 if (flags & update_gradients)
331 for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i)
332 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
333 output_data.shape_gradients[k][i] =
334 Point<2>(coeffs[k][0], coeffs[k][1]);
335}
336
337
338
339void
341{
342 // coefficient relation between children and mother
343 interface_constraints.TableBase<2, double>::reinit(
345
346 interface_constraints(0, 0) = 0.5;
347 interface_constraints(0, 1) = 0.5;
348}
349
350
virtual void fill_fe_face_values(const Triangulation< 2, 2 >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< 1 > &quadrature, const Mapping< 2, 2 > &mapping, const Mapping< 2, 2 >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< 2, 2 > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< 2, 2 > &output_data) const override
Definition: fe_p1nc.cc:247
FE_P1NC()
Definition: fe_p1nc.cc:24
virtual std::unique_ptr< FiniteElement< 2, 2 >::InternalDataBase > get_face_data(const UpdateFlags update_flags, const Mapping< 2, 2 > &, const hp::QCollection< 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< 2, 2 > &output_data) const override
Definition: fe_p1nc.cc:160
static ndarray< double, 4, 3 > get_linear_shape_coefficients(const Triangulation< 2, 2 >::cell_iterator &cell)
Definition: fe_p1nc.cc:89
virtual std::unique_ptr< FiniteElement< 2, 2 > > clone() const override
Definition: fe_p1nc.cc:69
virtual void fill_fe_subface_values(const Triangulation< 2, 2 >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< 1 > &quadrature, const Mapping< 2, 2 > &mapping, const Mapping< 2, 2 >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< 2, 2 > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< 2, 2 > &output_data) const override
Definition: fe_p1nc.cc:294
virtual void fill_fe_values(const Triangulation< 2, 2 >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< 2 > &quadrature, const Mapping< 2, 2 > &mapping, const Mapping< 2, 2 >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< 2, 2 > &mapping_data, const FiniteElement< 2, 2 >::InternalDataBase &fe_internal, internal::FEValuesImplementation::FiniteElementRelatedData< 2, 2 > &output_data) const override
Definition: fe_p1nc.cc:210
void initialize_constraints()
Definition: fe_p1nc.cc:340
virtual std::unique_ptr< FiniteElement< 2, 2 >::InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< 2, 2 > &, const Quadrature< 2 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< 2, 2 > &output_data) const override
Definition: fe_p1nc.cc:136
virtual std::unique_ptr< FiniteElement< 2, 2 >::InternalDataBase > get_subface_data(const UpdateFlags update_flags, const Mapping< 2, 2 > &, const Quadrature< 1 > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< 2, 2 > &output_data) const override
Definition: fe_p1nc.cc:186
virtual UpdateFlags requires_update_flags(const UpdateFlags flags) const override
Definition: fe_p1nc.cc:50
static std::vector< unsigned int > get_dpo_vector()
Definition: fe_p1nc.cc:77
virtual std::string get_name() const override
Definition: fe_p1nc.cc:42
unsigned int n_dofs_per_cell() const
ReferenceCell reference_cell() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition: fe.h:2458
TableIndices< 2 > interface_constraints_size() const
FullMatrix< double > interface_constraints
Definition: fe.h:2439
Abstract base class for mapping classes.
Definition: mapping.h:311
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
Definition: tensor.h:503
unsigned int size() const
Definition: collection.h:264
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:2788
std::vector< Point< spacedim > > quadrature_points
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ update_default
No update.
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
STL namespace.
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108