Reference documentation for deal.II version 9.3.2
\(\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_nedelec_sz.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2021 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_nedelec_sz.h>
18 
19 #include <memory>
20 
21 DEAL_II_NAMESPACE_OPEN
22 
23 // Constructor:
24 template <int dim, int spacedim>
25 FE_NedelecSZ<dim, spacedim>::FE_NedelecSZ(const unsigned int order)
26  : FiniteElement<dim, dim>(
27  FiniteElementData<dim>(get_dpo_vector(order),
28  dim,
29  order + 1,
30  FiniteElementData<dim>::Hcurl),
31  std::vector<bool>(compute_num_dofs(order), true),
32  std::vector<ComponentMask>(compute_num_dofs(order),
33  std::vector<bool>(dim, true)))
34 {
35  Assert(dim >= 2, ExcImpossibleInDim(dim));
36 
37  this->mapping_kind = mapping_nedelec;
38  // Set up the table converting components to base components. Since we have
39  // only one base element, everything remains zero except the component in the
40  // base, which is the component itself.
41  for (unsigned int comp = 0; comp < this->n_components(); ++comp)
42  {
43  this->component_to_base_table[comp].first.second = comp;
44  }
45 
46  // Generate the 1-D polynomial basis.
47  create_polynomials(order);
48 }
49 
50 
51 
52 // Shape functions:
53 template <int dim, int spacedim>
54 double
55 FE_NedelecSZ<dim, spacedim>::shape_value(const unsigned int /*i*/,
56  const Point<dim> & /*p*/) const
57 {
58  Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
59  return 0.;
60 }
61 
62 
63 
64 template <int dim, int spacedim>
65 double
66 FE_NedelecSZ<dim, spacedim>::shape_value_component(
67  const unsigned int /*i*/,
68  const Point<dim> & /*p*/,
69  const unsigned int /*component*/) const
70 {
71  // Not implemented yet:
72  Assert(false, ExcNotImplemented());
73  return 0.;
74 }
75 
76 
77 
78 template <int dim, int spacedim>
79 Tensor<1, dim>
80 FE_NedelecSZ<dim, spacedim>::shape_grad(const unsigned int /*i*/,
81  const Point<dim> & /*p*/) const
82 {
83  Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
84  return Tensor<1, dim>();
85 }
86 
87 
88 
89 template <int dim, int spacedim>
90 Tensor<1, dim>
91 FE_NedelecSZ<dim, spacedim>::shape_grad_component(
92  const unsigned int /*i*/,
93  const Point<dim> & /*p*/,
94  const unsigned int /*component*/) const
95 {
96  Assert(false, ExcNotImplemented());
97  return Tensor<1, dim>();
98 }
99 
100 
101 
102 template <int dim, int spacedim>
103 Tensor<2, dim>
104 FE_NedelecSZ<dim, spacedim>::shape_grad_grad(const unsigned int /*i*/,
105  const Point<dim> & /*p*/) const
106 {
107  Assert(false, (typename FiniteElement<dim, spacedim>::ExcFENotPrimitive()));
108  return Tensor<2, dim>();
109 }
110 
111 
112 
113 template <int dim, int spacedim>
114 Tensor<2, dim>
115 FE_NedelecSZ<dim, spacedim>::shape_grad_grad_component(
116  const unsigned int /*i*/,
117  const Point<dim> & /*p*/,
118  const unsigned int /*component*/) const
119 {
120  Assert(false, ExcNotImplemented());
121  return Tensor<2, dim>();
122 }
123 
124 
125 
126 template <int dim, int spacedim>
127 std::unique_ptr<typename ::FiniteElement<dim, spacedim>::InternalDataBase>
128 FE_NedelecSZ<dim, spacedim>::get_data(
129  const UpdateFlags update_flags,
130  const Mapping<dim, spacedim> & /*mapping*/,
131  const Quadrature<dim> &quadrature,
132  ::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
133  spacedim>
134  & /*output_data*/) const
135 {
136  std::unique_ptr<
137  typename ::FiniteElement<dim, spacedim>::InternalDataBase>
138  data_ptr = std::make_unique<InternalData>();
139  auto &data = dynamic_cast<InternalData &>(*data_ptr);
140  data.update_each = requires_update_flags(update_flags);
141 
142  // Useful quantities:
143  const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
144 
145  const unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
146  const unsigned int lines_per_cell = GeometryInfo<dim>::lines_per_cell;
147  const unsigned int faces_per_cell = GeometryInfo<dim>::faces_per_cell;
148 
149  const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
150 
151  // we assume that all quads have the same number of dofs
152  const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
153 
154  const UpdateFlags flags(data.update_each);
155  const unsigned int n_q_points = quadrature.size();
156 
157  // Resize the internal data storage:
158  data.sigma_imj_values.resize(
159  n_q_points,
160  std::vector<std::vector<double>>(vertices_per_cell,
161  std::vector<double>(vertices_per_cell)));
162 
163  data.sigma_imj_grads.resize(vertices_per_cell,
164  std::vector<std::vector<double>>(
165  vertices_per_cell, std::vector<double>(dim)));
166 
167  // Resize shape function arrays according to update flags:
168  if (flags & update_values)
169  {
170  data.shape_values.resize(this->n_dofs_per_cell(),
171  std::vector<Tensor<1, dim>>(n_q_points));
172  }
173 
174  if (flags & update_gradients)
175  {
176  data.shape_grads.resize(this->n_dofs_per_cell(),
177  std::vector<DerivativeForm<1, dim, dim>>(
178  n_q_points));
179  }
180  // Not implementing second derivatives yet:
181  if (flags & update_hessians)
182  {
183  Assert(false, ExcNotImplemented());
184  }
185 
186  std::vector<Point<dim>> p_list(n_q_points);
187  p_list = quadrature.get_points();
188 
189 
190  switch (dim)
191  {
192  case 2:
193  {
194  // Compute values of sigma & lambda and the sigma differences and
195  // lambda additions.
196  std::vector<std::vector<double>> sigma(
197  n_q_points, std::vector<double>(lines_per_cell));
198  std::vector<std::vector<double>> lambda(
199  n_q_points, std::vector<double>(lines_per_cell));
200 
201  for (unsigned int q = 0; q < n_q_points; ++q)
202  {
203  sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
204  sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
205  sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
206  sigma[q][3] = p_list[q][0] + p_list[q][1];
207 
208  lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
209  lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
210  lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
211  lambda[q][3] = p_list[q][0] * p_list[q][1];
212  for (unsigned int i = 0; i < vertices_per_cell; ++i)
213  {
214  for (unsigned int j = 0; j < vertices_per_cell; ++j)
215  {
216  data.sigma_imj_values[q][i][j] =
217  sigma[q][i] - sigma[q][j];
218  }
219  }
220  }
221 
222  // Calculate the gradient of sigma_imj_values[q][i][j] =
223  // sigma[q][i]-sigma[q][j]
224  // - this depends on the component and the direction of the
225  // corresponding edge.
226  // - the direction of the edge is determined by
227  // sigma_imj_sign[i][j].
228  // Helper arrays:
229  const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2] = {
230  {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
231  int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
232  unsigned int sigma_imj_component[vertices_per_cell]
233  [vertices_per_cell];
234 
235  for (unsigned int i = 0; i < vertices_per_cell; ++i)
236  {
237  for (unsigned int j = 0; j < vertices_per_cell; ++j)
238  {
239  // sigma_imj_sign is the sign (+/-) of the coefficient of
240  // x/y/z in sigma_imj_values Due to the numbering of vertices
241  // on the reference element it is easy to find edges in the
242  // positive direction are from smaller to higher local vertex
243  // numbering.
244  sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
245  sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
246 
247  // Now store the component which the sigma_i - sigma_j
248  // corresponds to:
249  sigma_imj_component[i][j] = 0;
250  for (unsigned int d = 0; d < dim; ++d)
251  {
252  int temp_imj =
253  sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
254  // Only interested in the first non-zero
255  // as if there is a second, it can not be a valid edge.
256  if (temp_imj != 0)
257  {
258  sigma_imj_component[i][j] = d;
259  break;
260  }
261  }
262  // Can now calculate the gradient, only non-zero in the
263  // component given: Note some i,j combinations will be
264  // incorrect, but only on invalid edges.
265  data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
266  2.0 * sigma_imj_sign[i][j];
267  }
268  }
269 
270  // Now compute the edge parameterisations for a single element
271  // with global numbering matching that of the reference element:
272 
273  // Resize the edge parameterisations
274  data.edge_sigma_values.resize(lines_per_cell);
275  data.edge_sigma_grads.resize(lines_per_cell);
276  for (unsigned int m = 0; m < lines_per_cell; ++m)
277  {
278  data.edge_sigma_values[m].resize(n_q_points);
279 
280  // sigma grads are constant in a cell (no need for quad points)
281  data.edge_sigma_grads[m].resize(dim);
282  }
283 
284  // Fill the values for edge lambda and edge sigma:
285  const unsigned int
286  edge_sigma_direction[GeometryInfo<2>::lines_per_cell] = {1,
287  1,
288  0,
289  0};
290 
291  data.edge_lambda_values.resize(lines_per_cell,
292  std::vector<double>(n_q_points));
293  data.edge_lambda_grads_2d.resize(lines_per_cell,
294  std::vector<double>(dim));
295  for (unsigned int m = 0; m < lines_per_cell; ++m)
296  {
297  // e1=max(reference vertex numbering on this edge)
298  // e2=min(reference vertex numbering on this edge)
299  // Which is guaranteed to be:
300  const unsigned int e1(
301  GeometryInfo<dim>::line_to_cell_vertices(m, 1));
302  const unsigned int e2(
303  GeometryInfo<dim>::line_to_cell_vertices(m, 0));
304  for (unsigned int q = 0; q < n_q_points; ++q)
305  {
306  data.edge_sigma_values[m][q] =
307  data.sigma_imj_values[q][e2][e1];
308  data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
309  }
310 
311  data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
312  }
313  data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
314  data.edge_lambda_grads_2d[1] = {1.0, 0.0};
315  data.edge_lambda_grads_2d[2] = {0.0, -1.0};
316  data.edge_lambda_grads_2d[3] = {0.0, 1.0};
317 
318  // If the polynomial order is 0, then no more work to do:
319  if (degree < 1)
320  {
321  break;
322  }
323 
324  // Otherwise, we can compute the non-cell dependent shape functions.
325  //
326  // Note: the local dof numberings follow the usual order of lines ->
327  // faces -> cells
328  // (we have no vertex-based DoFs in this element).
329  // For a given cell we have:
330  // n_line_dofs = dofs_per_line*lines_per_cell.
331  // n_face_dofs = dofs_per_face*faces_per_cell.
332  // n_cell_dofs = dofs_per_quad (2D)
333  // = dofs_per_hex (3D)
334  //
335  // i.e. For the local dof numbering:
336  // the first line dof is 0,
337  // the first face dof is n_line_dofs,
338  // the first cell dof is n_line_dofs + n_face_dofs.
339  //
340  // On a line, DoFs are ordered first by line_dof and then line_index:
341  // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
342  //
343  // and similarly for faces:
344  // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
345  //
346  // HOWEVER, we have different types of DoFs on a line/face/cell.
347  // On a line we have two types, lowest order and higher order
348  // gradients.
349  // - The numbering is such the lowest order is first, then higher
350  // order.
351  // This is simple enough as there is only 1 lowest order and
352  // degree higher orders DoFs per line.
353  //
354  // On a 2D cell, we have 3 types: Type 1/2/3:
355  // - The ordering done by type:
356  // - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
357  // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index
358  // = ij1.
359  // - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
360  // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index
361  // = degree^2 + ij2
362  // - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
363  // Numbered: ij3 = i3. i.e. cell_dof_index
364  // = 2*(degree^2) + ij3.
365  //
366  // These then fit into the local dof numbering described above:
367  // - local dof numberings are:
368  // line_dofs: local_dof = line_dof_index. 0 <= local_dof <
369  // dofs_per_line*lines_per_cell face_dofs: local_dof =
370  // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
371  // = n_lines_dof + n_face_dofs + cell_dof_index.
372  //
373  // The cell-based shape functions are:
374  //
375  // Type 1 (gradients):
376  // \phi^{C_{1}}_{ij) = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
377  //
378  // 0 <= i,j < degree.
379  //
380  // NOTE: The derivative produced by IntegratedLegendrePolynomials does
381  // not account for the
382  // (2*x-1) or (2*y-1) so we must take this into account when
383  // taking derivatives.
384  const unsigned int cell_type1_offset = n_line_dofs;
385 
386  // Type 2:
387  // \phi^{C_{2}}_{ij) = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
388  // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
389  //
390  // 0 <= i,j < degree.
391  const unsigned int cell_type2_offset =
392  cell_type1_offset + degree * degree;
393 
394  // Type 3 (two subtypes):
395  // \phi^{C_{3}}_{j) = L_{j+2}(2y-1) \mathbf{e}_{x}
396  //
397  // \phi^{C_{3}}_{j+degree) = L_{j+2}(2x-1) \mathbf{e}_{y}
398  //
399  // 0 <= j < degree
400  const unsigned int cell_type3_offset1 =
401  cell_type2_offset + degree * degree;
402  const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
403 
404  if (flags & (update_values | update_gradients))
405  {
406  // compute all points we must evaluate the 1d polynomials at:
407  std::vector<Point<dim>> cell_points(n_q_points);
408  for (unsigned int q = 0; q < n_q_points; ++q)
409  {
410  for (unsigned int d = 0; d < dim; ++d)
411  {
412  cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
413  }
414  }
415 
416  // Loop through quad points:
417  for (unsigned int q = 0; q < n_q_points; ++q)
418  {
419  // pre-compute values & required derivatives at this quad
420  // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
421  //
422  // for each polyc[d], c=x,y, contains the d-th derivative with
423  // respect to the coordinate c.
424 
425  // We only need poly values and 1st derivative for
426  // update_values, but need the 2nd derivative too for
427  // update_gradients.
428  //
429  // Note that this will need to be updated if we're supporting
430  // update_hessians.
431  const unsigned int poly_length(
432  (flags & update_gradients) ? 3 : 2);
433 
434  std::vector<std::vector<double>> polyx(
435  degree, std::vector<double>(poly_length));
436  std::vector<std::vector<double>> polyy(
437  degree, std::vector<double>(poly_length));
438  for (unsigned int i = 0; i < degree; ++i)
439  {
440  // Compute all required 1d polynomials and their
441  // derivatives, starting at degree 2. e.g. to access
442  // L'_{3}(2x-1) use polyx[1][1].
443  IntegratedLegendrePolynomials[i + 2].value(
444  cell_points[q][0], polyx[i]);
445  IntegratedLegendrePolynomials[i + 2].value(
446  cell_points[q][1], polyy[i]);
447  }
448  // Now use these to compute the shape functions:
449  if (flags & update_values)
450  {
451  for (unsigned int j = 0; j < degree; ++j)
452  {
453  const unsigned int shift_j(j * degree);
454  for (unsigned int i = 0; i < degree; ++i)
455  {
456  const unsigned int shift_ij(i + shift_j);
457 
458  // Type 1:
459  const unsigned int dof_index1(cell_type1_offset +
460  shift_ij);
461  data.shape_values[dof_index1][q][0] =
462  2.0 * polyx[i][1] * polyy[j][0];
463  data.shape_values[dof_index1][q][1] =
464  2.0 * polyx[i][0] * polyy[j][1];
465 
466  // Type 2:
467  const unsigned int dof_index2(cell_type2_offset +
468  shift_ij);
469  data.shape_values[dof_index2][q][0] =
470  data.shape_values[dof_index1][q][0];
471  data.shape_values[dof_index2][q][1] =
472  -1.0 * data.shape_values[dof_index1][q][1];
473  }
474  // Type 3:
475  const unsigned int dof_index3_1(cell_type3_offset1 +
476  j);
477  data.shape_values[dof_index3_1][q][0] = polyy[j][0];
478  data.shape_values[dof_index3_1][q][1] = 0.0;
479 
480  const unsigned int dof_index3_2(cell_type3_offset2 +
481  j);
482  data.shape_values[dof_index3_2][q][0] = 0.0;
483  data.shape_values[dof_index3_2][q][1] = polyx[j][0];
484  }
485  }
486  if (flags & update_gradients)
487  {
488  for (unsigned int j = 0; j < degree; ++j)
489  {
490  const unsigned int shift_j(j * degree);
491  for (unsigned int i = 0; i < degree; ++i)
492  {
493  const unsigned int shift_ij(i + shift_j);
494 
495  // Type 1:
496  const unsigned int dof_index1(cell_type1_offset +
497  shift_ij);
498  data.shape_grads[dof_index1][q][0][0] =
499  4.0 * polyx[i][2] * polyy[j][0];
500  data.shape_grads[dof_index1][q][0][1] =
501  4.0 * polyx[i][1] * polyy[j][1];
502  data.shape_grads[dof_index1][q][1][0] =
503  data.shape_grads[dof_index1][q][0][1];
504  data.shape_grads[dof_index1][q][1][1] =
505  4.0 * polyx[i][0] * polyy[j][2];
506 
507  // Type 2:
508  const unsigned int dof_index2(cell_type2_offset +
509  shift_ij);
510  data.shape_grads[dof_index2][q][0][0] =
511  data.shape_grads[dof_index1][q][0][0];
512  data.shape_grads[dof_index2][q][0][1] =
513  data.shape_grads[dof_index1][q][0][1];
514  data.shape_grads[dof_index2][q][1][0] =
515  -1.0 * data.shape_grads[dof_index1][q][1][0];
516  data.shape_grads[dof_index2][q][1][1] =
517  -1.0 * data.shape_grads[dof_index1][q][1][1];
518  }
519  // Type 3:
520  const unsigned int dof_index3_1(cell_type3_offset1 +
521  j);
522  data.shape_grads[dof_index3_1][q][0][0] = 0.0;
523  data.shape_grads[dof_index3_1][q][0][1] =
524  2.0 * polyy[j][1];
525  data.shape_grads[dof_index3_1][q][1][0] = 0.0;
526  data.shape_grads[dof_index3_1][q][1][1] = 0.0;
527 
528  const unsigned int dof_index3_2(cell_type3_offset2 +
529  j);
530  data.shape_grads[dof_index3_2][q][0][0] = 0.0;
531  data.shape_grads[dof_index3_2][q][0][1] = 0.0;
532  data.shape_grads[dof_index3_2][q][1][0] =
533  2.0 * polyx[j][1];
534  data.shape_grads[dof_index3_2][q][1][1] = 0.0;
535  }
536  }
537  }
538  }
539  break;
540  }
541  case 3:
542  {
543  // Compute values of sigma & lambda and the sigma differences and
544  // lambda additions.
545  std::vector<std::vector<double>> sigma(
546  n_q_points, std::vector<double>(lines_per_cell));
547  std::vector<std::vector<double>> lambda(
548  n_q_points, std::vector<double>(lines_per_cell));
549  for (unsigned int q = 0; q < n_q_points; ++q)
550  {
551  sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
552  (1 - p_list[q][2]);
553  sigma[q][1] =
554  p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
555  sigma[q][2] =
556  (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
557  sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
558  sigma[q][4] =
559  (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
560  sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
561  sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
562  sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
563 
564  lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
565  (1.0 - p_list[q][2]);
566  lambda[q][1] =
567  p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
568  lambda[q][2] =
569  (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
570  lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
571  lambda[q][4] =
572  (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
573  lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
574  lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
575  lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
576 
577  // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
578  // and lambda_ipj = \lambda_{i} + \lambda_{j}.
579  for (unsigned int i = 0; i < vertices_per_cell; ++i)
580  {
581  for (unsigned int j = 0; j < vertices_per_cell; ++j)
582  {
583  data.sigma_imj_values[q][i][j] =
584  sigma[q][i] - sigma[q][j];
585  }
586  }
587  }
588 
589  // We now want some additional information about
590  // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
591  // calculate values & derivatives of the shape functions we need to
592  // know:
593  // - The component the sigma_imj value corresponds to - this varies
594  // with i & j.
595  // - The gradient of the sigma_imj value
596  // - this depends on the component and the direction of the
597  // corresponding edge.
598  // - the direction of the edge is determined by
599  // sigma_imj_sign[i][j].
600  //
601  // Note that not every i,j combination is a valid edge (there are only
602  // 12 valid edges in 3D), but we compute them all as it simplifies
603  // things.
604 
605  // store the sign of each component x, y, z in the sigma list.
606  // can use this to fill in the sigma_imj_component data.
607  const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3] = {
608  {-1, -1, -1},
609  {1, -1, -1},
610  {-1, 1, -1},
611  {1, 1, -1},
612  {-1, -1, 1},
613  {1, -1, 1},
614  {-1, 1, 1},
615  {1, 1, 1}};
616 
617  int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
618  unsigned int sigma_imj_component[vertices_per_cell]
619  [vertices_per_cell];
620 
621  for (unsigned int i = 0; i < vertices_per_cell; ++i)
622  {
623  for (unsigned int j = 0; j < vertices_per_cell; ++j)
624  {
625  // sigma_imj_sign is the sign (+/-) of the coefficient of
626  // x/y/z in sigma_imj. Due to the numbering of vertices on the
627  // reference element this is easy to work out because edges in
628  // the positive direction go from smaller to higher local
629  // vertex numbering.
630  sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
631  sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
632 
633  // Now store the component which the sigma_i - sigma_j
634  // corresponds to:
635  sigma_imj_component[i][j] = 0;
636  for (unsigned int d = 0; d < dim; ++d)
637  {
638  int temp_imj =
639  sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
640  // Only interested in the first non-zero
641  // as if there is a second, it will not be a valid edge.
642  if (temp_imj != 0)
643  {
644  sigma_imj_component[i][j] = d;
645  break;
646  }
647  }
648  // Can now calculate the gradient, only non-zero in the
649  // component given: Note some i,j combinations will be
650  // incorrect, but only on invalid edges.
651  data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
652  2.0 * sigma_imj_sign[i][j];
653  }
654  }
655  // Now compute the edge parameterisations for a single element
656  // with global numbering matching that of the reference element:
657 
658  // resize the edge parameterisations
659  data.edge_sigma_values.resize(lines_per_cell);
660  data.edge_lambda_values.resize(lines_per_cell);
661  data.edge_sigma_grads.resize(lines_per_cell);
662  data.edge_lambda_grads_3d.resize(lines_per_cell);
663  data.edge_lambda_gradgrads_3d.resize(lines_per_cell);
664  for (unsigned int m = 0; m < lines_per_cell; ++m)
665  {
666  data.edge_sigma_values[m].resize(n_q_points);
667  data.edge_lambda_values[m].resize(n_q_points);
668 
669  // sigma grads are constant in a cell (no need for quad points)
670  data.edge_sigma_grads[m].resize(dim);
671 
672  data.edge_lambda_grads_3d[m].resize(n_q_points);
673  for (unsigned int q = 0; q < n_q_points; ++q)
674  {
675  data.edge_lambda_grads_3d[m][q].resize(dim);
676  }
677  // lambda_gradgrads are constant in a cell (no need for quad
678  // points)
679  data.edge_lambda_gradgrads_3d[m].resize(dim);
680  for (unsigned int d = 0; d < dim; ++d)
681  {
682  data.edge_lambda_gradgrads_3d[m][d].resize(dim);
683  }
684  }
685 
686  // Fill the values:
687  const unsigned int
688  edge_sigma_direction[GeometryInfo<3>::lines_per_cell] = {
689  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
690 
691  for (unsigned int m = 0; m < lines_per_cell; ++m)
692  {
693  // e1=max(reference vertex numbering on this edge)
694  // e2=min(reference vertex numbering on this edge)
695  // Which is guaranteed to be:
696  const unsigned int e1(
697  GeometryInfo<dim>::line_to_cell_vertices(m, 1));
698  const unsigned int e2(
699  GeometryInfo<dim>::line_to_cell_vertices(m, 0));
700 
701  for (unsigned int q = 0; q < n_q_points; ++q)
702  {
703  data.edge_sigma_values[m][q] =
704  data.sigma_imj_values[q][e2][e1];
705  data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
706  }
707 
708  data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
709  }
710  // edge_lambda_grads
711  for (unsigned int q = 0; q < n_q_points; ++q)
712  {
713  double x(p_list[q][0]);
714  double y(p_list[q][1]);
715  double z(p_list[q][2]);
716  data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
717  data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
718  data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
719  data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
720  data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
721  data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
722  data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
723  data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
724  data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
725  data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
726  data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
727  data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
728  }
729  // edge_lambda gradgrads:
730  const int edge_lambda_sign[GeometryInfo<3>::lines_per_cell] = {
731  1,
732  -1,
733  1,
734  -1,
735  -1,
736  1,
737  -1,
738  1,
739  1,
740  -1,
741  -1,
742  1}; // sign of the 2nd derivative for each edge.
743  const unsigned int
744  edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2] = {
745  {0, 2},
746  {0, 2},
747  {1, 2},
748  {1, 2},
749  {0, 2},
750  {0, 2},
751  {1, 2},
752  {1, 2},
753  {0, 1},
754  {0, 1},
755  {0, 1},
756  {0, 1}}; // component which edge_lambda[m] depends on.
757  for (unsigned int m = 0; m < lines_per_cell; ++m)
758  {
759  data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
760  [edge_lambda_directions[m][1]] =
761  edge_lambda_sign[m];
762  data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
763  [edge_lambda_directions[m][0]] =
764  edge_lambda_sign[m];
765  }
766  // Precomputation for higher order shape functions,
767  // and the face parameterisation.
768  if (degree > 0)
769  {
770  // resize required data:
771  data.face_lambda_values.resize(faces_per_cell);
772  data.face_lambda_grads.resize(faces_per_cell);
773  // for face-based shape functions:
774  for (unsigned int m = 0; m < faces_per_cell; ++m)
775  {
776  data.face_lambda_values[m].resize(n_q_points);
777  data.face_lambda_grads[m].resize(3);
778  }
779  // Fill in the values (these don't change between cells).
780  for (unsigned int q = 0; q < n_q_points; ++q)
781  {
782  double x(p_list[q][0]);
783  double y(p_list[q][1]);
784  double z(p_list[q][2]);
785  data.face_lambda_values[0][q] = 1.0 - x;
786  data.face_lambda_values[1][q] = x;
787  data.face_lambda_values[2][q] = 1.0 - y;
788  data.face_lambda_values[3][q] = y;
789  data.face_lambda_values[4][q] = 1.0 - z;
790  data.face_lambda_values[5][q] = z;
791  }
792  // gradients are constant:
793  data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
794  data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
795  data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
796  data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
797  data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
798  data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
799 
800  // for cell-based shape functions:
801  // these don't depend on the cell, so can precompute all here:
802  if (flags & (update_values | update_gradients))
803  {
804  // Cell-based shape functions:
805  //
806  // Type-1 (gradients):
807  // \phi^{C_{1}}_{ijk} = grad(
808  // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
809  //
810  // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
811  const unsigned int cell_type1_offset(n_line_dofs +
812  n_face_dofs);
813  // Type-2:
814  //
815  // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
816  // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
817  // -1)\phi^{C_{1}}_{ijk}
818  //
819  // 0 <= i,j,k < degree. (subtypes in groups of
820  // degree*degree*degree)
821  //
822  // here we order so that all of subtype 1 comes first, then
823  // subtype 2.
824  const unsigned int cell_type2_offset1(
825  cell_type1_offset + degree * degree * degree);
826  const unsigned int cell_type2_offset2(
827  cell_type2_offset1 + degree * degree * degree);
828  // Type-3
829  // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
830  // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
831  // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
832  //
833  // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
834  //
835  // again we order so we compute all of subtype 1 first, then
836  // subtype 2, etc.
837  const unsigned int cell_type3_offset1(
838  cell_type2_offset2 + degree * degree * degree);
839  const unsigned int cell_type3_offset2(cell_type3_offset1 +
840  degree * degree);
841  const unsigned int cell_type3_offset3(cell_type3_offset2 +
842  degree * degree);
843 
844  // compute all points we must evaluate the 1d polynomials at:
845  std::vector<Point<dim>> cell_points(n_q_points);
846  for (unsigned int q = 0; q < n_q_points; ++q)
847  {
848  for (unsigned int d = 0; d < dim; ++d)
849  {
850  cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
851  }
852  }
853 
854 
855  // only need poly values and 1st derivative for update_values,
856  // but need 2nd derivative too for update_gradients.
857  const unsigned int poly_length(
858  (flags & update_gradients) ? 3 : 2);
859  // Loop through quad points:
860  for (unsigned int q = 0; q < n_q_points; ++q)
861  {
862  // pre-compute values & required derivatives at this quad
863  // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
864  // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
865  //
866  // for each polyc[d], c=x,y,z, contains the d-th
867  // derivative with respect to the coordinate c.
868  std::vector<std::vector<double>> polyx(
869  degree, std::vector<double>(poly_length));
870  std::vector<std::vector<double>> polyy(
871  degree, std::vector<double>(poly_length));
872  std::vector<std::vector<double>> polyz(
873  degree, std::vector<double>(poly_length));
874  for (unsigned int i = 0; i < degree; ++i)
875  {
876  // compute all required 1d polynomials for i
877  IntegratedLegendrePolynomials[i + 2].value(
878  cell_points[q][0], polyx[i]);
879  IntegratedLegendrePolynomials[i + 2].value(
880  cell_points[q][1], polyy[i]);
881  IntegratedLegendrePolynomials[i + 2].value(
882  cell_points[q][2], polyz[i]);
883  }
884  // Now use these to compute the shape functions:
885  if (flags & update_values)
886  {
887  for (unsigned int k = 0; k < degree; ++k)
888  {
889  const unsigned int shift_k(k * degree * degree);
890  const unsigned int shift_j(
891  k * degree); // Used below when subbing k for j
892  // (type 3)
893  for (unsigned int j = 0; j < degree; ++j)
894  {
895  const unsigned int shift_jk(j * degree +
896  shift_k);
897  for (unsigned int i = 0; i < degree; ++i)
898  {
899  const unsigned int shift_ijk(shift_jk +
900  i);
901 
902  // Type 1:
903  const unsigned int dof_index1(
904  cell_type1_offset + shift_ijk);
905 
906  data.shape_values[dof_index1][q][0] =
907  2.0 * polyx[i][1] * polyy[j][0] *
908  polyz[k][0];
909  data.shape_values[dof_index1][q][1] =
910  2.0 * polyx[i][0] * polyy[j][1] *
911  polyz[k][0];
912  data.shape_values[dof_index1][q][2] =
913  2.0 * polyx[i][0] * polyy[j][0] *
914  polyz[k][1];
915 
916  // Type 2:
917  const unsigned int dof_index2_1(
918  cell_type2_offset1 + shift_ijk);
919  const unsigned int dof_index2_2(
920  cell_type2_offset2 + shift_ijk);
921 
922  data.shape_values[dof_index2_1][q][0] =
923  data.shape_values[dof_index1][q][0];
924  data.shape_values[dof_index2_1][q][1] =
925  -1.0 *
926  data.shape_values[dof_index1][q][1];
927  data.shape_values[dof_index2_1][q][2] =
928  data.shape_values[dof_index1][q][2];
929 
930  data.shape_values[dof_index2_2][q][0] =
931  data.shape_values[dof_index1][q][0];
932  data.shape_values[dof_index2_2][q][1] =
933  -1.0 *
934  data.shape_values[dof_index1][q][1];
935  data.shape_values[dof_index2_2][q][2] =
936  -1.0 *
937  data.shape_values[dof_index1][q][2];
938  }
939  // Type 3: (note we re-use k and j for
940  // convenience):
941  const unsigned int shift_ij(
942  j + shift_j); // here we've subbed j for i,
943  // k for j.
944  const unsigned int dof_index3_1(
945  cell_type3_offset1 + shift_ij);
946  const unsigned int dof_index3_2(
947  cell_type3_offset2 + shift_ij);
948  const unsigned int dof_index3_3(
949  cell_type3_offset3 + shift_ij);
950 
951  data.shape_values[dof_index3_1][q][0] =
952  polyy[j][0] * polyz[k][0];
953  data.shape_values[dof_index3_1][q][1] = 0.0;
954  data.shape_values[dof_index3_1][q][2] = 0.0;
955 
956  data.shape_values[dof_index3_2][q][0] = 0.0;
957  data.shape_values[dof_index3_2][q][1] =
958  polyx[j][0] * polyz[k][0];
959  data.shape_values[dof_index3_2][q][2] = 0.0;
960 
961  data.shape_values[dof_index3_3][q][0] = 0.0;
962  data.shape_values[dof_index3_3][q][1] = 0.0;
963  data.shape_values[dof_index3_3][q][2] =
964  polyx[j][0] * polyy[k][0];
965  }
966  }
967  }
968  if (flags & update_gradients)
969  {
970  for (unsigned int k = 0; k < degree; ++k)
971  {
972  const unsigned int shift_k(k * degree * degree);
973  const unsigned int shift_j(
974  k * degree); // Used below when subbing k for j
975  // (type 3)
976  for (unsigned int j = 0; j < degree; ++j)
977  {
978  const unsigned int shift_jk(j * degree +
979  shift_k);
980  for (unsigned int i = 0; i < degree; ++i)
981  {
982  const unsigned int shift_ijk(shift_jk +
983  i);
984 
985  // Type 1:
986  const unsigned int dof_index1(
987  cell_type1_offset + shift_ijk);
988 
989  data.shape_grads[dof_index1][q][0][0] =
990  4.0 * polyx[i][2] * polyy[j][0] *
991  polyz[k][0];
992  data.shape_grads[dof_index1][q][0][1] =
993  4.0 * polyx[i][1] * polyy[j][1] *
994  polyz[k][0];
995  data.shape_grads[dof_index1][q][0][2] =
996  4.0 * polyx[i][1] * polyy[j][0] *
997  polyz[k][1];
998 
999  data.shape_grads[dof_index1][q][1][0] =
1000  data.shape_grads[dof_index1][q][0][1];
1001  data.shape_grads[dof_index1][q][1][1] =
1002  4.0 * polyx[i][0] * polyy[j][2] *
1003  polyz[k][0];
1004  data.shape_grads[dof_index1][q][1][2] =
1005  4.0 * polyx[i][0] * polyy[j][1] *
1006  polyz[k][1];
1007 
1008  data.shape_grads[dof_index1][q][2][0] =
1009  data.shape_grads[dof_index1][q][0][2];
1010  data.shape_grads[dof_index1][q][2][1] =
1011  data.shape_grads[dof_index1][q][1][2];
1012  data.shape_grads[dof_index1][q][2][2] =
1013  4.0 * polyx[i][0] * polyy[j][0] *
1014  polyz[k][2];
1015 
1016  // Type 2:
1017  const unsigned int dof_index2_1(
1018  cell_type2_offset1 + shift_ijk);
1019  const unsigned int dof_index2_2(
1020  cell_type2_offset2 + shift_ijk);
1021 
1022  for (unsigned int d = 0; d < dim; ++d)
1023  {
1024  data.shape_grads[dof_index2_1][q][0]
1025  [d] =
1026  data
1027  .shape_grads[dof_index1][q][0][d];
1028  data.shape_grads[dof_index2_1][q][1]
1029  [d] =
1030  -1.0 *
1031  data
1032  .shape_grads[dof_index1][q][1][d];
1033  data.shape_grads[dof_index2_1][q][2]
1034  [d] =
1035  data
1036  .shape_grads[dof_index1][q][2][d];
1037 
1038  data.shape_grads[dof_index2_2][q][0]
1039  [d] =
1040  data
1041  .shape_grads[dof_index1][q][0][d];
1042  data.shape_grads[dof_index2_2][q][1]
1043  [d] =
1044  -1.0 *
1045  data
1046  .shape_grads[dof_index1][q][1][d];
1047  data.shape_grads[dof_index2_2][q][2]
1048  [d] =
1049  -1.0 *
1050  data
1051  .shape_grads[dof_index1][q][2][d];
1052  }
1053  }
1054  // Type 3: (note we re-use k and j for
1055  // convenience):
1056  const unsigned int shift_ij(
1057  j + shift_j); // here we've subbed j for i,
1058  // k for j.
1059  const unsigned int dof_index3_1(
1060  cell_type3_offset1 + shift_ij);
1061  const unsigned int dof_index3_2(
1062  cell_type3_offset2 + shift_ij);
1063  const unsigned int dof_index3_3(
1064  cell_type3_offset3 + shift_ij);
1065  for (unsigned int d1 = 0; d1 < dim; ++d1)
1066  {
1067  for (unsigned int d2 = 0; d2 < dim; ++d2)
1068  {
1069  data.shape_grads[dof_index3_1][q][d1]
1070  [d2] = 0.0;
1071  data.shape_grads[dof_index3_2][q][d1]
1072  [d2] = 0.0;
1073  data.shape_grads[dof_index3_3][q][d1]
1074  [d2] = 0.0;
1075  }
1076  }
1077  data.shape_grads[dof_index3_1][q][0][1] =
1078  2.0 * polyy[j][1] * polyz[k][0];
1079  data.shape_grads[dof_index3_1][q][0][2] =
1080  2.0 * polyy[j][0] * polyz[k][1];
1081 
1082  data.shape_grads[dof_index3_2][q][1][0] =
1083  2.0 * polyx[j][1] * polyz[k][0];
1084  data.shape_grads[dof_index3_2][q][1][2] =
1085  2.0 * polyx[j][0] * polyz[k][1];
1086 
1087  data.shape_grads[dof_index3_3][q][2][0] =
1088  2.0 * polyx[j][1] * polyy[k][0];
1089  data.shape_grads[dof_index3_3][q][2][1] =
1090  2.0 * polyx[j][0] * polyy[k][1];
1091  }
1092  }
1093  }
1094  }
1095  }
1096  }
1097  break;
1098  }
1099  default:
1100  {
1101  Assert(false, ExcNotImplemented());
1102  }
1103  }
1104  return data_ptr;
1105 }
1106 
1107 
1108 
1109 template <int dim, int spacedim>
1110 void
1111 FE_NedelecSZ<dim, spacedim>::fill_edge_values(
1112  const typename Triangulation<dim, dim>::cell_iterator &cell,
1113  const Quadrature<dim> & quadrature,
1114  const InternalData & fe_data) const
1115 {
1116  // This function handles the cell-dependent construction of the EDGE-based
1117  // shape functions.
1118  //
1119  // Note it will handle both 2D and 3D, in 2D, the edges are faces, but we
1120  // handle them here.
1121  //
1122  // It will fill in the missing parts of fe_data which were not possible to
1123  // fill in the get_data routine, with respect to the edge-based shape
1124  // functions.
1125  //
1126  // It should be called by the fill_fe_*_values routines in order to complete
1127  // the basis set at quadrature points on the current cell for each edge.
1128 
1129  const UpdateFlags flags(fe_data.update_each);
1130  const unsigned int n_q_points = quadrature.size();
1131 
1132  Assert(!(flags & update_values) ||
1133  fe_data.shape_values.size() == this->n_dofs_per_cell(),
1134  ExcDimensionMismatch(fe_data.shape_values.size(),
1135  this->n_dofs_per_cell()));
1136  Assert(!(flags & update_values) ||
1137  fe_data.shape_values[0].size() == n_q_points,
1138  ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
1139 
1140  // Useful constants:
1141  const unsigned int degree(
1142  this->degree -
1143  1); // Note: constructor takes input degree + 1, so need to knock 1 off.
1144 
1145  // Useful geometry info:
1146  const unsigned int vertices_per_line(2);
1147  const unsigned int lines_per_cell(GeometryInfo<dim>::lines_per_cell);
1148 
1149  // Calculate edge orderings:
1150  std::vector<std::vector<unsigned int>> edge_order(
1151  lines_per_cell, std::vector<unsigned int>(vertices_per_line));
1152 
1153 
1154  switch (dim)
1155  {
1156  case 2:
1157  {
1158  if (flags & (update_values | update_gradients))
1159  {
1160  // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1},
1161  // e^{m}_{2}] e1 = higher global numbering of the two local
1162  // vertices e2 = lower global numbering of the two local vertices
1163  std::vector<int> edge_sign(lines_per_cell);
1164  for (unsigned int m = 0; m < lines_per_cell; ++m)
1165  {
1166  unsigned int v0_loc =
1167  GeometryInfo<dim>::line_to_cell_vertices(m, 0);
1168  unsigned int v1_loc =
1169  GeometryInfo<dim>::line_to_cell_vertices(m, 1);
1170  unsigned int v0_glob = cell->vertex_index(v0_loc);
1171  unsigned int v1_glob = cell->vertex_index(v1_loc);
1172 
1173  if (v0_glob > v1_glob)
1174  {
1175  // Opposite to global numbering on our reference element
1176  edge_sign[m] = -1.0;
1177  }
1178  else
1179  {
1180  // Aligns with global numbering on our reference element.
1181  edge_sign[m] = 1.0;
1182  }
1183  }
1184 
1185  // Define \sigma_{m} = sigma_{e^{m}_{2}} - sigma_{e^{m}_{2}}
1186  // \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
1187  //
1188  // To help things, in fe_data, we have precomputed (sigma_{i} -
1189  // sigma_{j}) and (lambda_{i} + lambda_{j}) for 0<= i,j <
1190  // lines_per_cell.
1191  //
1192  // There are two types:
1193  // - lower order (1 per edge, m):
1194  // \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
1195  //
1196  // - higher order (degree per edge, m):
1197  // \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
1198  //
1199  // NOTE: sigma_{m} and lambda_{m} are either a function of x OR
1200  // y
1201  // and if sigma is of x, then lambda is of y, and vice
1202  // versa. This means that grad(\sigma) requires
1203  // multiplication by d(sigma)/dx_{i} for the i^th comp of
1204  // grad(sigma) and similarly when taking derivatives of
1205  // lambda.
1206  //
1207  // First handle the lowest order edges (dofs 0 to 3)
1208  // 0 and 1 are the edges in the y dir. (sigma is function of y,
1209  // lambda is function of x). 2 and 3 are the edges in the x dir.
1210  // (sigma is function of x, lambda is function of y).
1211  //
1212  // More more info: see GeometryInfo for picture of the standard
1213  // element.
1214  //
1215  // Fill edge-based points:
1216  // std::vector<std::vector< Point<dim> > >
1217  // edge_points(lines_per_cell, std::vector<Point<dim>>
1218  // (n_q_points));
1219 
1220  std::vector<std::vector<double>> edge_sigma_values(
1221  fe_data.edge_sigma_values);
1222  std::vector<std::vector<double>> edge_sigma_grads(
1223  fe_data.edge_sigma_grads);
1224 
1225  std::vector<std::vector<double>> edge_lambda_values(
1226  fe_data.edge_lambda_values);
1227  std::vector<std::vector<double>> edge_lambda_grads(
1228  fe_data.edge_lambda_grads_2d);
1229 
1230  // Adjust the edge_sigma_* for the current cell:
1231  for (unsigned int m = 0; m < lines_per_cell; ++m)
1232  {
1233  std::transform(edge_sigma_values[m].begin(),
1234  edge_sigma_values[m].end(),
1235  edge_sigma_values[m].begin(),
1236  [&](const double edge_sigma_value) {
1237  return edge_sign[m] * edge_sigma_value;
1238  });
1239 
1240  std::transform(edge_sigma_grads[m].begin(),
1241  edge_sigma_grads[m].end(),
1242  edge_sigma_grads[m].begin(),
1243  [&](const double edge_sigma_grad) {
1244  return edge_sign[m] * edge_sigma_grad;
1245  });
1246  }
1247 
1248  // If we want to generate shape gradients then we need second
1249  // derivatives of the 1d polynomials, but only first derivatives
1250  // for the shape values.
1251  const unsigned int poly_length((flags & update_gradients) ? 3 :
1252  2);
1253 
1254  for (unsigned int m = 0; m < lines_per_cell; ++m)
1255  {
1256  const unsigned int shift_m(m * this->n_dofs_per_line());
1257  for (unsigned int q = 0; q < n_q_points; ++q)
1258  {
1259  // Only compute 1d polynomials if degree>0.
1260  std::vector<std::vector<double>> poly(
1261  degree, std::vector<double>(poly_length));
1262  for (unsigned int i = 1; i < degree + 1; ++i)
1263  {
1264  // Compute all required 1d polynomials and their
1265  // derivatives, starting at degree 2. e.g. to access
1266  // L'_{3}(2x-1) use polyx[1][1].
1267  IntegratedLegendrePolynomials[i + 1].value(
1268  edge_sigma_values[m][q], poly[i - 1]);
1269  }
1270  if (flags & update_values)
1271  {
1272  // Lowest order edge shape functions:
1273  for (unsigned int d = 0; d < dim; ++d)
1274  {
1275  fe_data.shape_values[shift_m][q][d] =
1276  0.5 * edge_sigma_grads[m][d] *
1277  edge_lambda_values[m][q];
1278  }
1279  // Higher order edge shape functions:
1280  for (unsigned int i = 1; i < degree + 1; ++i)
1281  {
1282  const unsigned int poly_index = i - 1;
1283  const unsigned int dof_index(i + shift_m);
1284  for (unsigned int d = 0; d < dim; ++d)
1285  {
1286  fe_data.shape_values[dof_index][q][d] =
1287  edge_sigma_grads[m][d] *
1288  poly[poly_index][1] *
1289  edge_lambda_values[m][q] +
1290  poly[poly_index][0] *
1291  edge_lambda_grads[m][d];
1292  }
1293  }
1294  }
1295  if (flags & update_gradients)
1296  {
1297  // Lowest order edge shape functions:
1298  for (unsigned int d1 = 0; d1 < dim; ++d1)
1299  {
1300  for (unsigned int d2 = 0; d2 < dim; ++d2)
1301  {
1302  // Note: gradient is constant for a given
1303  // edge.
1304  fe_data.shape_grads[shift_m][q][d1][d2] =
1305  0.5 * edge_sigma_grads[m][d1] *
1306  edge_lambda_grads[m][d2];
1307  }
1308  }
1309  // Higher order edge shape functions:
1310  for (unsigned int i = 1; i < degree + 1; ++i)
1311  {
1312  const unsigned int poly_index = i - 1;
1313  const unsigned int dof_index(i + shift_m);
1314 
1315  fe_data.shape_grads[dof_index][q][0][0] =
1316  edge_sigma_grads[m][0] *
1317  edge_sigma_grads[m][0] *
1318  edge_lambda_values[m][q] * poly[poly_index][2];
1319 
1320  fe_data.shape_grads[dof_index][q][0][1] =
1321  (edge_sigma_grads[m][0] *
1322  edge_lambda_grads[m][1] +
1323  edge_sigma_grads[m][1] *
1324  edge_lambda_grads[m][0]) *
1325  poly[poly_index][1];
1326 
1327  fe_data.shape_grads[dof_index][q][1][0] =
1328  fe_data.shape_grads[dof_index][q][0][1];
1329 
1330  fe_data.shape_grads[dof_index][q][1][1] =
1331  edge_sigma_grads[m][1] *
1332  edge_sigma_grads[m][1] *
1333  edge_lambda_values[m][q] * poly[poly_index][2];
1334  }
1335  }
1336  }
1337  }
1338  }
1339  break;
1340  }
1341  case 3:
1342  {
1343  if (flags & (update_values | update_gradients))
1344  {
1345  // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1},
1346  // e^{m}_{2}] e1 = higher global numbering of the two local
1347  // vertices e2 = lower global numbering of the two local vertices
1348  std::vector<int> edge_sign(lines_per_cell);
1349  for (unsigned int m = 0; m < lines_per_cell; ++m)
1350  {
1351  unsigned int v0_loc =
1352  GeometryInfo<dim>::line_to_cell_vertices(m, 0);
1353  unsigned int v1_loc =
1354  GeometryInfo<dim>::line_to_cell_vertices(m, 1);
1355  unsigned int v0_glob = cell->vertex_index(v0_loc);
1356  unsigned int v1_glob = cell->vertex_index(v1_loc);
1357 
1358  if (v0_glob > v1_glob)
1359  {
1360  // Opposite to global numbering on our reference element
1361  edge_sign[m] = -1.0;
1362  }
1363  else
1364  {
1365  // Aligns with global numbering on our reference element.
1366  edge_sign[m] = 1.0;
1367  }
1368  }
1369 
1370  // Define \sigma_{m} = sigma_{e^{m}_{2}} - sigma_{e^{m}_{2}}
1371  // \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
1372  //
1373  // To help things, in fe_data, we have precomputed (sigma_{i} -
1374  // sigma_{j}) and (lambda_{i} + lambda_{j}) for 0<= i,j <
1375  // lines_per_cell.
1376  //
1377  // There are two types:
1378  // - lower order (1 per edge, m):
1379  // \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
1380  //
1381  // - higher order (degree per edge, m):
1382  // \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
1383  //
1384  // NOTE: In the ref cell, sigma_{m} is a function of x OR y OR Z
1385  // and lambda_{m} a function of the remaining co-ords.
1386  // for example, if sigma is of x, then lambda is of y AND
1387  // z, and so on. This means that grad(\sigma) requires
1388  // multiplication by d(sigma)/dx_{i} for the i^th comp of
1389  // grad(sigma) and similarly when taking derivatives of
1390  // lambda.
1391  //
1392  // First handle the lowest order edges (dofs 0 to 11)
1393  // 0 and 1 are the edges in the y dir at z=0. (sigma is a fn of y,
1394  // lambda is a fn of x & z). 2 and 3 are the edges in the x dir at
1395  // z=0. (sigma is a fn of x, lambda is a fn of y & z). 4 and 5 are
1396  // the edges in the y dir at z=1. (sigma is a fn of y, lambda is a
1397  // fn of x & z). 6 and 7 are the edges in the x dir at z=1. (sigma
1398  // is a fn of x, lambda is a fn of y & z). 8 and 9 are the edges
1399  // in the z dir at y=0. (sigma is a fn of z, lambda is a fn of x &
1400  // y). 10 and 11 are the edges in the z dir at y=1. (sigma is a fn
1401  // of z, lambda is a fn of x & y).
1402  //
1403  // For more info: see GeometryInfo for picture of the standard
1404  // element.
1405 
1406  // Copy over required edge-based data:
1407  std::vector<std::vector<double>> edge_sigma_values(
1408  fe_data.edge_sigma_values);
1409  std::vector<std::vector<double>> edge_lambda_values(
1410  fe_data.edge_lambda_values);
1411  std::vector<std::vector<double>> edge_sigma_grads(
1412  fe_data.edge_sigma_grads);
1413  std::vector<std::vector<std::vector<double>>> edge_lambda_grads(
1414  fe_data.edge_lambda_grads_3d);
1415  std::vector<std::vector<std::vector<double>>>
1416  edge_lambda_gradgrads_3d(fe_data.edge_lambda_gradgrads_3d);
1417 
1418  // Adjust the edge_sigma_* for the current cell:
1419  for (unsigned int m = 0; m < lines_per_cell; ++m)
1420  {
1421  std::transform(edge_sigma_values[m].begin(),
1422  edge_sigma_values[m].end(),
1423  edge_sigma_values[m].begin(),
1424  [&](const double edge_sigma_value) {
1425  return edge_sign[m] * edge_sigma_value;
1426  });
1427  std::transform(edge_sigma_grads[m].begin(),
1428  edge_sigma_grads[m].end(),
1429  edge_sigma_grads[m].begin(),
1430  [&](const double edge_sigma_grad) {
1431  return edge_sign[m] * edge_sigma_grad;
1432  });
1433  }
1434 
1435  // Now calculate the edge-based shape functions:
1436  // If we want to generate shape gradients then we need second
1437  // derivatives of the 1d polynomials, but only first derivatives
1438  // for the shape values.
1439  const unsigned int poly_length((flags & update_gradients) ? 3 :
1440  2);
1441  std::vector<std::vector<double>> poly(
1442  degree, std::vector<double>(poly_length));
1443  for (unsigned int m = 0; m < lines_per_cell; ++m)
1444  {
1445  const unsigned int shift_m(m * this->n_dofs_per_line());
1446  for (unsigned int q = 0; q < n_q_points; ++q)
1447  {
1448  // precompute values of all 1d polynomials required:
1449  if (degree > 0)
1450  {
1451  for (unsigned int i = 0; i < degree; ++i)
1452  {
1453  IntegratedLegendrePolynomials[i + 2].value(
1454  edge_sigma_values[m][q], poly[i]);
1455  }
1456  }
1457  if (flags & update_values)
1458  {
1459  // Lowest order shape functions:
1460  for (unsigned int d = 0; d < dim; ++d)
1461  {
1462  fe_data.shape_values[shift_m][q][d] =
1463  0.5 * edge_sigma_grads[m][d] *
1464  edge_lambda_values[m][q];
1465  }
1466  if (degree > 0)
1467  {
1468  for (unsigned int i = 0; i < degree; ++i)
1469  {
1470  const unsigned int dof_index(i + 1 + shift_m);
1471  for (unsigned int d = 0; d < dim; ++d)
1472  {
1473  fe_data.shape_values[dof_index][q][d] =
1474  edge_sigma_grads[m][d] * poly[i][1] *
1475  edge_lambda_values[m][q] +
1476  poly[i][0] * edge_lambda_grads[m][q][d];
1477  }
1478  }
1479  }
1480  }
1481  if (flags & update_gradients)
1482  {
1483  // Lowest order shape functions:
1484  for (unsigned int d1 = 0; d1 < dim; ++d1)
1485  {
1486  for (unsigned int d2 = 0; d2 < dim; ++d2)
1487  {
1488  fe_data.shape_grads[shift_m][q][d1][d2] =
1489  0.5 * edge_sigma_grads[m][d1] *
1490  edge_lambda_grads[m][q][d2];
1491  }
1492  }
1493  if (degree > 0)
1494  {
1495  for (unsigned int i = 0; i < degree; ++i)
1496  {
1497  const unsigned int dof_index(i + 1 + shift_m);
1498 
1499  for (unsigned int d1 = 0; d1 < dim; ++d1)
1500  {
1501  for (unsigned int d2 = 0; d2 < dim; ++d2)
1502  {
1503  fe_data
1504  .shape_grads[dof_index][q][d1][d2] =
1505  edge_sigma_grads[m][d1] *
1506  edge_sigma_grads[m][d2] *
1507  poly[i][2] *
1508  edge_lambda_values[m][q] +
1509  (edge_sigma_grads[m][d1] *
1510  edge_lambda_grads[m][q][d2] +
1511  edge_sigma_grads[m][d2] *
1512  edge_lambda_grads[m][q][d1]) *
1513  poly[i][1] +
1514  edge_lambda_gradgrads_3d[m][d1]
1515  [d2] *
1516  poly[i][0];
1517  }
1518  }
1519  }
1520  }
1521  }
1522  }
1523  }
1524  }
1525  break;
1526  }
1527  default:
1528  {
1529  Assert(false, ExcNotImplemented());
1530  }
1531  }
1532 }
1533 
1534 
1535 
1536 template <int dim, int spacedim>
1537 void
1538 FE_NedelecSZ<dim, spacedim>::fill_face_values(
1539  const typename Triangulation<dim, dim>::cell_iterator &cell,
1540  const Quadrature<dim> & quadrature,
1541  const InternalData & fe_data) const
1542 {
1543  // This function handles the cell-dependent construction of the FACE-based
1544  // shape functions.
1545  //
1546  // Note that it should only be called in 3D.
1547  Assert(dim == 3, ExcDimensionMismatch(dim, 3));
1548  //
1549  // It will fill in the missing parts of fe_data which were not possible to
1550  // fill in the get_data routine, with respect to face-based shape functions.
1551  //
1552  // It should be called by the fill_fe_*_values routines in order to complete
1553  // the basis set at quadrature points on the current cell for each face.
1554 
1555  // Useful constants:
1556  const unsigned int degree(
1557  this->degree -
1558  1); // Note: constructor takes input degree + 1, so need to knock 1 off.
1559 
1560  // Do nothing if FE degree is 0.
1561  if (degree > 0)
1562  {
1563  const UpdateFlags flags(fe_data.update_each);
1564 
1565  if (flags & (update_values | update_gradients))
1566  {
1567  const unsigned int n_q_points = quadrature.size();
1568 
1569  Assert(!(flags & update_values) ||
1570  fe_data.shape_values.size() == this->n_dofs_per_cell(),
1571  ExcDimensionMismatch(fe_data.shape_values.size(),
1572  this->n_dofs_per_cell()));
1573  Assert(!(flags & update_values) ||
1574  fe_data.shape_values[0].size() == n_q_points,
1575  ExcDimensionMismatch(fe_data.shape_values[0].size(),
1576  n_q_points));
1577 
1578  // Useful geometry info:
1579  const unsigned int vertices_per_face(
1580  GeometryInfo<dim>::vertices_per_face);
1581  const unsigned int faces_per_cell(GeometryInfo<dim>::faces_per_cell);
1582 
1583  // DoF info:
1584  const unsigned int n_line_dofs =
1585  this->n_dofs_per_line() * GeometryInfo<dim>::lines_per_cell;
1586 
1587  // First we find the global face orientations on the current cell.
1588  std::vector<std::vector<unsigned int>> face_orientation(
1589  faces_per_cell, std::vector<unsigned int>(vertices_per_face));
1590 
1591  const unsigned int
1592  vertex_opposite_on_face[GeometryInfo<3>::vertices_per_face] = {3,
1593  2,
1594  1,
1595  0};
1596 
1597  const unsigned int
1598  vertices_adjacent_on_face[GeometryInfo<3>::vertices_per_face][2] = {
1599  {1, 2}, {0, 3}, {0, 3}, {1, 2}};
1600 
1601  for (unsigned int m = 0; m < faces_per_cell; ++m)
1602  {
1603  // Find the local vertex on this face with the highest global
1604  // numbering. This is f^m_0.
1605  unsigned int current_max = 0;
1606  unsigned int current_glob = cell->vertex_index(
1607  GeometryInfo<dim>::face_to_cell_vertices(m, 0));
1608  for (unsigned int v = 1; v < vertices_per_face; ++v)
1609  {
1610  if (current_glob <
1611  cell->vertex_index(
1612  GeometryInfo<dim>::face_to_cell_vertices(m, v)))
1613  {
1614  current_max = v;
1615  current_glob = cell->vertex_index(
1616  GeometryInfo<dim>::face_to_cell_vertices(m, v));
1617  }
1618  }
1619  face_orientation[m][0] =
1620  GeometryInfo<dim>::face_to_cell_vertices(m, current_max);
1621 
1622  // f^m_2 is the vertex opposite f^m_0.
1623  face_orientation[m][2] = GeometryInfo<dim>::face_to_cell_vertices(
1624  m, vertex_opposite_on_face[current_max]);
1625 
1626  // Finally, f^m_1 is the vertex with the greater global numbering
1627  // of the remaining two local vertices. Then, f^m_3 is the other.
1628  if (cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(
1629  m, vertices_adjacent_on_face[current_max][0])) >
1630  cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(
1631  m, vertices_adjacent_on_face[current_max][1])))
1632  {
1633  face_orientation[m][1] =
1634  GeometryInfo<dim>::face_to_cell_vertices(
1635  m, vertices_adjacent_on_face[current_max][0]);
1636  face_orientation[m][3] =
1637  GeometryInfo<dim>::face_to_cell_vertices(
1638  m, vertices_adjacent_on_face[current_max][1]);
1639  }
1640  else
1641  {
1642  face_orientation[m][1] =
1643  GeometryInfo<dim>::face_to_cell_vertices(
1644  m, vertices_adjacent_on_face[current_max][1]);
1645  face_orientation[m][3] =
1646  GeometryInfo<dim>::face_to_cell_vertices(
1647  m, vertices_adjacent_on_face[current_max][0]);
1648  }
1649  }
1650 
1651  // Now we know the face orientation on the current cell, we can
1652  // generate the parameterisation:
1653  std::vector<std::vector<double>> face_xi_values(
1654  faces_per_cell, std::vector<double>(n_q_points));
1655  std::vector<std::vector<double>> face_xi_grads(
1656  faces_per_cell, std::vector<double>(dim));
1657  std::vector<std::vector<double>> face_eta_values(
1658  faces_per_cell, std::vector<double>(n_q_points));
1659  std::vector<std::vector<double>> face_eta_grads(
1660  faces_per_cell, std::vector<double>(dim));
1661 
1662  std::vector<std::vector<double>> face_lambda_values(
1663  faces_per_cell, std::vector<double>(n_q_points));
1664  std::vector<std::vector<double>> face_lambda_grads(
1665  faces_per_cell, std::vector<double>(dim));
1666  for (unsigned int m = 0; m < faces_per_cell; ++m)
1667  {
1668  for (unsigned int q = 0; q < n_q_points; ++q)
1669  {
1670  face_xi_values[m][q] =
1671  fe_data.sigma_imj_values[q][face_orientation[m][0]]
1672  [face_orientation[m][1]];
1673  face_eta_values[m][q] =
1674  fe_data.sigma_imj_values[q][face_orientation[m][0]]
1675  [face_orientation[m][3]];
1676  face_lambda_values[m][q] = fe_data.face_lambda_values[m][q];
1677  }
1678  for (unsigned int d = 0; d < dim; ++d)
1679  {
1680  face_xi_grads[m][d] =
1681  fe_data.sigma_imj_grads[face_orientation[m][0]]
1682  [face_orientation[m][1]][d];
1683  face_eta_grads[m][d] =
1684  fe_data.sigma_imj_grads[face_orientation[m][0]]
1685  [face_orientation[m][3]][d];
1686 
1687  face_lambda_grads[m][d] = fe_data.face_lambda_grads[m][d];
1688  }
1689  }
1690  // Now can generate the basis
1691  const unsigned int poly_length((flags & update_gradients) ? 3 : 2);
1692  std::vector<std::vector<double>> polyxi(
1693  degree, std::vector<double>(poly_length));
1694  std::vector<std::vector<double>> polyeta(
1695  degree, std::vector<double>(poly_length));
1696 
1697  // Loop through quad points:
1698  for (unsigned int m = 0; m < faces_per_cell; ++m)
1699  {
1700  // we assume that all quads have the same number of dofs
1701  const unsigned int shift_m(m * this->n_dofs_per_quad(0));
1702  // Calculate the offsets for each face-based shape function:
1703  //
1704  // Type-1 (gradients)
1705  // \phi^{F_m,1}_{ij} = \nabla( L_{i+2}(\xi_{F_{m}})
1706  // L_{j+2}(\eta_{F_{m}}) \lambda_{F_{m}} )
1707  //
1708  // 0 <= i,j < degree (in a group of degree*degree)
1709  const unsigned int face_type1_offset(n_line_dofs + shift_m);
1710  // Type-2:
1711  //
1712  // \phi^{F_m,2}_{ij} = ( L'_{i+2}(\xi_{F_{m}})
1713  // L_{j+2}(\eta_{F_{m}}) \nabla\xi_{F_{m}}
1714  // - L_{i+2}(\xi_{F_{m}})
1715  // L'_{j+2}(\eta_{F_{m}}) \nabla\eta_{F_{m}}
1716  // ) \lambda_{F_{m}}
1717  //
1718  // 0 <= i,j < degree (in a group of degree*degree)
1719  const unsigned int face_type2_offset(face_type1_offset +
1720  degree * degree);
1721  // Type-3:
1722  //
1723  // \phi^{F_m,3}_{i} = L_{i+2}(\eta_{F_{m}}) \lambda_{F_{m}}
1724  // \nabla\xi_{F_{m}} \phi^{F_m,3}_{i+p} = L_{i+2}(\xi_{F_{m}})
1725  // \lambda_{F_{m}} \nabla\eta_{F_{m}}
1726  //
1727  // 0 <= i < degree.
1728  //
1729  // here we order so that all of subtype 1 comes first, then
1730  // subtype 2.
1731  const unsigned int face_type3_offset1(face_type2_offset +
1732  degree * degree);
1733  const unsigned int face_type3_offset2(face_type3_offset1 +
1734  degree);
1735 
1736  // Loop over all faces:
1737  for (unsigned int q = 0; q < n_q_points; ++q)
1738  {
1739  // pre-compute values & required derivatives at this quad
1740  // point: polyxi = L_{i+2}(\xi_{F_{m}}), polyeta =
1741  // L_{j+2}(\eta_{F_{m}}),
1742  //
1743  // each polypoint[k][d], contains the dth derivative of
1744  // L_{k+2} at the point \xi or \eta. Note that this doesn't
1745  // include the derivative of xi/eta via the chain rule.
1746  for (unsigned int i = 0; i < degree; ++i)
1747  {
1748  // compute all required 1d polynomials:
1749  IntegratedLegendrePolynomials[i + 2].value(
1750  face_xi_values[m][q], polyxi[i]);
1751  IntegratedLegendrePolynomials[i + 2].value(
1752  face_eta_values[m][q], polyeta[i]);
1753  }
1754  // Now use these to compute the shape functions:
1755  if (flags & update_values)
1756  {
1757  for (unsigned int j = 0; j < degree; ++j)
1758  {
1759  const unsigned int shift_j(j * degree);
1760  for (unsigned int i = 0; i < degree; ++i)
1761  {
1762  const unsigned int shift_ij(shift_j + i);
1763  // Type 1:
1764  const unsigned int dof_index1(face_type1_offset +
1765  shift_ij);
1766  for (unsigned int d = 0; d < dim; ++d)
1767  {
1768  fe_data.shape_values[dof_index1][q][d] =
1769  (face_xi_grads[m][d] * polyxi[i][1] *
1770  polyeta[j][0] +
1771  face_eta_grads[m][d] * polyxi[i][0] *
1772  polyeta[j][1]) *
1773  face_lambda_values[m][q] +
1774  face_lambda_grads[m][d] * polyxi[i][0] *
1775  polyeta[j][0];
1776  }
1777  // Type 2:
1778  const unsigned int dof_index2(face_type2_offset +
1779  shift_ij);
1780  for (unsigned int d = 0; d < dim; ++d)
1781  {
1782  fe_data.shape_values[dof_index2][q][d] =
1783  (face_xi_grads[m][d] * polyxi[i][1] *
1784  polyeta[j][0] -
1785  face_eta_grads[m][d] * polyxi[i][0] *
1786  polyeta[j][1]) *
1787  face_lambda_values[m][q];
1788  }
1789  }
1790  // Type 3:
1791  const unsigned int dof_index3_1(face_type3_offset1 +
1792  j);
1793  const unsigned int dof_index3_2(face_type3_offset2 +
1794  j);
1795  for (unsigned int d = 0; d < dim; ++d)
1796  {
1797  fe_data.shape_values[dof_index3_1][q][d] =
1798  face_xi_grads[m][d] * polyeta[j][0] *
1799  face_lambda_values[m][q];
1800 
1801  fe_data.shape_values[dof_index3_2][q][d] =
1802  face_eta_grads[m][d] * polyxi[j][0] *
1803  face_lambda_values[m][q];
1804  }
1805  }
1806  }
1807  if (flags & update_gradients)
1808  {
1809  for (unsigned int j = 0; j < degree; ++j)
1810  {
1811  const unsigned int shift_j(j * degree);
1812  for (unsigned int i = 0; i < degree; ++i)
1813  {
1814  const unsigned int shift_ij(shift_j + i);
1815  // Type 1:
1816  const unsigned int dof_index1(face_type1_offset +
1817  shift_ij);
1818  for (unsigned int d1 = 0; d1 < dim; ++d1)
1819  {
1820  for (unsigned int d2 = 0; d2 < dim; ++d2)
1821  {
1822  fe_data
1823  .shape_grads[dof_index1][q][d1][d2] =
1824  (face_xi_grads[m][d1] *
1825  face_xi_grads[m][d2] * polyxi[i][2] *
1826  polyeta[j][0] +
1827  (face_xi_grads[m][d1] *
1828  face_eta_grads[m][d2] +
1829  face_xi_grads[m][d2] *
1830  face_eta_grads[m][d1]) *
1831  polyxi[i][1] * polyeta[j][1] +
1832  face_eta_grads[m][d1] *
1833  face_eta_grads[m][d2] *
1834  polyxi[i][0] * polyeta[j][2]) *
1835  face_lambda_values[m][q] +
1836  (face_xi_grads[m][d2] * polyxi[i][1] *
1837  polyeta[j][0] +
1838  face_eta_grads[m][d2] * polyxi[i][0] *
1839  polyeta[j][1]) *
1840  face_lambda_grads[m][d1] +
1841  (face_xi_grads[m][d1] * polyxi[i][1] *
1842  polyeta[j][0] +
1843  face_eta_grads[m][d1] * polyxi[i][0] *
1844  polyeta[j][1]) *
1845  face_lambda_grads[m][d2];
1846  }
1847  }
1848  // Type 2:
1849  const unsigned int dof_index2(face_type2_offset +
1850  shift_ij);
1851  for (unsigned int d1 = 0; d1 < dim; ++d1)
1852  {
1853  for (unsigned int d2 = 0; d2 < dim; ++d2)
1854  {
1855  fe_data
1856  .shape_grads[dof_index2][q][d1][d2] =
1857  (face_xi_grads[m][d1] *
1858  face_xi_grads[m][d2] * polyxi[i][2] *
1859  polyeta[j][0] +
1860  (face_xi_grads[m][d1] *
1861  face_eta_grads[m][d2] -
1862  face_xi_grads[m][d2] *
1863  face_eta_grads[m][d1]) *
1864  polyxi[i][1] * polyeta[j][1] -
1865  face_eta_grads[m][d1] *
1866  face_eta_grads[m][d2] *
1867  polyxi[i][0] * polyeta[j][2]) *
1868  face_lambda_values[m][q] +
1869  (face_xi_grads[m][d1] * polyxi[i][1] *
1870  polyeta[j][0] -
1871  face_eta_grads[m][d1] * polyxi[i][0] *
1872  polyeta[j][1]) *
1873  face_lambda_grads[m][d2];
1874  }
1875  }
1876  }
1877  // Type 3:
1878  const unsigned int dof_index3_1(face_type3_offset1 +
1879  j);
1880  const unsigned int dof_index3_2(face_type3_offset2 +
1881  j);
1882  for (unsigned int d1 = 0; d1 < dim; ++d1)
1883  {
1884  for (unsigned int d2 = 0; d2 < dim; ++d2)
1885  {
1886  fe_data.shape_grads[dof_index3_1][q][d1][d2] =
1887  face_xi_grads[m][d1] *
1888  (face_eta_grads[m][d2] * polyeta[j][1] *
1889  face_lambda_values[m][q] +
1890  face_lambda_grads[m][d2] * polyeta[j][0]);
1891 
1892  fe_data.shape_grads[dof_index3_2][q][d1][d2] =
1893  face_eta_grads[m][d1] *
1894  (face_xi_grads[m][d2] * polyxi[j][1] *
1895  face_lambda_values[m][q] +
1896  face_lambda_grads[m][d2] * polyxi[j][0]);
1897  }
1898  }
1899  }
1900  }
1901  }
1902  }
1903  }
1904  if (flags & update_hessians)
1905  {
1906  Assert(false, ExcNotImplemented());
1907  }
1908  }
1909 }
1910 
1911 
1912 
1913 template <int dim, int spacedim>
1914 void
1915 FE_NedelecSZ<dim, spacedim>::fill_fe_values(
1916  const typename Triangulation<dim, dim>::cell_iterator &cell,
1917  const CellSimilarity::Similarity /*cell_similarity*/,
1918  const Quadrature<dim> & quadrature,
1919  const Mapping<dim, dim> & mapping,
1920  const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
1921  const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
1922  & mapping_data,
1923  const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
1924  ::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim>
1925  &data) const
1926 {
1927  // Convert to the correct internal data class for this FE class.
1928  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
1929  ExcInternalError());
1930  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
1931 
1932  // Now update the edge-based DoFs, which depend on the cell.
1933  // This will fill in the missing items in the InternalData
1934  // (fe_internal/fe_data) which was not filled in by get_data.
1935  fill_edge_values(cell, quadrature, fe_data);
1936  if (dim == 3 && this->degree > 1)
1937  {
1938  fill_face_values(cell, quadrature, fe_data);
1939  }
1940 
1941  const UpdateFlags flags(fe_data.update_each);
1942  const unsigned int n_q_points = quadrature.size();
1943 
1944  Assert(!(flags & update_values) ||
1945  fe_data.shape_values.size() == this->n_dofs_per_cell(),
1946  ExcDimensionMismatch(fe_data.shape_values.size(),
1947  this->n_dofs_per_cell()));
1948  Assert(!(flags & update_values) ||
1949  fe_data.shape_values[0].size() == n_q_points,
1950  ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
1951 
1952  if (flags & update_values)
1953  {
1954  // Now have all shape_values stored on the reference cell.
1955  // Must now transform to the physical cell.
1956  std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
1957  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1958  {
1959  const unsigned int first =
1960  data.shape_function_to_row_table[dof * this->n_components() +
1961  this->get_nonzero_components(dof)
1962  .first_selected_component()];
1963 
1964  mapping.transform(make_array_view(fe_data.shape_values[dof]),
1965  mapping_covariant,
1966  mapping_internal,
1967  make_array_view(transformed_shape_values));
1968  for (unsigned int q = 0; q < n_q_points; ++q)
1969  {
1970  for (unsigned int d = 0; d < dim; ++d)
1971  {
1972  data.shape_values(first + d, q) =
1973  transformed_shape_values[q][d];
1974  }
1975  }
1976  }
1977  }
1978  if (flags & update_gradients)
1979  {
1980  // Now have all shape_grads stored on the reference cell.
1981  // Must now transform to the physical cell.
1982  std::vector<Tensor<2, dim>> input(n_q_points);
1983  std::vector<Tensor<2, dim>> transformed_shape_grads(n_q_points);
1984  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1985  {
1986  for (unsigned int q = 0; q < n_q_points; ++q)
1987  {
1988  input[q] = fe_data.shape_grads[dof][q];
1989  }
1990  mapping.transform(make_array_view(input),
1991  mapping_covariant_gradient,
1992  mapping_internal,
1993  make_array_view(transformed_shape_grads));
1994 
1995  const unsigned int first =
1996  data.shape_function_to_row_table[dof * this->n_components() +
1997  this->get_nonzero_components(dof)
1998  .first_selected_component()];
1999 
2000  for (unsigned int q = 0; q < n_q_points; ++q)
2001  {
2002  for (unsigned int d1 = 0; d1 < dim; ++d1)
2003  {
2004  for (unsigned int d2 = 0; d2 < dim; ++d2)
2005  {
2006  transformed_shape_grads[q][d1] -=
2007  data.shape_values(first + d2, q) *
2008  mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
2009  }
2010  }
2011  }
2012 
2013  for (unsigned int q = 0; q < n_q_points; ++q)
2014  {
2015  for (unsigned int d = 0; d < dim; ++d)
2016  {
2017  data.shape_gradients[first + d][q] =
2018  transformed_shape_grads[q][d];
2019  }
2020  }
2021  }
2022  }
2023 }
2024 
2025 
2026 
2027 template <int dim, int spacedim>
2028 void
2029 FE_NedelecSZ<dim, spacedim>::fill_fe_face_values(
2030  const typename Triangulation<dim, dim>::cell_iterator &cell,
2031  const unsigned int face_no,
2032  const hp::QCollection<dim - 1> & quadrature,
2033  const Mapping<dim, dim> & mapping,
2034  const typename Mapping<dim, dim>::InternalDataBase & mapping_internal,
2035  const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2036  & mapping_data,
2037  const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
2038  ::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim>
2039  &data) const
2040 {
2041  AssertDimension(quadrature.size(), 1);
2042 
2043  // Note for future improvement:
2044  // We don't have the full quadrature - should use QProjector to create the 2D
2045  // quadrature.
2046  //
2047  // For now I am effectively generating all of the shape function vals/grads,
2048  // etc. On all quad points on all faces and then only using them for one face.
2049  // This is obviously inefficient. I should cache the cell number and cache
2050  // all of the shape_values/gradients etc and then reuse them for each face.
2051 
2052  // convert data object to internal
2053  // data for this class. fails with
2054  // an exception if that is not
2055  // possible
2056  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
2057  ExcInternalError());
2058  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
2059 
2060  // Now update the edge-based DoFs, which depend on the cell.
2061  // This will fill in the missing items in the InternalData
2062  // (fe_internal/fe_data) which was not filled in by get_data.
2063  fill_edge_values(cell,
2064  QProjector<dim>::project_to_all_faces(this->reference_cell(),
2065  quadrature[0]),
2066  fe_data);
2067  if (dim == 3 && this->degree > 1)
2068  {
2069  fill_face_values(cell,
2070  QProjector<dim>::project_to_all_faces(
2071  this->reference_cell(), quadrature[0]),
2072  fe_data);
2073  }
2074 
2075  const UpdateFlags flags(fe_data.update_each);
2076  const unsigned int n_q_points = quadrature[0].size();
2077  const auto offset =
2078  QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
2079  face_no,
2080  cell->face_orientation(face_no),
2081  cell->face_flip(face_no),
2082  cell->face_rotation(face_no),
2083  n_q_points);
2084 
2085  if (flags & update_values)
2086  {
2087  // Now have all shape_values stored on the reference cell.
2088  // Must now transform to the physical cell.
2089  std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
2090  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2091  {
2092  mapping.transform(make_array_view(fe_data.shape_values[dof],
2093  offset,
2094  n_q_points),
2095  mapping_covariant,
2096  mapping_internal,
2097  make_array_view(transformed_shape_values));
2098 
2099  const unsigned int first =
2100  data.shape_function_to_row_table[dof * this->n_components() +
2101  this->get_nonzero_components(dof)
2102  .first_selected_component()];
2103 
2104  for (unsigned int q = 0; q < n_q_points; ++q)
2105  {
2106  for (unsigned int d = 0; d < dim; ++d)
2107  {
2108  data.shape_values(first + d, q) =
2109  transformed_shape_values[q][d];
2110  }
2111  }
2112  }
2113  }
2114  if (flags & update_gradients)
2115  {
2116  // Now have all shape_grads stored on the reference cell.
2117  // Must now transform to the physical cell.
2118  std::vector<Tensor<2, dim>> input(n_q_points);
2119  std::vector<Tensor<2, dim>> transformed_shape_grads(n_q_points);
2120  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2121  {
2122  for (unsigned int q = 0; q < n_q_points; ++q)
2123  {
2124  input[q] = fe_data.shape_grads[dof][offset + q];
2125  }
2126  mapping.transform(input,
2127  mapping_covariant_gradient,
2128  mapping_internal,
2129  make_array_view(transformed_shape_grads));
2130 
2131  const unsigned int first =
2132  data.shape_function_to_row_table[dof * this->n_components() +
2133  this->get_nonzero_components(dof)
2134  .first_selected_component()];
2135 
2136  for (unsigned int q = 0; q < n_q_points; ++q)
2137  {
2138  for (unsigned int d1 = 0; d1 < dim; ++d1)
2139  {
2140  for (unsigned int d2 = 0; d2 < dim; ++d2)
2141  {
2142  transformed_shape_grads[q][d1] -=
2143  data.shape_values(first + d2, q) *
2144  mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
2145  }
2146  }
2147  }
2148 
2149  for (unsigned int q = 0; q < n_q_points; ++q)
2150  {
2151  for (unsigned int d = 0; d < dim; ++d)
2152  {
2153  data.shape_gradients[first + d][q] =
2154  transformed_shape_grads[q][d];
2155  }
2156  }
2157  }
2158  }
2159 }
2160 
2161 
2162 
2163 template <int dim, int spacedim>
2164 void
2165 FE_NedelecSZ<dim, spacedim>::fill_fe_subface_values(
2166  const typename Triangulation<dim, dim>::cell_iterator & /*cell*/,
2167  const unsigned int /*face_no*/,
2168  const unsigned int /*sub_no*/,
2169  const Quadrature<dim - 1> & /*quadrature*/,
2170  const Mapping<dim, dim> & /*mapping*/,
2171  const typename Mapping<dim, dim>::InternalDataBase & /*mapping_internal*/,
2172  const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2173  & /*mapping_data*/,
2174  const typename FiniteElement<dim, dim>::InternalDataBase & /*fe_internal*/,
2175  ::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim>
2176  & /*data*/) const
2177 {
2178  Assert(false, ExcNotImplemented());
2179 }
2180 
2181 
2182 
2183 template <int dim, int spacedim>
2184 UpdateFlags
2185 FE_NedelecSZ<dim, spacedim>::requires_update_flags(
2186  const UpdateFlags flags) const
2187 {
2188  UpdateFlags out = update_default;
2189 
2190  if (flags & update_values)
2191  out |= update_values | update_covariant_transformation;
2192 
2193  if (flags & update_gradients)
2194  out |= update_gradients | update_values |
2195  update_jacobian_pushed_forward_grads |
2196  update_covariant_transformation;
2197 
2198  if (flags & update_hessians)
2199  // Assert (false, ExcNotImplemented());
2200  out |= update_hessians | update_values | update_gradients |
2201  update_jacobian_pushed_forward_grads |
2202  update_jacobian_pushed_forward_2nd_derivatives |
2203  update_covariant_transformation;
2204 
2205  return out;
2206 }
2207 
2208 
2209 
2210 template <int dim, int spacedim>
2211 std::string
2212 FE_NedelecSZ<dim, spacedim>::get_name() const
2213 {
2214  // note that the FETools::get_fe_by_name function depends on the particular
2215  // format of the string this function returns, so they have to be kept in sync
2216  std::ostringstream namebuf;
2217  namebuf << "FE_NedelecSZ<" << dim << ">(" << this->degree - 1 << ")";
2218 
2219  return namebuf.str();
2220 }
2221 
2222 
2223 
2224 template <int dim, int spacedim>
2225 std::unique_ptr<FiniteElement<dim, dim>>
2226 FE_NedelecSZ<dim, spacedim>::clone() const
2227 {
2228  return std::make_unique<FE_NedelecSZ<dim, spacedim>>(*this);
2229 }
2230 
2231 
2232 
2233 template <int dim, int spacedim>
2234 std::vector<unsigned int>
2235 FE_NedelecSZ<dim, spacedim>::get_dpo_vector(const unsigned int degree)
2236 {
2237  // internal function to return a vector of "dofs per object"
2238  // where the objects inside the vector refer to:
2239  // 0 = vertex
2240  // 1 = edge
2241  // 2 = face (which is a cell in 2D)
2242  // 3 = cell
2243  std::vector<unsigned int> dpo;
2244 
2245  dpo.push_back(0);
2246  dpo.push_back(degree + 1);
2247  if (dim > 1)
2248  dpo.push_back(2 * degree * (degree + 1));
2249  if (dim > 2)
2250  dpo.push_back(3 * degree * degree * (degree + 1));
2251 
2252  return dpo;
2253 }
2254 
2255 
2256 
2257 template <int dim, int spacedim>
2258 unsigned int
2259 FE_NedelecSZ<dim, spacedim>::compute_num_dofs(const unsigned int degree) const
2260 {
2261  // Internal function to compute the number of DoFs
2262  // for a given dimension & polynomial order.
2263  switch (dim)
2264  {
2265  case 2:
2266  return 2 * (degree + 1) * (degree + 2);
2267 
2268  case 3:
2269  return 3 * (degree + 1) * (degree + 2) * (degree + 2);
2270 
2271  default:
2272  {
2273  Assert(false, ExcNotImplemented());
2274  return 0;
2275  }
2276  }
2277 }
2278 
2279 
2280 
2281 template <int dim, int spacedim>
2282 void
2283 FE_NedelecSZ<dim, spacedim>::create_polynomials(const unsigned int degree)
2284 {
2285  // fill the 1d polynomials vector:
2286  IntegratedLegendrePolynomials =
2287  IntegratedLegendreSZ::generate_complete_basis(degree + 1);
2288 }
2289 
2290 
2291 
2292 // explicit instantiations
2293 #include "fe_nedelec_sz.inst"
2294 
2295 DEAL_II_NAMESPACE_CLOSE
ArrayView::make_array_view
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:697
FE_NedelecSZ::InternalData
Definition: fe_nedelec_sz.h:239
FE_NedelecSZ::InternalData::edge_lambda_grads_2d
std::vector< std::vector< double > > edge_lambda_grads_2d
Definition: fe_nedelec_sz.h:344
FE_NedelecSZ::InternalData::shape_values
std::vector< std::vector< Tensor< 1, dim > > > shape_values
Definition: fe_nedelec_sz.h:247
FE_NedelecSZ::InternalData::sigma_imj_values
std::vector< std::vector< std::vector< double > > > sigma_imj_values
Definition: fe_nedelec_sz.h:273
FE_NedelecSZ::InternalData::edge_sigma_grads
std::vector< std::vector< double > > edge_sigma_grads
Definition: fe_nedelec_sz.h:316
FE_NedelecSZ::InternalData::edge_sigma_values
std::vector< std::vector< double > > edge_sigma_values
Definition: fe_nedelec_sz.h:303
FE_NedelecSZ::InternalData::sigma_imj_grads
std::vector< std::vector< std::vector< double > > > sigma_imj_grads
Definition: fe_nedelec_sz.h:290
FE_NedelecSZ::InternalData::face_lambda_values
std::vector< std::vector< double > > face_lambda_values
Definition: fe_nedelec_sz.h:385
FE_NedelecSZ::InternalData::shape_grads
std::vector< std::vector< DerivativeForm< 1, dim, dim > > > shape_grads
Definition: fe_nedelec_sz.h:256
FE_NedelecSZ::InternalData::edge_lambda_grads_3d
std::vector< std::vector< std::vector< double > > > edge_lambda_grads_3d
Definition: fe_nedelec_sz.h:355
FE_NedelecSZ::InternalData::edge_lambda_values
std::vector< std::vector< double > > edge_lambda_values
Definition: fe_nedelec_sz.h:333
FE_NedelecSZ::InternalData::edge_lambda_gradgrads_3d
std::vector< std::vector< std::vector< double > > > edge_lambda_gradgrads_3d
Definition: fe_nedelec_sz.h:367
FE_NedelecSZ::InternalData::face_lambda_grads
std::vector< std::vector< double > > face_lambda_grads
Definition: fe_nedelec_sz.h:395
FE_NedelecSZ::fill_fe_values
virtual void fill_fe_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const Mapping< dim, dim > &mapping, const typename Mapping< dim, dim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, dim > &mapping_data, const typename FiniteElement< dim, dim >::InternalDataBase &fedata, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, dim > &data) const override
Definition: fe_nedelec_sz.cc:1915
FE_NedelecSZ::compute_num_dofs
unsigned int compute_num_dofs(const unsigned int degree) const
Definition: fe_nedelec_sz.cc:2259
FE_NedelecSZ::clone
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
Definition: fe_nedelec_sz.cc:2226
FE_NedelecSZ::fill_fe_face_values
virtual void fill_fe_face_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const Mapping< dim, dim > &mapping, const typename Mapping< dim, dim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, dim > &mapping_data, const typename FiniteElement< dim, dim >::InternalDataBase &fedata, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, dim > &data) const override
Definition: fe_nedelec_sz.cc:2029
FE_NedelecSZ::FE_NedelecSZ
FE_NedelecSZ(const unsigned int order)
Definition: fe_nedelec_sz.cc:25
FE_NedelecSZ::requires_update_flags
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: fe_nedelec_sz.cc:2185
FE_NedelecSZ::get_dpo_vector
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_nedelec_sz.cc:2235
FE_NedelecSZ::fill_fe_subface_values
virtual void fill_fe_subface_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim - 1 > &quadrature, const Mapping< dim, dim > &mapping, const typename Mapping< dim, dim >::InternalDataBase &mapping_internal, const ::internal::FEValuesImplementation::MappingRelatedData< dim, dim > &mapping_data, const typename FiniteElement< dim, dim >::InternalDataBase &fedata, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, dim > &data) const override
Definition: fe_nedelec_sz.cc:2165
FE_NedelecSZ::create_polynomials
void create_polynomials(const unsigned int degree)
Definition: fe_nedelec_sz.cc:2283
FE_NedelecSZ::get_data
virtual std::unique_ptr< typename ::FiniteElement< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags update_flags, const Mapping< dim, spacedim > &mapping, const Quadrature< dim > &quadrature, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
Definition: fe_nedelec_sz.cc:128
FE_NedelecSZ::get_name
virtual std::string get_name() const override
Definition: fe_nedelec_sz.cc:2212
FE_NedelecSZ::shape_grad_grad_component
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_nedelec_sz.cc:115
FE_NedelecSZ::shape_grad_component
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_nedelec_sz.cc:91
FE_NedelecSZ::shape_grad
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
Definition: fe_nedelec_sz.cc:80
FE_NedelecSZ::shape_grad_grad
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
Definition: fe_nedelec_sz.cc:104
FE_NedelecSZ::shape_value_component
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
Definition: fe_nedelec_sz.cc:66
FE_NedelecSZ::fill_face_values
void fill_face_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const Quadrature< dim > &quadrature, const InternalData &fedata) const
Definition: fe_nedelec_sz.cc:1538
FE_NedelecSZ::fill_edge_values
void fill_edge_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const Quadrature< dim > &quadrature, const InternalData &fedata) const
Definition: fe_nedelec_sz.cc:1111
FE_NedelecSZ::mapping_kind
MappingKind mapping_kind
Definition: fe_nedelec_sz.h:152
FE_NedelecSZ::shape_value
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
Definition: fe_nedelec_sz.cc:55
FiniteElementData::n_components
unsigned int n_components() const
FiniteElement< dim, dim >::component_to_base_table
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition: fe.h:2565
IntegratedLegendreSZ::generate_complete_basis
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomials_integrated_legendre_sz.cc:76
Mapping::transform
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const =0
Point
Definition: point.h:111
QProjector::DataSetDescriptor::face
static DataSetDescriptor face(const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Definition: qprojector.cc:1365
Quadrature::get_points
const std::vector< Point< dim > > & get_points() const
Quadrature::size
unsigned int size() const
hp::Collection::size
unsigned int size() const
Definition: collection.h:109
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:396
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:397
fe_nedelec_sz.h
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:67
update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_2nd_derivatives
Definition: fe_update_flags.h:198
update_jacobian_pushed_forward_grads
@ update_jacobian_pushed_forward_grads
Definition: fe_update_flags.h:189
update_hessians
@ update_hessians
Second derivatives of shape functions.
Definition: fe_update_flags.h:90
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
update_covariant_transformation
@ update_covariant_transformation
Covariant transformation.
Definition: fe_update_flags.h:160
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
update_default
@ update_default
No update.
Definition: fe_update_flags.h:69
first
Point< 2 > first
Definition: grid_out.cc:4587
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
mapping_covariant_gradient
@ mapping_covariant_gradient
Definition: mapping.h:86
mapping_covariant
@ mapping_covariant
Definition: mapping.h:75
mapping_nedelec
@ mapping_nedelec
Definition: mapping.h:115
FE_P_BubblesImplementation::get_dpo_vector
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_simplex_p_bubbles.cc:118
GridGenerator::reference_cell
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:53
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:67
parallel::transform
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:213
GeometryInfo::face_to_cell_vertices
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
GeometryInfo::line_to_cell_vertices
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
Generated by doxygen 1.9.1