Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55: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_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 
41 std::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 
68 std::unique_ptr<FiniteElement<2, 2>>
70 {
71  return std::make_unique<FE_P1NC>(*this);
72 }
73 
74 
75 
76 std::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 
105  ndarray<double, 4, 3> coeffs;
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 
135 std::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 
159 std::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 
185 std::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 
209 void
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 
246 void
249  const unsigned int face_no,
250  const hp::QCollection<1> & quadrature,
251  const Mapping<2, 2> & mapping,
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 
293 void
296  const unsigned int face_no,
297  const unsigned int sub_no,
298  const Quadrature<1> & quadrature,
299  const Mapping<2, 2> & mapping,
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 
339 void
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 
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 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 std::unique_ptr< FiniteElement< 2, 2 > > clone() const override
Definition: fe_p1nc.cc:69
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
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:2447
TableIndices< 2 > interface_constraints_size() const
FullMatrix< double > interface_constraints
Definition: fe.h:2428
Abstract base class for mapping classes.
Definition: mapping.h:314
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 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
Definition: tensor.h:516
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:2702
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
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.
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
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108