Reference documentation for deal.II version 8.5.1
fe_p1nc.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2016 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/fe/fe_p1nc.h>
18 
19 DEAL_II_NAMESPACE_OPEN
20 
21 
23  :
24  FiniteElement<2,2>(FiniteElementData<2>(get_dpo_vector(),
25  1,
26  1,
27  FiniteElementData<2>::L2),
28  std::vector<bool> (1, false),
29  std::vector<ComponentMask>(4, ComponentMask(1,true)))
30 {
31  // face support points: 2 end vertices
32  unit_face_support_points.resize(2);
33  unit_face_support_points[0][0] = 0.0;
34  unit_face_support_points[1][0] = 1.0;
35 
36  // initialize constraints matrix
38 }
39 
40 
41 
42 std::string FE_P1NC::get_name () const
43 {
44  return "FE_P1NC";
45 }
46 
47 
48 
50 {
52 
53  if (flags & update_values)
55  if (flags & update_gradients)
56  out |= update_gradients;
57  if (flags & update_cell_normal_vectors)
59  if (flags & update_hessians)
60  out |= update_hessians;
61 
62  return out;
63 }
64 
65 
66 
68 {
69  return new FE_P1NC(*this);
70 }
71 
72 
73 
75 
76 
77 
78 std::vector<unsigned int>
80 {
81  std::vector<unsigned int> dpo(3);
82  dpo[0] = 1; // dofs per object: vertex
83  dpo[1] = 0; // line
84  dpo[2] = 0; // quad
85  return dpo;
86 }
87 
88 
89 
90 std_cxx11::array<std_cxx11::array<double,3>,4>
92 {
93  // edge midpoints
94  const Point<2> mpt[4] = { (cell->vertex(0) + cell->vertex(2)) / 2,
95  (cell->vertex(1) + cell->vertex(3)) / 2,
96  (cell->vertex(0) + cell->vertex(1)) / 2,
97  (cell->vertex(2) + cell->vertex(3)) / 2
98  };
99 
100  // center point
101  const Point<2> cpt = (cell->vertex(0) +
102  cell->vertex(1) +
103  cell->vertex(2) +
104  cell->vertex(3)) / 4;
105 
106  const double det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1));
107 
108  std_cxx11::array<std_cxx11::array<double,3>,4> coeffs;
109  coeffs[0][0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det;
110  coeffs[1][0] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det;
111  coeffs[2][0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det;
112  coeffs[3][0] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det;
113 
114  coeffs[0][1] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det;
115  coeffs[1][1] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det;
116  coeffs[2][1] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det;
117  coeffs[3][1] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det;
118 
119  coeffs[0][2] = 0.25 - cpt(0)*coeffs[0][0] - cpt(1)*coeffs[0][1];
120  coeffs[1][2] = 0.25 - cpt(0)*coeffs[1][0] - cpt(1)*coeffs[1][1];
121  coeffs[2][2] = 0.25 - cpt(0)*coeffs[2][0] - cpt(1)*coeffs[2][1];
122  coeffs[3][2] = 0.25 - cpt(0)*coeffs[3][0] - cpt(1)*coeffs[3][1];
123 
124  return coeffs;
125 }
126 
127 
128 
130 FE_P1NC::get_data (const UpdateFlags update_flags,
131  const Mapping<2,2> &,
132  const Quadrature<2> &quadrature,
134 {
136 
137  data->update_each = requires_update_flags(update_flags);
138 
139  const unsigned int n_q_points = quadrature.size();
140  output_data.initialize (n_q_points, FE_P1NC(), data->update_each);
141 
142  // this is a linear element, so its second derivatives are zero
143  if (data->update_each & update_hessians)
144  output_data.shape_hessians.fill(Tensor<2,2>());
145 
146  return data;
147 }
148 
149 
150 
152 FE_P1NC::get_face_data (const UpdateFlags update_flags,
153  const Mapping<2,2> &,
154  const Quadrature<1> &quadrature,
156 {
158 
159  data->update_each = requires_update_flags(update_flags);
160 
161  const unsigned int n_q_points = quadrature.size();
162  output_data.initialize (n_q_points, FE_P1NC(), data->update_each);
163 
164  // this is a linear element, so its second derivatives are zero
165  if (data->update_each & update_hessians)
166  output_data.shape_hessians.fill(Tensor<2,2>());
167 
168  return data;
169 }
170 
171 
172 
174 FE_P1NC::get_subface_data (const UpdateFlags update_flags,
175  const Mapping<2,2> &,
176  const Quadrature<1> &quadrature,
178 {
180 
181  data->update_each = requires_update_flags(update_flags);
182 
183  const unsigned int n_q_points = quadrature.size();
184  output_data.initialize (n_q_points, FE_P1NC(), data->update_each);
185 
186  // this is a linear element, so its second derivatives are zero
187  if (data->update_each & update_hessians)
188  output_data.shape_hessians.fill(Tensor<2,2>());
189 
190  return data;
191 }
192 
193 
194 
195 void
198  const Quadrature<2> &,
199  const Mapping<2,2> &,
202  const FiniteElement<2,2>::InternalDataBase &fe_internal,
204 {
205  const UpdateFlags flags(fe_internal.update_each);
206 
207  const unsigned int n_q_points = mapping_data.quadrature_points.size();
208 
209  // linear shape functions
210  std_cxx11::array<std_cxx11::array<double,3>,4> coeffs = get_linear_shape_coefficients (cell);
211 
212  // compute on the cell
213  if (flags & update_values)
214  for (unsigned int i=0; i<n_q_points; ++i)
215  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
216  output_data.shape_values[k][i] = (coeffs[k][0]*mapping_data.quadrature_points[i](0) +
217  coeffs[k][1]*mapping_data.quadrature_points[i](1) +
218  coeffs[k][2]);
219 
220  if (flags & update_gradients)
221  for (unsigned int i=0; i<n_q_points; ++i)
222  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
223  output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0],
224  coeffs[k][1]);
225 }
226 
227 
228 
229 void
231  const unsigned int face_no,
232  const Quadrature<1> &quadrature,
233  const Mapping<2,2> &mapping,
235  const ::internal::FEValues::MappingRelatedData<2,2> &,
236  const InternalDataBase &fe_internal,
238 {
239  const UpdateFlags flags(fe_internal.update_each);
240 
241  // linear shape functions
242  const std_cxx11::array<std_cxx11::array<double,3>,4> coeffs
244 
245  // compute on the face
246  const Quadrature<2> quadrature_on_face = QProjector<2>::project_to_face(quadrature, face_no);
247 
248  if (flags & update_values)
249  for (unsigned int i=0; i<quadrature_on_face.size(); ++i)
250  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
251  {
252  const Point<2> quadrature_point
253  = mapping.transform_unit_to_real_cell(cell, quadrature_on_face.point(i));
254 
255  output_data.shape_values[k][i] = (coeffs[k][0]*quadrature_point(0) +
256  coeffs[k][1]*quadrature_point(1) +
257  coeffs[k][2]);
258  }
259 
260  if (flags & update_gradients)
261  for (unsigned int i=0; i<quadrature_on_face.size(); ++i)
262  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
263  output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0],
264  coeffs[k][1]);
265 }
266 
267 
268 
269 void
271  const unsigned int face_no,
272  const unsigned int sub_no,
273  const Quadrature<1> &quadrature,
274  const Mapping<2,2> &mapping,
276  const ::internal::FEValues::MappingRelatedData<2,2> &,
277  const InternalDataBase &fe_internal,
279 {
280  const UpdateFlags flags(fe_internal.update_each);
281 
282  // linear shape functions
283  const std_cxx11::array<std_cxx11::array<double,3>,4> coeffs
285 
286  // compute on the subface
287  const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(quadrature, face_no, sub_no);
288 
289  if (flags & update_values)
290  for (unsigned int i=0; i<quadrature_on_subface.size(); ++i)
291  {
292  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
293  {
294  const Point<2> quadrature_point
295  = mapping.transform_unit_to_real_cell(cell, quadrature_on_subface.point(i));
296 
297  output_data.shape_values[k][i] = (coeffs[k][0]*quadrature_point(0) +
298  coeffs[k][1]*quadrature_point(1) +
299  coeffs[k][2]);
300  }
301  }
302 
303  if (flags & update_gradients)
304  for (unsigned int i=0; i<quadrature_on_subface.size(); ++i)
305  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
306  output_data.shape_gradients[k][i] = Point<2>(coeffs[k][0],
307  coeffs[k][1]);
308 }
309 
310 
311 
313 {
314  // coefficient relation between children and mother
316  .TableBase<2,double>::reinit (interface_constraints_size());
317 
318  interface_constraints(0,0) = 0.5;
319  interface_constraints(0,1) = 0.5;
320 }
321 
322 
323 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
Shape function values.
FE_P1NC()
Definition: fe_p1nc.cc:22
FullMatrix< double > interface_constraints
Definition: fe.h:2211
virtual FiniteElement< 2, 2 >::InternalDataBase * get_data(const UpdateFlags update_flags, const Mapping< 2, 2 > &, const Quadrature< 2 > &quadrature, ::internal::FEValues::FiniteElementRelatedData< 2, 2 > &output_data) const
Definition: fe_p1nc.cc:130
const Point< dim > & point(const unsigned int i) const
STL namespace.
Transformed quadrature points.
void initialize(const unsigned int n_quadrature_points, const FiniteElement< dim, spacedim > &fe, const UpdateFlags flags)
Definition: fe_values.cc:2181
virtual ~FE_P1NC()
Definition: fe_p1nc.cc:74
static std_cxx11::array< std_cxx11::array< double, 3 >, 4 > get_linear_shape_coefficients(const Triangulation< 2, 2 >::cell_iterator &cell)
Definition: fe_p1nc.cc:91
static std::vector< unsigned int > get_dpo_vector()
Definition: fe_p1nc.cc:79
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::FEValues::MappingRelatedData< 2, 2 > &mapping_data, const FiniteElement< 2, 2 >::InternalDataBase &fe_internal, internal::FEValues::FiniteElementRelatedData< 2, 2 > &output_data) const
Definition: fe_p1nc.cc:196
No update.
UpdateFlags
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)
virtual UpdateFlags requires_update_flags(const UpdateFlags flags) const
Definition: fe_p1nc.cc:49
Abstract base class for mapping classes.
Definition: dof_tools.h:46
std::vector< Point< dim-1 > > unit_face_support_points
Definition: fe.h:2230
Second derivatives of shape functions.
virtual FiniteElement< 2, 2 > * clone() const
Definition: fe_p1nc.cc:67
virtual std::string get_name() const
Definition: fe_p1nc.cc:42
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:295
Definition: mpi.h:41
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::FEValues::MappingRelatedData< 2, 2 > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValues::FiniteElementRelatedData< 2, 2 > &output_data) const
Definition: fe_p1nc.cc:230
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
Shape function gradients.
std::vector< Point< spacedim > > quadrature_points
void initialize_constraints()
Definition: fe_p1nc.cc:312
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:823
void fill(InputIterator entries, const bool C_style_indexing=true)
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::FEValues::MappingRelatedData< 2, 2 > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValues::FiniteElementRelatedData< 2, 2 > &output_data) const
Definition: fe_p1nc.cc:270