Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_p1nc.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2019 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/base/std_cxx14/memory.h>
18 
19 #include <deal.II/fe/fe_p1nc.h>
20 
21 DEAL_II_NAMESPACE_OPEN
22 
23 
25  : FiniteElement<2, 2>(
26  FiniteElementData<2>(get_dpo_vector(), 1, 1, FiniteElementData<2>::L2),
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.resize(2);
32  unit_face_support_points[0][0] = 0.0;
33  unit_face_support_points[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_cell_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_cxx14::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 
88 std::array<std::array<double, 3>, 4>
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  std::array<std::array<double, 3>, 4> 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 =
144  std_cxx14::make_unique<FiniteElement<2, 2>::InternalDataBase>();
145 
146  data_ptr->update_each = requires_update_flags(update_flags);
147 
148  const unsigned int n_q_points = quadrature.size();
149  output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
150 
151  // this is a linear element, so its second derivatives are zero
152  if (data_ptr->update_each & update_hessians)
153  output_data.shape_hessians.fill(Tensor<2, 2>());
154 
155  return data_ptr;
156 }
157 
158 
159 
160 std::unique_ptr<FiniteElement<2, 2>::InternalDataBase>
161 FE_P1NC::get_face_data(
162  const UpdateFlags update_flags,
163  const Mapping<2, 2> &,
164  const Quadrature<1> &quadrature,
166  &output_data) const
167 {
168  auto data_ptr =
169  std_cxx14::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.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>
186 FE_P1NC::get_subface_data(
187  const UpdateFlags update_flags,
188  const Mapping<2, 2> &,
189  const Quadrature<1> &quadrature,
191  &output_data) const
192 {
193  auto data_ptr =
194  std_cxx14::make_unique<FiniteElement<2, 2>::InternalDataBase>();
195 
196  data_ptr->update_each = requires_update_flags(update_flags);
197 
198  const unsigned int n_q_points = quadrature.size();
199  output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each);
200 
201  // this is a linear element, so its second derivatives are zero
202  if (data_ptr->update_each & update_hessians)
203  output_data.shape_hessians.fill(Tensor<2, 2>());
204 
205  return data_ptr;
206 }
207 
208 
209 
210 void
214  const Quadrature<2> &,
215  const Mapping<2, 2> &,
218  & mapping_data,
219  const FiniteElement<2, 2>::InternalDataBase &fe_internal,
221  const
222 {
223  const UpdateFlags flags(fe_internal.update_each);
224 
225  const unsigned int n_q_points = mapping_data.quadrature_points.size();
226 
227  // linear shape functions
228  std::array<std::array<double, 3>, 4> coeffs =
230 
231  // compute on the cell
232  if (flags & update_values)
233  for (unsigned int i = 0; i < n_q_points; ++i)
234  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
235  output_data.shape_values[k][i] =
236  (coeffs[k][0] * mapping_data.quadrature_points[i](0) +
237  coeffs[k][1] * mapping_data.quadrature_points[i](1) + coeffs[k][2]);
238 
239  if (flags & update_gradients)
240  for (unsigned int i = 0; i < n_q_points; ++i)
241  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
242  output_data.shape_gradients[k][i] =
243  Point<2>(coeffs[k][0], coeffs[k][1]);
244 }
245 
246 
247 
248 void
251  const unsigned int face_no,
252  const Quadrature<1> & quadrature,
253  const Mapping<2, 2> & mapping,
255  const ::internal::FEValuesImplementation::MappingRelatedData<2, 2> &,
256  const InternalDataBase &fe_internal,
258  &output_data) const
259 {
260  const UpdateFlags flags(fe_internal.update_each);
261 
262  // linear shape functions
263  const std::array<std::array<double, 3>, 4> coeffs =
265 
266  // compute on the face
267  const Quadrature<2> quadrature_on_face =
268  QProjector<2>::project_to_face(quadrature, 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->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->dofs_per_cell; ++k)
286  output_data.shape_gradients[k][i] =
287  Point<2>(coeffs[k][0], coeffs[k][1]);
288 }
289 
290 
291 
292 void
295  const unsigned int face_no,
296  const unsigned int sub_no,
297  const Quadrature<1> & quadrature,
298  const Mapping<2, 2> & mapping,
300  const ::internal::FEValuesImplementation::MappingRelatedData<2, 2> &,
301  const InternalDataBase &fe_internal,
303  &output_data) const
304 {
305  const UpdateFlags flags(fe_internal.update_each);
306 
307  // linear shape functions
308  const std::array<std::array<double, 3>, 4> coeffs =
310 
311  // compute on the subface
312  const Quadrature<2> quadrature_on_subface =
313  QProjector<2>::project_to_subface(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->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->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 
351 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
Shape function values.
virtual void fill_fe_face_values(const Triangulation< 2, 2 >::cell_iterator &cell, const unsigned int face_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:249
FE_P1NC()
Definition: fe_p1nc.cc:24
FullMatrix< double > interface_constraints
Definition: fe.h:2468
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)
const Point< dim > & point(const unsigned int i) const
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
std::vector< Point< dim - 1 > > unit_face_support_points
Definition: fe.h:2487
STL namespace.
Transformed quadrature points.
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:293
static std::vector< unsigned int > get_dpo_vector()
Definition: fe_p1nc.cc:77
No update.
virtual UpdateFlags requires_update_flags(const UpdateFlags flags) const override
Definition: fe_p1nc.cc:50
static std::array< std::array< double, 3 >, 4 > get_linear_shape_coefficients(const Triangulation< 2, 2 >::cell_iterator &cell)
Definition: fe_p1nc.cc:89
UpdateFlags
Abstract base class for mapping classes.
Definition: dof_tools.h:57
virtual std::string get_name() const override
Definition: fe_p1nc.cc:42
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:3007
Second derivatives of shape functions.
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
std::vector< Point< spacedim > > quadrature_points
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
Definition: mpi.h:90
Shape function gradients.
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:211
void initialize_constraints()
Definition: fe_p1nc.cc:340
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
TableIndices< 2 > interface_constraints_size() const
Definition: fe.cc:857
void fill(InputIterator entries, const bool C_style_indexing=true)