Reference documentation for deal.II version GIT relicensing-870-g40dcf97519 2024-06-20 23:40: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\}}\)
Loading...
Searching...
No Matches
fe_p1nc.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2024 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
16#include <deal.II/fe/fe_p1nc.h>
17
18#include <memory>
19
21
22
24 : FiniteElement<2, 2>(
25 FiniteElementData<2>(get_dpo_vector(), 1, 1, FiniteElementData<2>::L2),
26 std::vector<bool>(1, false),
27 std::vector<ComponentMask>(4, ComponentMask(1, true)))
28{
29 // face support points: 2 end vertices
30 unit_face_support_points[0].resize(2);
31 unit_face_support_points[0][0][0] = 0.0;
32 unit_face_support_points[0][1][0] = 1.0;
33
34 // initialize constraints matrix
36}
37
38
39
40std::string
42{
43 return "FE_P1NC";
44}
45
46
47
50{
52
53 if (flags & update_values)
55 if (flags & update_gradients)
56 out |= update_gradients;
57 if (flags & update_normal_vectors)
59 if (flags & update_hessians)
60 out |= update_hessians;
61
62 return out;
63}
64
65
66
67std::unique_ptr<FiniteElement<2, 2>>
69{
70 return std::make_unique<FE_P1NC>(*this);
71}
72
73
74
75std::vector<unsigned int>
77{
78 std::vector<unsigned int> dpo(3);
79 dpo[0] = 1; // dofs per object: vertex
80 dpo[1] = 0; // line
81 dpo[2] = 0; // quad
82 return dpo;
83}
84
85
86
90{
91 // edge midpoints
92 const Point<2> mpt[4] = {(cell->vertex(0) + cell->vertex(2)) / 2,
93 (cell->vertex(1) + cell->vertex(3)) / 2,
94 (cell->vertex(0) + cell->vertex(1)) / 2,
95 (cell->vertex(2) + cell->vertex(3)) / 2};
96
97 // center point
98 const Point<2> cpt =
99 (cell->vertex(0) + cell->vertex(1) + cell->vertex(2) + cell->vertex(3)) / 4;
100
101 const double det = (mpt[0][0] - mpt[1][0]) * (mpt[2][1] - mpt[3][1]) -
102 (mpt[2][0] - mpt[3][0]) * (mpt[0][1] - mpt[1][1]);
103
105 coeffs[0][0] =
106 ((mpt[2][1] - mpt[3][1]) * (0.5) - (mpt[0][1] - mpt[1][1]) * (0.5)) / det;
107 coeffs[1][0] =
108 ((mpt[2][1] - mpt[3][1]) * (-0.5) - (mpt[0][1] - mpt[1][1]) * (0.5)) / det;
109 coeffs[2][0] =
110 ((mpt[2][1] - mpt[3][1]) * (0.5) - (mpt[0][1] - mpt[1][1]) * (-0.5)) / det;
111 coeffs[3][0] =
112 ((mpt[2][1] - mpt[3][1]) * (-0.5) - (mpt[0][1] - mpt[1][1]) * (-0.5)) / det;
113
114 coeffs[0][1] =
115 (-(mpt[2][0] - mpt[3][0]) * (0.5) + (mpt[0][0] - mpt[1][0]) * (0.5)) / det;
116 coeffs[1][1] =
117 (-(mpt[2][0] - mpt[3][0]) * (-0.5) + (mpt[0][0] - mpt[1][0]) * (0.5)) / det;
118 coeffs[2][1] =
119 (-(mpt[2][0] - mpt[3][0]) * (0.5) + (mpt[0][0] - mpt[1][0]) * (-0.5)) / det;
120 coeffs[3][1] =
121 (-(mpt[2][0] - mpt[3][0]) * (-0.5) + (mpt[0][0] - mpt[1][0]) * (-0.5)) /
122 det;
123
124 coeffs[0][2] = 0.25 - cpt[0] * coeffs[0][0] - cpt[1] * coeffs[0][1];
125 coeffs[1][2] = 0.25 - cpt[0] * coeffs[1][0] - cpt[1] * coeffs[1][1];
126 coeffs[2][2] = 0.25 - cpt[0] * coeffs[2][0] - cpt[1] * coeffs[2][1];
127 coeffs[3][2] = 0.25 - cpt[0] * coeffs[3][0] - cpt[1] * coeffs[3][1];
128
129 return coeffs;
130}
131
132
133
134std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
136 const UpdateFlags update_flags,
137 const Mapping<2, 2> &,
138 const Quadrature<2> &quadrature,
140 &output_data) const
141{
142 auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
143
144 data_ptr->update_each = requires_update_flags(update_flags);
145
146 const unsigned int n_q_points = quadrature.size();
147 output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
148
149 // this is a linear element, so its second derivatives are zero
150 if (data_ptr->update_each & update_hessians)
151 output_data.shape_hessians.fill(Tensor<2, 2>());
152
153 return data_ptr;
154}
155
156
157
158std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
160 const UpdateFlags update_flags,
161 const Mapping<2, 2> &,
162 const hp::QCollection<1> &quadrature,
164 &output_data) const
165{
166 AssertDimension(quadrature.size(), 1);
167
168 auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
169
170 data_ptr->update_each = requires_update_flags(update_flags);
171
172 const unsigned int n_q_points = quadrature[0].size();
173 output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
174
175 // this is a linear element, so its second derivatives are zero
176 if (data_ptr->update_each & update_hessians)
177 output_data.shape_hessians.fill(Tensor<2, 2>());
178
179 return data_ptr;
180}
181
182
183
184std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
186 const UpdateFlags update_flags,
187 const Mapping<2, 2> &,
188 const Quadrature<1> &quadrature,
190 &output_data) const
191{
192 auto data_ptr = std::make_unique<FiniteElement<2, 2>::InternalDataBase>();
193
194 data_ptr->update_each = requires_update_flags(update_flags);
195
196 const unsigned int n_q_points = quadrature.size();
197 output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
198
199 // this is a linear element, so its second derivatives are zero
200 if (data_ptr->update_each & update_hessians)
201 output_data.shape_hessians.fill(Tensor<2, 2>());
202
203 return data_ptr;
204}
205
206
207
208void
212 const Quadrature<2> &,
213 const Mapping<2, 2> &,
216 &mapping_data,
217 const FiniteElement<2, 2>::InternalDataBase &fe_internal,
219 const
220{
221 const UpdateFlags flags(fe_internal.update_each);
222
223 const unsigned int n_q_points = mapping_data.quadrature_points.size();
224
225 // linear shape functions
227
228 // compute on the cell
229 if (flags & update_values)
230 for (unsigned int i = 0; i < n_q_points; ++i)
231 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
232 output_data.shape_values[k][i] =
233 (coeffs[k][0] * mapping_data.quadrature_points[i][0] +
234 coeffs[k][1] * mapping_data.quadrature_points[i][1] + coeffs[k][2]);
235
236 if (flags & update_gradients)
237 for (unsigned int i = 0; i < n_q_points; ++i)
238 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
239 output_data.shape_gradients[k][i] =
240 Point<2>(coeffs[k][0], coeffs[k][1]);
241}
242
243
244
245void
248 const unsigned int face_no,
249 const hp::QCollection<1> &quadrature,
250 const Mapping<2, 2> &mapping,
253 const InternalDataBase &fe_internal,
255 &output_data) const
256{
257 AssertDimension(quadrature.size(), 1);
258
259 const UpdateFlags flags(fe_internal.update_each);
260
261 // linear shape functions
263
264 // compute on the face
265 const Quadrature<2> quadrature_on_face =
267 quadrature[0],
268 face_no);
269
270 if (flags & update_values)
271 for (unsigned int i = 0; i < quadrature_on_face.size(); ++i)
272 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
273 {
274 const Point<2> quadrature_point =
275 mapping.transform_unit_to_real_cell(cell,
276 quadrature_on_face.point(i));
277
278 output_data.shape_values[k][i] =
279 (coeffs[k][0] * quadrature_point[0] +
280 coeffs[k][1] * quadrature_point[1] + coeffs[k][2]);
281 }
282
283 if (flags & update_gradients)
284 for (unsigned int i = 0; i < quadrature_on_face.size(); ++i)
285 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
286 output_data.shape_gradients[k][i] =
287 Point<2>(coeffs[k][0], coeffs[k][1]);
288}
289
290
291
292void
295 const unsigned int face_no,
296 const unsigned int sub_no,
297 const Quadrature<1> &quadrature,
298 const Mapping<2, 2> &mapping,
301 const InternalDataBase &fe_internal,
303 &output_data) const
304{
305 const UpdateFlags flags(fe_internal.update_each);
306
307 // linear shape functions
309
310 // compute on the subface
311 const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(
312 this->reference_cell(), quadrature, face_no, sub_no);
313
314 if (flags & update_values)
315 for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i)
316 {
317 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
318 {
319 const Point<2> quadrature_point =
321 cell, quadrature_on_subface.point(i));
322
323 output_data.shape_values[k][i] =
324 (coeffs[k][0] * quadrature_point[0] +
325 coeffs[k][1] * quadrature_point[1] + coeffs[k][2]);
326 }
327 }
328
329 if (flags & update_gradients)
330 for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i)
331 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
332 output_data.shape_gradients[k][i] =
333 Point<2>(coeffs[k][0], coeffs[k][1]);
334}
335
336
337
338void
340{
341 // coefficient relation between children and mother
342 interface_constraints.TableBase<2, double>::reinit(
344
345 interface_constraints(0, 0) = 0.5;
346 interface_constraints(0, 1) = 0.5;
347}
348
349
FE_P1NC()
Definition fe_p1nc.cc:23
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:159
static ndarray< double, 4, 3 > get_linear_shape_coefficients(const Triangulation< 2, 2 >::cell_iterator &cell)
Definition fe_p1nc.cc:88
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:293
virtual std::unique_ptr< FiniteElement< 2, 2 > > clone() const override
Definition fe_p1nc.cc:68
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:246
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:209
void initialize_constraints()
Definition fe_p1nc.cc:339
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:135
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:185
virtual UpdateFlags requires_update_flags(const UpdateFlags flags) const override
Definition fe_p1nc.cc:49
static std::vector< unsigned int > get_dpo_vector()
Definition fe_p1nc.cc:76
virtual std::string get_name() const override
Definition fe_p1nc.cc:41
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:318
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
Definition point.h:111
static void project_to_subface(const ReferenceCell &reference_cell, 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 ReferenceCell &reference_cell, 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
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:83
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define AssertDimension(dim1, dim2)
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.
STL namespace.
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107