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