Reference documentation for deal.II version GIT defb42778c 2022-12-05 01:15:02+00:00
fe_nedelec_sz.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
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
18
19 #include <memory>
20
22
23 // Constructor:
24 template <int dim, int spacedim>
26  : FiniteElement<dim, dim>(
28  dim,
29  order + 1,
30  FiniteElementData<dim>::Hcurl),
31  std::vector<bool>(compute_num_dofs(order), true),
33  std::vector<bool>(dim, true)))
34 {
35  Assert(dim >= 2, ExcImpossibleInDim(dim));
36
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
56  const Point<dim> & /*p*/) const
57 {
59  return 0.;
60 }
61
62
63
64 template <int dim, int spacedim>
65 double
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>
80 FE_NedelecSZ<dim, spacedim>::shape_grad(const unsigned int /*i*/,
81  const Point<dim> & /*p*/) const
82 {
84  return Tensor<1, dim>();
85 }
86
87
88
89 template <int dim, int spacedim>
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>
105  const Point<dim> & /*p*/) const
106 {
108  return Tensor<2, dim>();
109 }
110
111
112
113 template <int dim, int spacedim>
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>
129  const UpdateFlags update_flags,
130  const Mapping<dim, spacedim> & /*mapping*/,
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
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
175  {
177  std::vector<DerivativeForm<1, dim, dim>>(
178  n_q_points));
179  }
180
181  if (flags & update_hessians)
182  {
183  data.shape_hessians.resize(this->n_dofs_per_cell(),
184  std::vector<DerivativeForm<2, dim, dim>>(
185  n_q_points));
186  }
187
188  std::vector<Point<dim>> p_list(n_q_points);
190
191
192  switch (dim)
193  {
194  case 2:
195  {
196  // Compute values of sigma & lambda and the sigma differences and
198  std::vector<std::vector<double>> sigma(
199  n_q_points, std::vector<double>(lines_per_cell));
200  std::vector<std::vector<double>> lambda(
201  n_q_points, std::vector<double>(lines_per_cell));
202
203  for (unsigned int q = 0; q < n_q_points; ++q)
204  {
205  sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
206  sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
207  sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
208  sigma[q][3] = p_list[q][0] + p_list[q][1];
209
210  lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
211  lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
212  lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
213  lambda[q][3] = p_list[q][0] * p_list[q][1];
214  for (unsigned int i = 0; i < vertices_per_cell; ++i)
215  {
216  for (unsigned int j = 0; j < vertices_per_cell; ++j)
217  {
218  data.sigma_imj_values[q][i][j] =
219  sigma[q][i] - sigma[q][j];
220  }
221  }
222  }
223
224  // Calculate the gradient of sigma_imj_values[q][i][j] =
225  // sigma[q][i]-sigma[q][j]
226  // - this depends on the component and the direction of the
227  // corresponding edge.
228  // - the direction of the edge is determined by
229  // sigma_imj_sign[i][j].
230  // Helper arrays:
231  const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2] = {
232  {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
233  int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
234  unsigned int sigma_imj_component[vertices_per_cell]
235  [vertices_per_cell];
236
237  for (unsigned int i = 0; i < vertices_per_cell; ++i)
238  {
239  for (unsigned int j = 0; j < vertices_per_cell; ++j)
240  {
241  // sigma_imj_sign is the sign (+/-) of the coefficient of
242  // x/y/z in sigma_imj_values Due to the numbering of vertices
243  // on the reference element it is easy to find edges in the
244  // positive direction are from smaller to higher local vertex
245  // numbering.
246  sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
247  sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
248
249  // Now store the component which the sigma_i - sigma_j
250  // corresponds to:
251  sigma_imj_component[i][j] = 0;
252  for (unsigned int d = 0; d < dim; ++d)
253  {
254  int temp_imj =
255  sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
256  // Only interested in the first non-zero
257  // as if there is a second, it can not be a valid edge.
258  if (temp_imj != 0)
259  {
260  sigma_imj_component[i][j] = d;
261  break;
262  }
263  }
264  // Can now calculate the gradient, only non-zero in the
265  // component given: Note some i,j combinations will be
266  // incorrect, but only on invalid edges.
268  2.0 * sigma_imj_sign[i][j];
269  }
270  }
271
272  // Now compute the edge parameterisations for a single element
273  // with global numbering matching that of the reference element:
274
275  // Resize the edge parameterisations
276  data.edge_sigma_values.resize(lines_per_cell);
278  for (unsigned int m = 0; m < lines_per_cell; ++m)
279  {
280  data.edge_sigma_values[m].resize(n_q_points);
281
282  // sigma grads are constant in a cell (no need for quad points)
284  }
285
286  // Fill the values for edge lambda and edge sigma:
287  const unsigned int
288  edge_sigma_direction[GeometryInfo<2>::lines_per_cell] = {1,
289  1,
290  0,
291  0};
292
293  data.edge_lambda_values.resize(lines_per_cell,
294  std::vector<double>(n_q_points));
296  std::vector<double>(dim));
297  for (unsigned int m = 0; m < lines_per_cell; ++m)
298  {
299  // e1=max(reference vertex numbering on this edge)
300  // e2=min(reference vertex numbering on this edge)
301  // Which is guaranteed to be:
302  const unsigned int e1(
304  const unsigned int e2(
306  for (unsigned int q = 0; q < n_q_points; ++q)
307  {
308  data.edge_sigma_values[m][q] =
309  data.sigma_imj_values[q][e2][e1];
310  data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
311  }
312
314  }
319
320  // If the polynomial order is 0, then no more work to do:
321  if (degree < 1)
322  {
323  break;
324  }
325
326  // Otherwise, we can compute the non-cell dependent shape functions.
327  //
328  // Note: the local dof numberings follow the usual order of lines ->
329  // faces -> cells
330  // (we have no vertex-based DoFs in this element).
331  // For a given cell we have:
332  // n_line_dofs = dofs_per_line*lines_per_cell.
333  // n_face_dofs = dofs_per_face*faces_per_cell.
334  // n_cell_dofs = dofs_per_quad (2D)
335  // = dofs_per_hex (3D)
336  //
337  // i.e. For the local dof numbering:
338  // the first line dof is 0,
339  // the first face dof is n_line_dofs,
340  // the first cell dof is n_line_dofs + n_face_dofs.
341  //
342  // On a line, DoFs are ordered first by line_dof and then line_index:
343  // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
344  //
345  // and similarly for faces:
346  // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
347  //
348  // HOWEVER, we have different types of DoFs on a line/face/cell.
349  // On a line we have two types, lowest order and higher order
351  // - The numbering is such the lowest order is first, then higher
352  // order.
353  // This is simple enough as there is only 1 lowest order and
354  // degree higher orders DoFs per line.
355  //
356  // On a 2D cell, we have 3 types: Type 1/2/3:
357  // - The ordering done by type:
358  // - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
359  // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index
360  // = ij1.
361  // - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
362  // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index
363  // = degree^2 + ij2
364  // - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
365  // Numbered: ij3 = i3. i.e. cell_dof_index
366  // = 2*(degree^2) + ij3.
367  //
368  // These then fit into the local dof numbering described above:
369  // - local dof numberings are:
370  // line_dofs: local_dof = line_dof_index. 0 <= local_dof <
371  // dofs_per_line*lines_per_cell face_dofs: local_dof =
372  // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
373  // = n_lines_dof + n_face_dofs + cell_dof_index.
374  //
375  // The cell-based shape functions are:
376  //
378  // \phi^{C_{1}}_{ij) = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
379  //
380  // 0 <= i,j < degree.
381  //
382  // NOTE: The derivative produced by IntegratedLegendrePolynomials does
383  // not account for the
384  // (2*x-1) or (2*y-1) so we must take this into account when
385  // taking derivatives.
386  const unsigned int cell_type1_offset = n_line_dofs;
387
388  // Type 2:
389  // \phi^{C_{2}}_{ij) = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
390  // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
391  //
392  // 0 <= i,j < degree.
393  const unsigned int cell_type2_offset =
394  cell_type1_offset + degree * degree;
395
396  // Type 3 (two subtypes):
397  // \phi^{C_{3}}_{j) = L_{j+2}(2y-1) \mathbf{e}_{x}
398  //
399  // \phi^{C_{3}}_{j+degree) = L_{j+2}(2x-1) \mathbf{e}_{y}
400  //
401  // 0 <= j < degree
402  const unsigned int cell_type3_offset1 =
403  cell_type2_offset + degree * degree;
404  const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
405
407  {
408  // compute all points we must evaluate the 1d polynomials at:
409  std::vector<Point<dim>> cell_points(n_q_points);
410  for (unsigned int q = 0; q < n_q_points; ++q)
411  {
412  for (unsigned int d = 0; d < dim; ++d)
413  {
414  cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
415  }
416  }
417
418  // Loop through quad points:
419  for (unsigned int q = 0; q < n_q_points; ++q)
420  {
421  // pre-compute values & required derivatives at this quad
422  // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
423  //
424  // for each polyc[d], c=x,y, contains the d-th derivative with
425  // respect to the coordinate c.
426
427  // We only need poly values and 1st derivative for
428  // update_values, but need the 2nd derivative too for
429  // update_gradients. For update_hessians we also need the 3rd
430  // derivatives.
431  const unsigned int poly_length =
432  (flags & update_hessians) ?
433  4 :
434  ((flags & update_gradients) ? 3 : 2);
435
436  std::vector<std::vector<double>> polyx(
437  degree, std::vector<double>(poly_length));
438  std::vector<std::vector<double>> polyy(
439  degree, std::vector<double>(poly_length));
440  for (unsigned int i = 0; i < degree; ++i)
441  {
442  // Compute all required 1d polynomials and their
443  // derivatives, starting at degree 2. e.g. to access
444  // L'_{3}(2x-1) use polyx[1][1].
445  IntegratedLegendrePolynomials[i + 2].value(
446  cell_points[q][0], polyx[i]);
447  IntegratedLegendrePolynomials[i + 2].value(
448  cell_points[q][1], polyy[i]);
449  }
450  // Now use these to compute the shape functions:
451  if (flags & update_values)
452  {
453  for (unsigned int j = 0; j < degree; ++j)
454  {
455  const unsigned int shift_j(j * degree);
456  for (unsigned int i = 0; i < degree; ++i)
457  {
458  const unsigned int shift_ij(i + shift_j);
459
460  // Type 1:
461  const unsigned int dof_index1(cell_type1_offset +
462  shift_ij);
463  data.shape_values[dof_index1][q][0] =
464  2.0 * polyx[i][1] * polyy[j][0];
465  data.shape_values[dof_index1][q][1] =
466  2.0 * polyx[i][0] * polyy[j][1];
467
468  // Type 2:
469  const unsigned int dof_index2(cell_type2_offset +
470  shift_ij);
471  data.shape_values[dof_index2][q][0] =
472  data.shape_values[dof_index1][q][0];
473  data.shape_values[dof_index2][q][1] =
474  -1.0 * data.shape_values[dof_index1][q][1];
475  }
476  // Type 3:
477  const unsigned int dof_index3_1(cell_type3_offset1 +
478  j);
479  data.shape_values[dof_index3_1][q][0] = polyy[j][0];
480  data.shape_values[dof_index3_1][q][1] = 0.0;
481
482  const unsigned int dof_index3_2(cell_type3_offset2 +
483  j);
484  data.shape_values[dof_index3_2][q][0] = 0.0;
485  data.shape_values[dof_index3_2][q][1] = polyx[j][0];
486  }
487  }
489  {
490  for (unsigned int j = 0; j < degree; ++j)
491  {
492  const unsigned int shift_j(j * degree);
493  for (unsigned int i = 0; i < degree; ++i)
494  {
495  const unsigned int shift_ij(i + shift_j);
496
497  // Type 1:
498  const unsigned int dof_index1(cell_type1_offset +
499  shift_ij);
501  4.0 * polyx[i][2] * polyy[j][0];
503  4.0 * polyx[i][1] * polyy[j][1];
507  4.0 * polyx[i][0] * polyy[j][2];
508
509  // Type 2:
510  const unsigned int dof_index2(cell_type2_offset +
511  shift_ij);
520  }
521  // Type 3:
522  const unsigned int dof_index3_1(cell_type3_offset1 +
523  j);
526  2.0 * polyy[j][1];
529
530  const unsigned int dof_index3_2(cell_type3_offset2 +
531  j);
535  2.0 * polyx[j][1];
537  }
538  }
539  if (flags & update_hessians)
540  {
541  for (unsigned int j = 0; j < degree; ++j)
542  {
543  const unsigned int shift_j(j * degree);
544  for (unsigned int i = 0; i < degree; ++i)
545  {
546  const unsigned int shift_ij(i + shift_j);
547
548  // Type 1:
549  const unsigned int dof_index1(cell_type1_offset +
550  shift_ij);
551  data.shape_hessians[dof_index1][q][0][0][0] =
552  8.0 * polyx[i][3] * polyy[j][0];
553  data.shape_hessians[dof_index1][q][1][0][0] =
554  8.0 * polyx[i][2] * polyy[j][1];
555
556  data.shape_hessians[dof_index1][q][0][1][0] =
557  data.shape_hessians[dof_index1][q][1][0][0];
558  data.shape_hessians[dof_index1][q][1][1][0] =
559  8.0 * polyx[i][1] * polyy[j][2];
560
561  data.shape_hessians[dof_index1][q][0][0][1] =
562  data.shape_hessians[dof_index1][q][1][0][0];
563  data.shape_hessians[dof_index1][q][1][0][1] =
564  data.shape_hessians[dof_index1][q][1][1][0];
565
566  data.shape_hessians[dof_index1][q][0][1][1] =
567  data.shape_hessians[dof_index1][q][1][1][0];
568  data.shape_hessians[dof_index1][q][1][1][1] =
569  8.0 * polyx[i][0] * polyy[j][3];
570
571
572
573  // Type 2:
574  const unsigned int dof_index2(cell_type2_offset +
575  shift_ij);
576  for (unsigned int d = 0; d < dim; ++d)
577  {
578  data.shape_hessians[dof_index2][q][0][0][d] =
579  data.shape_hessians[dof_index1][q][0][0][d];
580  data.shape_hessians[dof_index2][q][0][1][d] =
581  data.shape_hessians[dof_index1][q][0][1][d];
582  data.shape_hessians[dof_index2][q][1][0][d] =
583  -1.0 *
584  data.shape_hessians[dof_index1][q][1][0][d];
585  data.shape_hessians[dof_index2][q][1][1][d] =
586  -1.0 *
587  data.shape_hessians[dof_index1][q][1][1][d];
588  }
589  }
590  // Type 3:
591  const unsigned int dof_index3_1(cell_type3_offset1 +
592  j);
593  data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
594  data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
595  data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
596  data.shape_hessians[dof_index3_1][q][0][1][1] =
597  4.0 * polyy[j][2];
598  data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
599  data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
600  data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
601  data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
602
603  const unsigned int dof_index3_2(cell_type3_offset2 +
604  j);
605  data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
606  data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
607  data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
608  data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
609  data.shape_hessians[dof_index3_2][q][1][0][0] =
610  4.0 * polyx[j][2];
611  data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
612  data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
613  data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
614  }
615  }
616  }
617  }
618  break;
619  }
620  case 3:
621  {
622  // Compute values of sigma & lambda and the sigma differences and
624  std::vector<std::vector<double>> sigma(
625  n_q_points, std::vector<double>(lines_per_cell));
626  std::vector<std::vector<double>> lambda(
627  n_q_points, std::vector<double>(lines_per_cell));
628  for (unsigned int q = 0; q < n_q_points; ++q)
629  {
630  sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
631  (1 - p_list[q][2]);
632  sigma[q][1] =
633  p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
634  sigma[q][2] =
635  (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
636  sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
637  sigma[q][4] =
638  (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
639  sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
640  sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
641  sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
642
643  lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
644  (1.0 - p_list[q][2]);
645  lambda[q][1] =
646  p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
647  lambda[q][2] =
648  (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
649  lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
650  lambda[q][4] =
651  (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
652  lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
653  lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
654  lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
655
656  // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
657  // and lambda_ipj = \lambda_{i} + \lambda_{j}.
658  for (unsigned int i = 0; i < vertices_per_cell; ++i)
659  {
660  for (unsigned int j = 0; j < vertices_per_cell; ++j)
661  {
662  data.sigma_imj_values[q][i][j] =
663  sigma[q][i] - sigma[q][j];
664  }
665  }
666  }
667
669  // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
670  // calculate values & derivatives of the shape functions we need to
671  // know:
672  // - The component the sigma_imj value corresponds to - this varies
673  // with i & j.
674  // - The gradient of the sigma_imj value
675  // - this depends on the component and the direction of the
676  // corresponding edge.
677  // - the direction of the edge is determined by
678  // sigma_imj_sign[i][j].
679  //
680  // Note that not every i,j combination is a valid edge (there are only
681  // 12 valid edges in 3D), but we compute them all as it simplifies
682  // things.
683
684  // store the sign of each component x, y, z in the sigma list.
685  // can use this to fill in the sigma_imj_component data.
686  const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3] = {
687  {-1, -1, -1},
688  {1, -1, -1},
689  {-1, 1, -1},
690  {1, 1, -1},
691  {-1, -1, 1},
692  {1, -1, 1},
693  {-1, 1, 1},
694  {1, 1, 1}};
695
696  int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
697  unsigned int sigma_imj_component[vertices_per_cell]
698  [vertices_per_cell];
699
700  for (unsigned int i = 0; i < vertices_per_cell; ++i)
701  {
702  for (unsigned int j = 0; j < vertices_per_cell; ++j)
703  {
704  // sigma_imj_sign is the sign (+/-) of the coefficient of
705  // x/y/z in sigma_imj. Due to the numbering of vertices on the
706  // reference element this is easy to work out because edges in
707  // the positive direction go from smaller to higher local
708  // vertex numbering.
709  sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
710  sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
711
712  // Now store the component which the sigma_i - sigma_j
713  // corresponds to:
714  sigma_imj_component[i][j] = 0;
715  for (unsigned int d = 0; d < dim; ++d)
716  {
717  int temp_imj =
718  sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
719  // Only interested in the first non-zero
720  // as if there is a second, it will not be a valid edge.
721  if (temp_imj != 0)
722  {
723  sigma_imj_component[i][j] = d;
724  break;
725  }
726  }
727  // Can now calculate the gradient, only non-zero in the
728  // component given: Note some i,j combinations will be
729  // incorrect, but only on invalid edges.
731  2.0 * sigma_imj_sign[i][j];
732  }
733  }
734  // Now compute the edge parameterisations for a single element
735  // with global numbering matching that of the reference element:
736
737  // resize the edge parameterisations
738  data.edge_sigma_values.resize(lines_per_cell);
739  data.edge_lambda_values.resize(lines_per_cell);
743  for (unsigned int m = 0; m < lines_per_cell; ++m)
744  {
745  data.edge_sigma_values[m].resize(n_q_points);
746  data.edge_lambda_values[m].resize(n_q_points);
747
748  // sigma grads are constant in a cell (no need for quad points)
750
752  for (unsigned int q = 0; q < n_q_points; ++q)
753  {
755  }
757  // points)
759  for (unsigned int d = 0; d < dim; ++d)
760  {
762  }
763  }
764
765  // Fill the values:
766  const unsigned int
767  edge_sigma_direction[GeometryInfo<3>::lines_per_cell] = {
768  1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
769
770  for (unsigned int m = 0; m < lines_per_cell; ++m)
771  {
772  // e1=max(reference vertex numbering on this edge)
773  // e2=min(reference vertex numbering on this edge)
774  // Which is guaranteed to be:
775  const unsigned int e1(
777  const unsigned int e2(
779
780  for (unsigned int q = 0; q < n_q_points; ++q)
781  {
782  data.edge_sigma_values[m][q] =
783  data.sigma_imj_values[q][e2][e1];
784  data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
785  }
786
788  }
790  for (unsigned int q = 0; q < n_q_points; ++q)
791  {
792  double x(p_list[q][0]);
793  double y(p_list[q][1]);
794  double z(p_list[q][2]);
795  data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
796  data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
797  data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
798  data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
799  data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
800  data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
801  data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
802  data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
803  data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
804  data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
805  data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
806  data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
807  }
809  const int edge_lambda_sign[GeometryInfo<3>::lines_per_cell] = {
810  1,
811  -1,
812  1,
813  -1,
814  -1,
815  1,
816  -1,
817  1,
818  1,
819  -1,
820  -1,
821  1}; // sign of the 2nd derivative for each edge.
822  const unsigned int
823  edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2] = {
824  {0, 2},
825  {0, 2},
826  {1, 2},
827  {1, 2},
828  {0, 2},
829  {0, 2},
830  {1, 2},
831  {1, 2},
832  {0, 1},
833  {0, 1},
834  {0, 1},
835  {0, 1}}; // component which edge_lambda[m] depends on.
836  for (unsigned int m = 0; m < lines_per_cell; ++m)
837  {
839  [edge_lambda_directions[m][1]] =
840  edge_lambda_sign[m];
842  [edge_lambda_directions[m][0]] =
843  edge_lambda_sign[m];
844  }
845  // Precomputation for higher order shape functions,
846  // and the face parameterisation.
847  if (degree > 0)
848  {
849  // resize required data:
850  data.face_lambda_values.resize(faces_per_cell);
852  // for face-based shape functions:
853  for (unsigned int m = 0; m < faces_per_cell; ++m)
854  {
855  data.face_lambda_values[m].resize(n_q_points);
857  }
858  // Fill in the values (these don't change between cells).
859  for (unsigned int q = 0; q < n_q_points; ++q)
860  {
861  double x(p_list[q][0]);
862  double y(p_list[q][1]);
863  double z(p_list[q][2]);
864  data.face_lambda_values[0][q] = 1.0 - x;
865  data.face_lambda_values[1][q] = x;
866  data.face_lambda_values[2][q] = 1.0 - y;
867  data.face_lambda_values[3][q] = y;
868  data.face_lambda_values[4][q] = 1.0 - z;
869  data.face_lambda_values[5][q] = z;
870  }
872  data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
873  data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
874  data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
875  data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
876  data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
877  data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
878
879  // for cell-based shape functions:
880  // these don't depend on the cell, so can precompute all here:
882  {
883  // Cell-based shape functions:
884  //
887  // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
888  //
889  // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
890  const unsigned int cell_type1_offset(n_line_dofs +
891  n_face_dofs);
892  // Type-2:
893  //
894  // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
895  // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
896  // -1)\phi^{C_{1}}_{ijk}
897  //
898  // 0 <= i,j,k < degree. (subtypes in groups of
899  // degree*degree*degree)
900  //
901  // here we order so that all of subtype 1 comes first, then
902  // subtype 2.
903  const unsigned int cell_type2_offset1(
904  cell_type1_offset + degree * degree * degree);
905  const unsigned int cell_type2_offset2(
906  cell_type2_offset1 + degree * degree * degree);
907  // Type-3
908  // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
909  // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
910  // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
911  //
912  // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
913  //
914  // again we order so we compute all of subtype 1 first, then
915  // subtype 2, etc.
916  const unsigned int cell_type3_offset1(
917  cell_type2_offset2 + degree * degree * degree);
918  const unsigned int cell_type3_offset2(cell_type3_offset1 +
919  degree * degree);
920  const unsigned int cell_type3_offset3(cell_type3_offset2 +
921  degree * degree);
922
923  // compute all points we must evaluate the 1d polynomials at:
924  std::vector<Point<dim>> cell_points(n_q_points);
925  for (unsigned int q = 0; q < n_q_points; ++q)
926  {
927  for (unsigned int d = 0; d < dim; ++d)
928  {
929  cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
930  }
931  }
932
933  // We only need poly values and 1st derivative for
934  // update_values, but need the 2nd derivative too for
935  // update_gradients. For update_hessians we also need 3rd
936  // derivative.
937  const unsigned int poly_length =
938  (flags & update_hessians) ?
939  4 :
940  ((flags & update_gradients) ? 3 : 2);
941
942  // Loop through quad points:
943  for (unsigned int q = 0; q < n_q_points; ++q)
944  {
945  // pre-compute values & required derivatives at this quad
946  // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
947  // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
948  //
949  // for each polyc[d], c=x,y,z, contains the d-th
950  // derivative with respect to the coordinate c.
951  std::vector<std::vector<double>> polyx(
952  degree, std::vector<double>(poly_length));
953  std::vector<std::vector<double>> polyy(
954  degree, std::vector<double>(poly_length));
955  std::vector<std::vector<double>> polyz(
956  degree, std::vector<double>(poly_length));
957  for (unsigned int i = 0; i < degree; ++i)
958  {
959  // compute all required 1d polynomials for i
960  IntegratedLegendrePolynomials[i + 2].value(
961  cell_points[q][0], polyx[i]);
962  IntegratedLegendrePolynomials[i + 2].value(
963  cell_points[q][1], polyy[i]);
964  IntegratedLegendrePolynomials[i + 2].value(
965  cell_points[q][2], polyz[i]);
966  }
967  // Now use these to compute the shape functions:
968  if (flags & update_values)
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_values[dof_index1][q][0] =
990  2.0 * polyx[i][1] * polyy[j][0] *
991  polyz[k][0];
992  data.shape_values[dof_index1][q][1] =
993  2.0 * polyx[i][0] * polyy[j][1] *
994  polyz[k][0];
995  data.shape_values[dof_index1][q][2] =
996  2.0 * polyx[i][0] * polyy[j][0] *
997  polyz[k][1];
998
999  // Type 2:
1000  const unsigned int dof_index2_1(
1001  cell_type2_offset1 + shift_ijk);
1002  const unsigned int dof_index2_2(
1003  cell_type2_offset2 + shift_ijk);
1004
1005  data.shape_values[dof_index2_1][q][0] =
1006  data.shape_values[dof_index1][q][0];
1007  data.shape_values[dof_index2_1][q][1] =
1008  -1.0 *
1009  data.shape_values[dof_index1][q][1];
1010  data.shape_values[dof_index2_1][q][2] =
1011  data.shape_values[dof_index1][q][2];
1012
1013  data.shape_values[dof_index2_2][q][0] =
1014  data.shape_values[dof_index1][q][0];
1015  data.shape_values[dof_index2_2][q][1] =
1016  -1.0 *
1017  data.shape_values[dof_index1][q][1];
1018  data.shape_values[dof_index2_2][q][2] =
1019  -1.0 *
1020  data.shape_values[dof_index1][q][2];
1021  }
1022  // Type 3: (note we re-use k and j for
1023  // convenience):
1024  const unsigned int shift_ij(
1025  j + shift_j); // here we've subbed j for i,
1026  // k for j.
1027  const unsigned int dof_index3_1(
1028  cell_type3_offset1 + shift_ij);
1029  const unsigned int dof_index3_2(
1030  cell_type3_offset2 + shift_ij);
1031  const unsigned int dof_index3_3(
1032  cell_type3_offset3 + shift_ij);
1033
1034  data.shape_values[dof_index3_1][q][0] =
1035  polyy[j][0] * polyz[k][0];
1036  data.shape_values[dof_index3_1][q][1] = 0.0;
1037  data.shape_values[dof_index3_1][q][2] = 0.0;
1038
1039  data.shape_values[dof_index3_2][q][0] = 0.0;
1040  data.shape_values[dof_index3_2][q][1] =
1041  polyx[j][0] * polyz[k][0];
1042  data.shape_values[dof_index3_2][q][2] = 0.0;
1043
1044  data.shape_values[dof_index3_3][q][0] = 0.0;
1045  data.shape_values[dof_index3_3][q][1] = 0.0;
1046  data.shape_values[dof_index3_3][q][2] =
1047  polyx[j][0] * polyy[k][0];
1048  }
1049  }
1050  }
1052  {
1053  for (unsigned int k = 0; k < degree; ++k)
1054  {
1055  const unsigned int shift_k(k * degree * degree);
1056  const unsigned int shift_j(
1057  k * degree); // Used below when subbing k for j
1058  // (type 3)
1059  for (unsigned int j = 0; j < degree; ++j)
1060  {
1061  const unsigned int shift_jk(j * degree +
1062  shift_k);
1063  for (unsigned int i = 0; i < degree; ++i)
1064  {
1065  const unsigned int shift_ijk(shift_jk +
1066  i);
1067
1068  // Type 1:
1069  const unsigned int dof_index1(
1070  cell_type1_offset + shift_ijk);
1071
1073  4.0 * polyx[i][2] * polyy[j][0] *
1074  polyz[k][0];
1076  4.0 * polyx[i][1] * polyy[j][1] *
1077  polyz[k][0];
1079  4.0 * polyx[i][1] * polyy[j][0] *
1080  polyz[k][1];
1081
1085  4.0 * polyx[i][0] * polyy[j][2] *
1086  polyz[k][0];
1088  4.0 * polyx[i][0] * polyy[j][1] *
1089  polyz[k][1];
1090
1096  4.0 * polyx[i][0] * polyy[j][0] *
1097  polyz[k][2];
1098
1099  // Type 2:
1100  const unsigned int dof_index2_1(
1101  cell_type2_offset1 + shift_ijk);
1102  const unsigned int dof_index2_2(
1103  cell_type2_offset2 + shift_ijk);
1104
1105  for (unsigned int d = 0; d < dim; ++d)
1106  {
1108  [d] =
1109  data
1112  [d] =
1113  -1.0 *
1114  data
1117  [d] =
1118  data
1120
1122  [d] =
1123  data
1126  [d] =
1127  -1.0 *
1128  data
1131  [d] =
1132  -1.0 *
1133  data
1135  }
1136  }
1137  // Type 3: (note we re-use k and j for
1138  // convenience):
1139  const unsigned int shift_ij(
1140  j + shift_j); // here we've subbed j for i,
1141  // k for j.
1142  const unsigned int dof_index3_1(
1143  cell_type3_offset1 + shift_ij);
1144  const unsigned int dof_index3_2(
1145  cell_type3_offset2 + shift_ij);
1146  const unsigned int dof_index3_3(
1147  cell_type3_offset3 + shift_ij);
1148  for (unsigned int d1 = 0; d1 < dim; ++d1)
1149  {
1150  for (unsigned int d2 = 0; d2 < dim; ++d2)
1151  {
1153  [d2] = 0.0;
1155  [d2] = 0.0;
1157  [d2] = 0.0;
1158  }
1159  }
1161  2.0 * polyy[j][1] * polyz[k][0];
1163  2.0 * polyy[j][0] * polyz[k][1];
1164
1166  2.0 * polyx[j][1] * polyz[k][0];
1168  2.0 * polyx[j][0] * polyz[k][1];
1169
1171  2.0 * polyx[j][1] * polyy[k][0];
1173  2.0 * polyx[j][0] * polyy[k][1];
1174  }
1175  }
1176  }
1177  if (flags & update_hessians)
1178  {
1179  for (unsigned int k = 0; k < degree; ++k)
1180  {
1181  const unsigned int shift_k(k * degree * degree);
1182  const unsigned int shift_j(
1183  k * degree); // Used below when subbing k for j
1184  // type 3
1185
1186  for (unsigned int j = 0; j < degree; ++j)
1187  {
1188  const unsigned int shift_jk(j * degree +
1189  shift_k);
1190  for (unsigned int i = 0; i < degree; ++i)
1191  {
1192  const unsigned int shift_ijk(shift_jk +
1193  i);
1194
1195  // Type 1:
1196  const unsigned int dof_index1(
1197  cell_type1_offset + shift_ijk);
1198
1199  data.shape_hessians[dof_index1][q][0][0]
1200  [0] =
1201  8.0 * polyx[i][3] * polyy[j][0] *
1202  polyz[k][0];
1203  data.shape_hessians[dof_index1][q][1][0]
1204  [0] =
1205  8.0 * polyx[i][2] * polyy[j][1] *
1206  polyz[k][0];
1207  data.shape_hessians[dof_index1][q][2][0]
1208  [0] =
1209  8.0 * polyx[i][2] * polyy[j][0] *
1210  polyz[k][1];
1211
1212  data.shape_hessians[dof_index1][q][0][1]
1213  [0] =
1214  data.shape_hessians[dof_index1][q][1][0]
1215  [0];
1216  data.shape_hessians[dof_index1][q][1][1]
1217  [0] =
1218  8.0 * polyx[i][1] * polyy[j][2] *
1219  polyz[k][0];
1220  data.shape_hessians[dof_index1][q][2][1]
1221  [0] =
1222  8.0 * polyx[i][1] * polyy[j][1] *
1223  polyz[k][1];
1224
1225  data.shape_hessians[dof_index1][q][0][2]
1226  [0] =
1227  data.shape_hessians[dof_index1][q][2][0]
1228  [0];
1229  data.shape_hessians[dof_index1][q][1][2]
1230  [0] =
1231  data.shape_hessians[dof_index1][q][2][1]
1232  [0];
1233  data.shape_hessians[dof_index1][q][2][2]
1234  [0] =
1235  8.0 * polyx[i][1] * polyy[j][0] *
1236  polyz[k][2];
1237
1238
1239  data.shape_hessians[dof_index1][q][0][0]
1240  [1] =
1241  data.shape_hessians[dof_index1][q][1][0]
1242  [0];
1243  data.shape_hessians[dof_index1][q][1][0]
1244  [1] =
1245  data.shape_hessians[dof_index1][q][1][1]
1246  [0];
1247  data.shape_hessians[dof_index1][q][2][0]
1248  [1] =
1249  data.shape_hessians[dof_index1][q][2][1]
1250  [0];
1251
1252  data.shape_hessians[dof_index1][q][0][1]
1253  [1] =
1254  data.shape_hessians[dof_index1][q][1][1]
1255  [0];
1256  data.shape_hessians[dof_index1][q][1][1]
1257  [1] =
1258  8.0 * polyx[i][0] * polyy[j][3] *
1259  polyz[k][0];
1260  data.shape_hessians[dof_index1][q][2][1]
1261  [1] =
1262  8.0 * polyx[i][0] * polyy[j][2] *
1263  polyz[k][1];
1264
1265  data.shape_hessians[dof_index1][q][0][2]
1266  [1] =
1267  data.shape_hessians[dof_index1][q][2][1]
1268  [0];
1269  data.shape_hessians[dof_index1][q][1][2]
1270  [1] =
1271  data.shape_hessians[dof_index1][q][2][1]
1272  [1];
1273  data.shape_hessians[dof_index1][q][2][2]
1274  [1] =
1275  8.0 * polyx[i][0] * polyy[j][1] *
1276  polyz[k][2];
1277
1278
1279  data.shape_hessians[dof_index1][q][0][0]
1280  [2] =
1281  data.shape_hessians[dof_index1][q][2][0]
1282  [0];
1283  data.shape_hessians[dof_index1][q][1][0]
1284  [2] =
1285  data.shape_hessians[dof_index1][q][2][1]
1286  [0];
1287  data.shape_hessians[dof_index1][q][2][0]
1288  [2] =
1289  data.shape_hessians[dof_index1][q][2][2]
1290  [0];
1291
1292  data.shape_hessians[dof_index1][q][0][1]
1293  [2] =
1294  data.shape_hessians[dof_index1][q][2][1]
1295  [0];
1296  data.shape_hessians[dof_index1][q][1][1]
1297  [2] =
1298  data.shape_hessians[dof_index1][q][2][1]
1299  [1];
1300  data.shape_hessians[dof_index1][q][2][1]
1301  [2] =
1302  data.shape_hessians[dof_index1][q][2][2]
1303  [1];
1304
1305  data.shape_hessians[dof_index1][q][0][2]
1306  [2] =
1307  data.shape_hessians[dof_index1][q][2][2]
1308  [0];
1309  data.shape_hessians[dof_index1][q][1][2]
1310  [2] =
1311  data.shape_hessians[dof_index1][q][2][2]
1312  [1];
1313  data.shape_hessians[dof_index1][q][2][2]
1314  [2] =
1315  8.0 * polyx[i][0] * polyy[j][0] *
1316  polyz[k][3];
1317
1318
1319  // Type 2:
1320  const unsigned int dof_index2_1(
1321  cell_type2_offset1 + shift_ijk);
1322  const unsigned int dof_index2_2(
1323  cell_type2_offset2 + shift_ijk);
1324
1325  for (unsigned int d1 = 0; d1 < dim; ++d1)
1326  {
1327  for (unsigned int d2 = 0; d2 < dim;
1328  ++d2)
1329  {
1330  data
1331  .shape_hessians[dof_index2_1][q]
1332  [0][d1][d2] =
1333  data
1334  .shape_hessians[dof_index1][q]
1335  [0][d1][d2];
1336  data
1337  .shape_hessians[dof_index2_1][q]
1338  [1][d1][d2] =
1339  -1.0 *
1340  data
1341  .shape_hessians[dof_index1][q]
1342  [1][d1][d2];
1343  data
1344  .shape_hessians[dof_index2_1][q]
1345  [2][d1][d2] =
1346  data
1347  .shape_hessians[dof_index1][q]
1348  [2][d1][d2];
1349
1350  data
1351  .shape_hessians[dof_index2_2][q]
1352  [0][d1][d2] =
1353  data
1354  .shape_hessians[dof_index1][q]
1355  [0][d1][d2];
1356  data
1357  .shape_hessians[dof_index2_2][q]
1358  [1][d1][d2] =
1359  -1.0 *
1360  data
1361  .shape_hessians[dof_index1][q]
1362  [1][d1][d2];
1363  data
1364  .shape_hessians[dof_index2_2][q]
1365  [2][d1][d2] =
1366  -1.0 *
1367  data
1368  .shape_hessians[dof_index1][q]
1369  [2][d1][d2];
1370  }
1371  }
1372  }
1373  // Type 3: (note we re-use k and j for
1374  // convenience):
1375  const unsigned int shift_ij(
1376  j + shift_j); // here we've subbed j for i,
1377  // k for j.
1378  const unsigned int dof_index3_1(
1379  cell_type3_offset1 + shift_ij);
1380  const unsigned int dof_index3_2(
1381  cell_type3_offset2 + shift_ij);
1382  const unsigned int dof_index3_3(
1383  cell_type3_offset3 + shift_ij);
1384  for (unsigned int d1 = 0; d1 < dim; ++d1)
1385  {
1386  for (unsigned int d2 = 0; d2 < dim; ++d2)
1387  {
1388  for (unsigned int d3 = 0; d3 < dim;
1389  ++d3)
1390  {
1391  data
1392  .shape_hessians[dof_index3_1][q]
1393  [d1][d2][d3] =
1394  0.0;
1395  data
1396  .shape_hessians[dof_index3_2][q]
1397  [d1][d2][d3] =
1398  0.0;
1399  data
1400  .shape_hessians[dof_index3_3][q]
1401  [d1][d2][d3] =
1402  0.0;
1403  }
1404  }
1405  }
1406  data
1407  .shape_hessians[dof_index3_1][q][0][1][1] =
1408  4.0 * polyy[j][2] * polyz[k][0];
1409  data
1410  .shape_hessians[dof_index3_1][q][0][1][2] =
1411  4.0 * polyy[j][1] * polyz[k][1];
1412
1413  data
1414  .shape_hessians[dof_index3_1][q][0][2][1] =
1415  data
1416  .shape_hessians[dof_index3_1][q][0][1][2];
1417  data
1418  .shape_hessians[dof_index3_1][q][0][2][2] =
1419  4.0 * polyy[j][0] * polyz[k][2];
1420
1421
1422  data
1423  .shape_hessians[dof_index3_2][q][1][0][0] =
1424  4.0 * polyx[j][2] * polyz[k][0];
1425  data
1426  .shape_hessians[dof_index3_2][q][1][0][2] =
1427  4.0 * polyx[j][1] * polyz[k][1];
1428
1429  data
1430  .shape_hessians[dof_index3_2][q][1][2][0] =
1431  data
1432  .shape_hessians[dof_index3_2][q][1][0][2];
1433  data
1434  .shape_hessians[dof_index3_2][q][1][2][2] =
1435  4.0 * polyx[j][0] * polyz[k][2];
1436
1437
1438  data
1439  .shape_hessians[dof_index3_3][q][2][0][0] =
1440  4.0 * polyx[j][2] * polyy[k][0];
1441  data
1442  .shape_hessians[dof_index3_3][q][2][0][1] =
1443  4.0 * polyx[j][1] * polyy[k][1];
1444
1445  data
1446  .shape_hessians[dof_index3_3][q][2][1][0] =
1447  data
1448  .shape_hessians[dof_index3_3][q][2][0][1];
1449  data
1450  .shape_hessians[dof_index3_3][q][2][1][1] =
1451  4.0 * polyx[j][0] * polyy[k][2];
1452  }
1453  }
1454  }
1455  }
1456  }
1457  }
1458  break;
1459  }
1460  default:
1461  {
1462  Assert(false, ExcNotImplemented());
1463  }
1464  }
1465  return data_ptr;
1466 }
1467
1468
1469
1470 template <int dim, int spacedim>
1471 void
1473  const typename Triangulation<dim, dim>::cell_iterator &cell,
1475  const InternalData & fe_data) const
1476 {
1477  // This function handles the cell-dependent construction of the EDGE-based
1478  // shape functions.
1479  //
1480  // Note it will handle both 2D and 3D, in 2D, the edges are faces, but we
1481  // handle them here.
1482  //
1483  // It will fill in the missing parts of fe_data which were not possible to
1484  // fill in the get_data routine, with respect to the edge-based shape
1485  // functions.
1486  //
1487  // It should be called by the fill_fe_*_values routines in order to complete
1488  // the basis set at quadrature points on the current cell for each edge.
1489
1490  const UpdateFlags flags(fe_data.update_each);
1491  const unsigned int n_q_points = quadrature.size();
1492
1493  Assert(!(flags & update_values) ||
1494  fe_data.shape_values.size() == this->n_dofs_per_cell(),
1495  ExcDimensionMismatch(fe_data.shape_values.size(),
1496  this->n_dofs_per_cell()));
1497  Assert(!(flags & update_values) ||
1498  fe_data.shape_values[0].size() == n_q_points,
1499  ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
1500
1501  // Useful constants:
1502  const unsigned int degree(
1503  this->degree -
1504  1); // Note: constructor takes input degree + 1, so need to knock 1 off.
1505
1506  // Useful geometry info:
1507  const unsigned int vertices_per_line(2);
1508  const unsigned int lines_per_cell(GeometryInfo<dim>::lines_per_cell);
1509
1510  // Calculate edge orderings:
1511  std::vector<std::vector<unsigned int>> edge_order(
1512  lines_per_cell, std::vector<unsigned int>(vertices_per_line));
1513
1514
1515  switch (dim)
1516  {
1517  case 2:
1518  {
1520  {
1521  // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1},
1522  // e^{m}_{2}] e1 = higher global numbering of the two local
1523  // vertices e2 = lower global numbering of the two local vertices
1524  std::vector<int> edge_sign(lines_per_cell);
1525  for (unsigned int m = 0; m < lines_per_cell; ++m)
1526  {
1527  unsigned int v0_loc =
1529  unsigned int v1_loc =
1531  unsigned int v0_glob = cell->vertex_index(v0_loc);
1532  unsigned int v1_glob = cell->vertex_index(v1_loc);
1533
1534  if (v0_glob > v1_glob)
1535  {
1536  // Opposite to global numbering on our reference element
1537  edge_sign[m] = -1.0;
1538  }
1539  else
1540  {
1541  // Aligns with global numbering on our reference element.
1542  edge_sign[m] = 1.0;
1543  }
1544  }
1545
1546  // Define \sigma_{m} = sigma_{e^{m}_{2}} - sigma_{e^{m}_{1}}
1547  // \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
1548  //
1549  // To help things, in fe_data, we have precomputed (sigma_{i} -
1550  // sigma_{j}) and (lambda_{i} + lambda_{j}) for 0<= i,j <
1551  // lines_per_cell.
1552  //
1553  // There are two types:
1554  // - lower order (1 per edge, m):
1555  // \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
1556  //
1557  // - higher order (degree per edge, m):
1558  // \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
1559  //
1560  // NOTE: sigma_{m} and lambda_{m} are either a function of x OR
1561  // y
1562  // and if sigma is of x, then lambda is of y, and vice
1563  // versa. This means that grad(\sigma) requires
1564  // multiplication by d(sigma)/dx_{i} for the i^th comp of
1565  // grad(sigma) and similarly when taking derivatives of
1566  // lambda.
1567  //
1568  // First handle the lowest order edges (dofs 0 to 3)
1569  // 0 and 1 are the edges in the y dir. (sigma is function of y,
1570  // lambda is function of x). 2 and 3 are the edges in the x dir.
1571  // (sigma is function of x, lambda is function of y).
1572  //
1573  // More more info: see GeometryInfo for picture of the standard
1574  // element.
1575  //
1576  // Fill edge-based points:
1577  // std::vector<std::vector< Point<dim> > >
1578  // edge_points(lines_per_cell, std::vector<Point<dim>>
1579  // (n_q_points));
1580
1581  std::vector<std::vector<double>> edge_sigma_values(
1582  fe_data.edge_sigma_values);
1585
1586  std::vector<std::vector<double>> edge_lambda_values(
1587  fe_data.edge_lambda_values);
1590
1591  // Adjust the edge_sigma_* for the current cell:
1592  for (unsigned int m = 0; m < lines_per_cell; ++m)
1593  {
1594  std::transform(edge_sigma_values[m].begin(),
1595  edge_sigma_values[m].end(),
1596  edge_sigma_values[m].begin(),
1597  [&](const double edge_sigma_value) {
1598  return edge_sign[m] * edge_sigma_value;
1599  });
1600
1606  });
1607  }
1608
1609  // If we want to generate shape gradients then we need second
1610  // derivatives of the 1d polynomials, but only first derivatives
1611  // for the shape values.
1612  const unsigned int poly_length =
1613  (flags & update_hessians) ?
1614  4 :
1615  ((flags & update_gradients) ? 3 : 2);
1616
1617
1618  for (unsigned int m = 0; m < lines_per_cell; ++m)
1619  {
1620  const unsigned int shift_m(m * this->n_dofs_per_line());
1621  for (unsigned int q = 0; q < n_q_points; ++q)
1622  {
1623  // Only compute 1d polynomials if degree>0.
1624  std::vector<std::vector<double>> poly(
1625  degree, std::vector<double>(poly_length));
1626  for (unsigned int i = 1; i < degree + 1; ++i)
1627  {
1628  // Compute all required 1d polynomials and their
1629  // derivatives, starting at degree 2. e.g. to access
1630  // L'_{i+2}(edge_sigma) use polyx[i][1].
1631  IntegratedLegendrePolynomials[i + 1].value(
1632  edge_sigma_values[m][q], poly[i - 1]);
1633  }
1634  if (flags & update_values)
1635  {
1636  // Lowest order edge shape functions:
1637  for (unsigned int d = 0; d < dim; ++d)
1638  {
1639  fe_data.shape_values[shift_m][q][d] =
1641  edge_lambda_values[m][q];
1642  }
1643  // Higher order edge shape functions:
1644  for (unsigned int i = 1; i < degree + 1; ++i)
1645  {
1646  const unsigned int poly_index = i - 1;
1647  const unsigned int dof_index(i + shift_m);
1648  for (unsigned int d = 0; d < dim; ++d)
1649  {
1650  fe_data.shape_values[dof_index][q][d] =
1652  poly[poly_index][1] *
1653  edge_lambda_values[m][q] +
1654  poly[poly_index][0] *
1656  }
1657  }
1658  }
1660  {
1661  // Lowest order edge shape functions:
1662  for (unsigned int d1 = 0; d1 < dim; ++d1)
1663  {
1664  for (unsigned int d2 = 0; d2 < dim; ++d2)
1665  {
1666  // Note: gradient is constant for a given
1667  // edge.
1671  }
1672  }
1673  // Higher order edge shape functions:
1674  for (unsigned int i = 1; i < degree + 1; ++i)
1675  {
1676  const unsigned int poly_index = i - 1;
1677  const unsigned int dof_index(i + shift_m);
1678
1682  edge_lambda_values[m][q] * poly[poly_index][2];
1683
1689  poly[poly_index][1];
1690
1693
1697  edge_lambda_values[m][q] * poly[poly_index][2];
1698  }
1699  }
1700  if (flags & update_hessians)
1701  {
1702  // Lowest order edge shape function
1703  for (unsigned int d1 = 0; d1 < dim; ++d1)
1704  {
1705  for (unsigned int d2 = 0; d2 < dim; ++d2)
1706  {
1707  for (unsigned int d3 = 0; d3 < dim; ++d3)
1708  {
1709  fe_data.shape_hessians[shift_m][q][d1][d2]
1710  [d3] = 0;
1711  }
1712  }
1713  }
1714
1715  // Higher order edge shape function
1716  for (unsigned int i = 0; i < degree; ++i)
1717  {
1718  const unsigned int dof_index(i + 1 + shift_m);
1719
1720  for (unsigned int d1 = 0; d1 < dim; ++d1)
1721  {
1722  for (unsigned int d2 = 0; d2 < dim; ++d2)
1723  {
1724  for (unsigned int d3 = 0; d3 < dim; ++d3)
1725  {
1726  fe_data.shape_hessians[dof_index][q]
1727  [d1][d2][d3] =
1731  poly[i][3] *
1732  edge_lambda_values[m][q] +
1733  poly[i][2] *
1743  }
1744  }
1745  }
1746  }
1747  }
1748  }
1749  }
1750  }
1751  break;
1752  }
1753  case 3:
1754  {
1756  {
1757  // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1},
1758  // e^{m}_{2}] e1 = higher global numbering of the two local
1759  // vertices e2 = lower global numbering of the two local vertices
1760  std::vector<int> edge_sign(lines_per_cell);
1761  for (unsigned int m = 0; m < lines_per_cell; ++m)
1762  {
1763  unsigned int v0_loc =
1765  unsigned int v1_loc =
1767  unsigned int v0_glob = cell->vertex_index(v0_loc);
1768  unsigned int v1_glob = cell->vertex_index(v1_loc);
1769
1770  if (v0_glob > v1_glob)
1771  {
1772  // Opposite to global numbering on our reference element
1773  edge_sign[m] = -1.0;
1774  }
1775  else
1776  {
1777  // Aligns with global numbering on our reference element.
1778  edge_sign[m] = 1.0;
1779  }
1780  }
1781
1782  // Define \sigma_{m} = sigma_{e^{m}_{1}} - sigma_{e^{m}_{2}}
1783  // \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
1784  //
1785  // To help things, in fe_data, we have precomputed (sigma_{i} -
1786  // sigma_{j}) and (lambda_{i} + lambda_{j}) for 0<= i,j <
1787  // lines_per_cell.
1788  //
1789  // There are two types:
1790  // - lower order (1 per edge, m):
1791  // \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
1792  //
1793  // - higher order (degree per edge, m):
1794  // \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
1795  //
1796  // NOTE: In the ref cell, sigma_{m} is a function of x OR y OR Z
1797  // and lambda_{m} a function of the remaining co-ords.
1798  // for example, if sigma is of x, then lambda is of y AND
1799  // z, and so on. This means that grad(\sigma) requires
1800  // multiplication by d(sigma)/dx_{i} for the i^th comp of
1801  // grad(sigma) and similarly when taking derivatives of
1802  // lambda.
1803  //
1804  // First handle the lowest order edges (dofs 0 to 11)
1805  // 0 and 1 are the edges in the y dir at z=0. (sigma is a fn of y,
1806  // lambda is a fn of x & z). 2 and 3 are the edges in the x dir at
1807  // z=0. (sigma is a fn of x, lambda is a fn of y & z). 4 and 5 are
1808  // the edges in the y dir at z=1. (sigma is a fn of y, lambda is a
1809  // fn of x & z). 6 and 7 are the edges in the x dir at z=1. (sigma
1810  // is a fn of x, lambda is a fn of y & z). 8 and 9 are the edges
1811  // in the z dir at y=0. (sigma is a fn of z, lambda is a fn of x &
1812  // y). 10 and 11 are the edges in the z dir at y=1. (sigma is a fn
1813  // of z, lambda is a fn of x & y).
1814  //
1815  // For more info: see GeometryInfo for picture of the standard
1816  // element.
1817
1818  // Copy over required edge-based data:
1819  std::vector<std::vector<double>> edge_sigma_values(
1820  fe_data.edge_sigma_values);
1821  std::vector<std::vector<double>> edge_lambda_values(
1822  fe_data.edge_lambda_values);
1827  std::vector<std::vector<std::vector<double>>>
1829
1830  // Adjust the edge_sigma_* for the current cell:
1831  for (unsigned int m = 0; m < lines_per_cell; ++m)
1832  {
1833  std::transform(edge_sigma_values[m].begin(),
1834  edge_sigma_values[m].end(),
1835  edge_sigma_values[m].begin(),
1836  [&](const double edge_sigma_value) {
1837  return edge_sign[m] * edge_sigma_value;
1838  });
1844  });
1845  }
1846
1847  // Now calculate the edge-based shape functions:
1848  // If we want to generate shape gradients then we need second
1849  // derivatives of the 1d polynomials, but only first derivatives
1850  // for the shape values.
1851  const unsigned int poly_length =
1852  (flags & update_hessians) ?
1853  4 :
1854  ((flags & update_gradients) ? 3 : 2);
1855
1856  std::vector<std::vector<double>> poly(
1857  degree, std::vector<double>(poly_length));
1858  for (unsigned int m = 0; m < lines_per_cell; ++m)
1859  {
1860  const unsigned int shift_m(m * this->n_dofs_per_line());
1861  for (unsigned int q = 0; q < n_q_points; ++q)
1862  {
1863  // precompute values of all 1d polynomials required:
1864  // for example poly[i][1] = L'_{i+2}(edge_sigma_values)
1865  if (degree > 0)
1866  {
1867  for (unsigned int i = 0; i < degree; ++i)
1868  {
1869  IntegratedLegendrePolynomials[i + 2].value(
1870  edge_sigma_values[m][q], poly[i]);
1871  }
1872  }
1873  if (flags & update_values)
1874  {
1875  // Lowest order edge shape functions:
1876  for (unsigned int d = 0; d < dim; ++d)
1877  {
1878  fe_data.shape_values[shift_m][q][d] =
1880  edge_lambda_values[m][q];
1881  }
1882  // Higher order edge shape functions
1883  if (degree > 0)
1884  {
1885  for (unsigned int i = 0; i < degree; ++i)
1886  {
1887  const unsigned int dof_index(i + 1 + shift_m);
1888  for (unsigned int d = 0; d < dim; ++d)
1889  {
1890  fe_data.shape_values[dof_index][q][d] =
1892  edge_lambda_values[m][q] +
1894  }
1895  }
1896  }
1897  }
1899  {
1900  // Lowest order edge shape functions:
1901  for (unsigned int d1 = 0; d1 < dim; ++d1)
1902  {
1903  for (unsigned int d2 = 0; d2 < dim; ++d2)
1904  {
1908  }
1909  }
1910  // Higher order edge shape functions
1911  if (degree > 0)
1912  {
1913  for (unsigned int i = 0; i < degree; ++i)
1914  {
1915  const unsigned int dof_index(i + 1 + shift_m);
1916
1917  for (unsigned int d1 = 0; d1 < dim; ++d1)
1918  {
1919  for (unsigned int d2 = 0; d2 < dim; ++d2)
1920  {
1921  fe_data
1925  poly[i][2] *
1926  edge_lambda_values[m][q] +
1931  poly[i][1] +
1933  [d2] *
1934  poly[i][0];
1935  }
1936  }
1937  }
1938  }
1939  }
1940  if (flags & update_hessians)
1941  {
1942  // Lowest order edge shape functions:
1943  for (unsigned int d1 = 0; d1 < dim; ++d1)
1944  {
1945  for (unsigned int d2 = 0; d2 < dim; ++d2)
1946  {
1947  for (unsigned int d3 = 0; d3 < dim; ++d3)
1948  {
1949  fe_data.shape_hessians[shift_m][q][d1][d2]
1950  [d3] =
1953  }
1954  }
1955  }
1956
1957  // Higher order edge shape functions
1958  if (degree > 0)
1959  {
1960  for (unsigned int i = 0; i < degree; ++i)
1961  {
1962  const unsigned int dof_index(i + 1 + shift_m);
1963
1964  for (unsigned int d1 = 0; d1 < dim; ++d1)
1965  {
1966  for (unsigned int d2 = 0; d2 < dim; ++d2)
1967  {
1968  for (unsigned int d3 = 0; d3 < dim;
1969  ++d3)
1970  {
1971  fe_data
1972  .shape_hessians[dof_index][q]
1973  [d1][d2][d3] =
1977  poly[i][3] *
1978  edge_lambda_values[m][q] +
1979  poly[i][2] *
1983  [d3] +
1987  [d2] +
1991  [d1]) +
1992  poly[i][1] *
1995  [m][d3][d2] +
1998  [m][d3][d1] +
2001  [m][d2][d1]);
2002  }
2003  }
2004  }
2005  }
2006  }
2007  }
2008  }
2009  }
2010  }
2011  break;
2012  }
2013  default:
2014  {
2015  Assert(false, ExcNotImplemented());
2016  }
2017  }
2018 }
2019
2020
2021
2022 template <int dim, int spacedim>
2023 void
2025  const typename Triangulation<dim, dim>::cell_iterator &cell,
2027  const InternalData & fe_data) const
2028 {
2029  // This function handles the cell-dependent construction of the FACE-based
2030  // shape functions.
2031  //
2032  // Note that it should only be called in 3D.
2033  Assert(dim == 3, ExcDimensionMismatch(dim, 3));
2034  //
2035  // It will fill in the missing parts of fe_data which were not possible to
2036  // fill in the get_data routine, with respect to face-based shape functions.
2037  //
2038  // It should be called by the fill_fe_*_values routines in order to complete
2039  // the basis set at quadrature points on the current cell for each face.
2040
2041  // Useful constants:
2042  const unsigned int degree(
2043  this->degree -
2044  1); // Note: constructor takes input degree + 1, so need to knock 1 off.
2045
2046  // Do nothing if FE degree is 0.
2047  if (degree > 0)
2048  {
2049  const UpdateFlags flags(fe_data.update_each);
2050
2052  {
2053  const unsigned int n_q_points = quadrature.size();
2054
2055  Assert(!(flags & update_values) ||
2056  fe_data.shape_values.size() == this->n_dofs_per_cell(),
2057  ExcDimensionMismatch(fe_data.shape_values.size(),
2058  this->n_dofs_per_cell()));
2059  Assert(!(flags & update_values) ||
2060  fe_data.shape_values[0].size() == n_q_points,
2061  ExcDimensionMismatch(fe_data.shape_values[0].size(),
2062  n_q_points));
2063
2064  // Useful geometry info:
2065  const unsigned int vertices_per_face(
2067  const unsigned int faces_per_cell(GeometryInfo<dim>::faces_per_cell);
2068
2069  // DoF info:
2070  const unsigned int n_line_dofs =
2071  this->n_dofs_per_line() * GeometryInfo<dim>::lines_per_cell;
2072
2073  // First we find the global face orientations on the current cell.
2074  std::vector<std::vector<unsigned int>> face_orientation(
2075  faces_per_cell, std::vector<unsigned int>(vertices_per_face));
2076
2077  const unsigned int
2078  vertex_opposite_on_face[GeometryInfo<3>::vertices_per_face] = {3,
2079  2,
2080  1,
2081  0};
2082
2083  const unsigned int
2085  {1, 2}, {0, 3}, {0, 3}, {1, 2}};
2086
2087  for (unsigned int m = 0; m < faces_per_cell; ++m)
2088  {
2089  // Find the local vertex on this face with the highest global
2090  // numbering. This is f^m_0.
2091  unsigned int current_max = 0;
2092  unsigned int current_glob = cell->vertex_index(
2094  for (unsigned int v = 1; v < vertices_per_face; ++v)
2095  {
2096  if (current_glob <
2097  cell->vertex_index(
2099  {
2100  current_max = v;
2101  current_glob = cell->vertex_index(
2103  }
2104  }
2105  face_orientation[m][0] =
2107
2108  // f^m_2 is the vertex opposite f^m_0.
2109  face_orientation[m][2] = GeometryInfo<dim>::face_to_cell_vertices(
2110  m, vertex_opposite_on_face[current_max]);
2111
2112  // Finally, f^m_1 is the vertex with the greater global numbering
2113  // of the remaining two local vertices. Then, f^m_3 is the other.
2114  if (cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(
2116  cell->vertex_index(GeometryInfo<dim>::face_to_cell_vertices(
2118  {
2119  face_orientation[m][1] =
2122  face_orientation[m][3] =
2125  }
2126  else
2127  {
2128  face_orientation[m][1] =
2131  face_orientation[m][3] =
2134  }
2135  }
2136
2137  // Now we know the face orientation on the current cell, we can
2138  // generate the parameterisation:
2139  std::vector<std::vector<double>> face_xi_values(
2140  faces_per_cell, std::vector<double>(n_q_points));
2142  faces_per_cell, std::vector<double>(dim));
2143  std::vector<std::vector<double>> face_eta_values(
2144  faces_per_cell, std::vector<double>(n_q_points));
2146  faces_per_cell, std::vector<double>(dim));
2147
2148  std::vector<std::vector<double>> face_lambda_values(
2149  faces_per_cell, std::vector<double>(n_q_points));
2151  faces_per_cell, std::vector<double>(dim));
2152  for (unsigned int m = 0; m < faces_per_cell; ++m)
2153  {
2154  for (unsigned int q = 0; q < n_q_points; ++q)
2155  {
2156  face_xi_values[m][q] =
2157  fe_data.sigma_imj_values[q][face_orientation[m][0]]
2158  [face_orientation[m][1]];
2159  face_eta_values[m][q] =
2160  fe_data.sigma_imj_values[q][face_orientation[m][0]]
2161  [face_orientation[m][3]];
2162  face_lambda_values[m][q] = fe_data.face_lambda_values[m][q];
2163  }
2164  for (unsigned int d = 0; d < dim; ++d)
2165  {
2168  [face_orientation[m][1]][d];
2171  [face_orientation[m][3]][d];
2172
2174  }
2175  }
2176  // Now can generate the basis
2177  const unsigned int poly_length =
2178  (flags & update_hessians) ? 4 :
2179  ((flags & update_gradients) ? 3 : 2);
2180
2181
2182  std::vector<std::vector<double>> polyxi(
2183  degree, std::vector<double>(poly_length));
2184  std::vector<std::vector<double>> polyeta(
2185  degree, std::vector<double>(poly_length));
2186
2187  // Loop through quad points:
2188  for (unsigned int m = 0; m < faces_per_cell; ++m)
2189  {
2190  // we assume that all quads have the same number of dofs
2191  const unsigned int shift_m(m * this->n_dofs_per_quad(0));
2192  // Calculate the offsets for each face-based shape function:
2193  //
2195  // \phi^{F_m,1}_{ij} = \nabla( L_{i+2}(\xi_{F_{m}})
2196  // L_{j+2}(\eta_{F_{m}}) \lambda_{F_{m}} )
2197  //
2198  // 0 <= i,j < degree (in a group of degree*degree)
2199  const unsigned int face_type1_offset(n_line_dofs + shift_m);
2200  // Type-2:
2201  //
2202  // \phi^{F_m,2}_{ij} = ( L'_{i+2}(\xi_{F_{m}})
2203  // L_{j+2}(\eta_{F_{m}}) \nabla\xi_{F_{m}}
2204  // - L_{i+2}(\xi_{F_{m}})
2205  // L'_{j+2}(\eta_{F_{m}}) \nabla\eta_{F_{m}}
2206  // ) \lambda_{F_{m}}
2207  //
2208  // 0 <= i,j < degree (in a group of degree*degree)
2209  const unsigned int face_type2_offset(face_type1_offset +
2210  degree * degree);
2211  // Type-3:
2212  //
2213  // \phi^{F_m,3}_{i} = L_{i+2}(\eta_{F_{m}}) \lambda_{F_{m}}
2214  // \nabla\xi_{F_{m}} \phi^{F_m,3}_{i+p} = L_{i+2}(\xi_{F_{m}})
2215  // \lambda_{F_{m}} \nabla\eta_{F_{m}}
2216  //
2217  // 0 <= i < degree.
2218  //
2219  // here we order so that all of subtype 1 comes first, then
2220  // subtype 2.
2221  const unsigned int face_type3_offset1(face_type2_offset +
2222  degree * degree);
2223  const unsigned int face_type3_offset2(face_type3_offset1 +
2224  degree);
2225
2226  // Loop over all faces:
2227  for (unsigned int q = 0; q < n_q_points; ++q)
2228  {
2229  // pre-compute values & required derivatives at this quad
2230  // point: polyxi = L_{i+2}(\xi_{F_{m}}), polyeta =
2231  // L_{j+2}(\eta_{F_{m}}),
2232  //
2233  // each polypoint[k][d], contains the dth derivative of
2234  // L_{k+2} at the point \xi or \eta. Note that this doesn't
2235  // include the derivative of xi/eta via the chain rule.
2236  for (unsigned int i = 0; i < degree; ++i)
2237  {
2238  // compute all required 1d polynomials:
2239  IntegratedLegendrePolynomials[i + 2].value(
2240  face_xi_values[m][q], polyxi[i]);
2241  IntegratedLegendrePolynomials[i + 2].value(
2242  face_eta_values[m][q], polyeta[i]);
2243  }
2244  // Now use these to compute the shape functions:
2245  if (flags & update_values)
2246  {
2247  for (unsigned int j = 0; j < degree; ++j)
2248  {
2249  const unsigned int shift_j(j * degree);
2250  for (unsigned int i = 0; i < degree; ++i)
2251  {
2252  const unsigned int shift_ij(shift_j + i);
2253  // Type 1:
2254  const unsigned int dof_index1(face_type1_offset +
2255  shift_ij);
2256  for (unsigned int d = 0; d < dim; ++d)
2257  {
2258  fe_data.shape_values[dof_index1][q][d] =
2260  polyeta[j][0] +
2262  polyeta[j][1]) *
2263  face_lambda_values[m][q] +
2265  polyeta[j][0];
2266  }
2267  // Type 2:
2268  const unsigned int dof_index2(face_type2_offset +
2269  shift_ij);
2270  for (unsigned int d = 0; d < dim; ++d)
2271  {
2272  fe_data.shape_values[dof_index2][q][d] =
2274  polyeta[j][0] -
2276  polyeta[j][1]) *
2277  face_lambda_values[m][q];
2278  }
2279  }
2280  // Type 3:
2281  const unsigned int dof_index3_1(face_type3_offset1 +
2282  j);
2283  const unsigned int dof_index3_2(face_type3_offset2 +
2284  j);
2285  for (unsigned int d = 0; d < dim; ++d)
2286  {
2287  fe_data.shape_values[dof_index3_1][q][d] =
2289  face_lambda_values[m][q];
2290
2291  fe_data.shape_values[dof_index3_2][q][d] =
2293  face_lambda_values[m][q];
2294  }
2295  }
2296  }
2298  {
2299  for (unsigned int j = 0; j < degree; ++j)
2300  {
2301  const unsigned int shift_j(j * degree);
2302  for (unsigned int i = 0; i < degree; ++i)
2303  {
2304  const unsigned int shift_ij(shift_j + i);
2305  // Type 1:
2306  const unsigned int dof_index1(face_type1_offset +
2307  shift_ij);
2308  for (unsigned int d1 = 0; d1 < dim; ++d1)
2309  {
2310  for (unsigned int d2 = 0; d2 < dim; ++d2)
2311  {
2312  fe_data
2316  polyeta[j][0] +
2321  polyxi[i][1] * polyeta[j][1] +
2324  polyxi[i][0] * polyeta[j][2]) *
2325  face_lambda_values[m][q] +
2327  polyeta[j][0] +
2329  polyeta[j][1]) *
2332  polyeta[j][0] +
2334  polyeta[j][1]) *
2336  }
2337  }
2338  // Type 2:
2339  const unsigned int dof_index2(face_type2_offset +
2340  shift_ij);
2341  for (unsigned int d1 = 0; d1 < dim; ++d1)
2342  {
2343  for (unsigned int d2 = 0; d2 < dim; ++d2)
2344  {
2345  fe_data
2349  polyeta[j][0] +
2354  polyxi[i][1] * polyeta[j][1] -
2357  polyxi[i][0] * polyeta[j][2]) *
2358  face_lambda_values[m][q] +
2360  polyeta[j][0] -
2362  polyeta[j][1]) *
2364  }
2365  }
2366  }
2367  // Type 3:
2368  const unsigned int dof_index3_1(face_type3_offset1 +
2369  j);
2370  const unsigned int dof_index3_2(face_type3_offset2 +
2371  j);
2372  for (unsigned int d1 = 0; d1 < dim; ++d1)
2373  {
2374  for (unsigned int d2 = 0; d2 < dim; ++d2)
2375  {
2379  face_lambda_values[m][q] +
2381
2385  face_lambda_values[m][q] +
2387  }
2388  }
2389  }
2390  }
2391  if (flags & update_hessians)
2392  {
2393  for (unsigned int j = 0; j < degree; ++j)
2394  {
2395  const unsigned int shift_j(j * degree);
2396  for (unsigned int i = 0; i < degree; ++i)
2397  {
2398  const unsigned int shift_ij(shift_j + i);
2399
2400  // Type 1:
2401  const unsigned int dof_index1(face_type1_offset +
2402  shift_ij);
2403  for (unsigned int d1 = 0; d1 < dim; ++d1)
2404  {
2405  for (unsigned int d2 = 0; d2 < dim; ++d2)
2406  {
2407  for (unsigned int d3 = 0; d3 < dim; ++d3)
2408  {
2409  fe_data.shape_hessians[dof_index1][q]
2410  [d1][d2][d3] =
2411  polyxi[i][1] *
2414  (polyeta[j][2] *
2416  face_lambda_values[m][q] +
2417  polyeta[j][1] *
2419  polyeta[j][1] *
2422  polyxi[i][0] *
2423  (polyeta[j][3] *
2427  face_lambda_values[m][q] +
2428  polyeta[j][2] *
2435  [d2] +
2438  [m][d1]))) +
2439  (polyxi[i][1] * polyeta[j][1] *
2441  polyxi[i][2] * polyeta[j][0] *
2448  (polyxi[i][2] * polyeta[j][0] *
2451  polyxi[i][1] * polyeta[j][1] *
2456  face_lambda_values[m][q] *
2457  (polyxi[i][3] * polyeta[j][0] *
2461  polyxi[i][1] * polyeta[j][2] *
2467  polyxi[i][2] * polyeta[j][1] *
2476  }
2477  }
2478  }
2479
2480  // Type 2:
2481  const unsigned int dof_index2(face_type2_offset +
2482  shift_ij);
2483  for (unsigned int d1 = 0; d1 < dim; ++d1)
2484  {
2485  for (unsigned int d2 = 0; d2 < dim; ++d2)
2486  {
2487  for (unsigned int d3 = 0; d3 < dim; ++d3)
2488  {
2489  fe_data.shape_hessians[dof_index2][q]
2490  [d1][d2][d3] =
2492  (polyxi[i][1] * polyeta[j][1] *
2497  polyxi[i][2] * polyeta[j][0] *
2502  face_lambda_values[m][q] *
2504  (polyxi[i][1] *
2505  polyeta[j][2] *
2507  polyxi[i][2] *
2508  polyeta[j][1] *
2511  (polyxi[i][2] *
2512  polyeta[j][1] *
2514  polyxi[i][3] *
2515  polyeta[j][0] *
2517  polyxi[i][0] *
2520  (polyeta[j][3] *
2522  face_lambda_values[m][q] +
2523  polyeta[j][2] *
2525  polyeta[j][2] *
2529  (polyxi[i][1] *
2531  (polyeta[j][2] *
2533  face_lambda_values[m][q] +
2534  polyeta[j][1] *
2537  (polyxi[i][1] *
2538  (polyeta[j][2] *
2540  face_lambda_values[m]
2541  [q] +
2542  polyeta[j][1] *
2544  [d3]) +
2545  polyxi[i][2] * polyeta[j][1] *
2547  face_lambda_values[m][q]));
2548  }
2549  }
2550  }
2551  }
2552  // Type 3:
2553  const unsigned int dof_index3_1(face_type3_offset1 +
2554  j);
2555  const unsigned int dof_index3_2(face_type3_offset2 +
2556  j);
2557  for (unsigned int d1 = 0; d1 < dim; ++d1)
2558  {
2559  for (unsigned int d2 = 0; d2 < dim; ++d2)
2560  {
2561  for (unsigned int d3 = 0; d3 < dim; ++d3)
2562  {
2563  fe_data.shape_hessians[dof_index3_1][q]
2564  [d1][d2][d3] =
2567  (polyeta[j][2] *
2569  face_lambda_values[m][q] +
2570  polyeta[j][1] *
2573  polyeta[j][1] *
2575
2576  fe_data.shape_hessians[dof_index3_2][q]
2577  [d1][d2][d3] =
2580  (polyxi[j][2] *
2582  face_lambda_values[m][q] +
2583  polyxi[j][1] *
2587  }
2588  }
2589  }
2590  }
2591  }
2592  }
2593  }
2594  }
2595  }
2596 }
2597
2598
2599
2600 template <int dim, int spacedim>
2601 void
2603  const typename Triangulation<dim, dim>::cell_iterator &cell,
2604  const CellSimilarity::Similarity /*cell_similarity*/,
2606  const Mapping<dim, dim> & mapping,
2607  const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
2608  const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2609  & mapping_data,
2610  const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
2612  &data) const
2613 {
2614  // Convert to the correct internal data class for this FE class.
2615  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
2616  ExcInternalError());
2617  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
2618
2619  // Now update the edge-based DoFs, which depend on the cell.
2620  // This will fill in the missing items in the InternalData
2621  // (fe_internal/fe_data) which was not filled in by get_data.
2623  if (dim == 3 && this->degree > 1)
2624  {
2626  }
2627
2628  const UpdateFlags flags(fe_data.update_each);
2629  const unsigned int n_q_points = quadrature.size();
2630
2631  Assert(!(flags & update_values) ||
2632  fe_data.shape_values.size() == this->n_dofs_per_cell(),
2633  ExcDimensionMismatch(fe_data.shape_values.size(),
2634  this->n_dofs_per_cell()));
2635  Assert(!(flags & update_values) ||
2636  fe_data.shape_values[0].size() == n_q_points,
2637  ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
2638
2639  if (flags & update_values)
2640  {
2641  // Now have all shape_values stored on the reference cell.
2642  // Must now transform to the physical cell.
2643  std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
2644  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2645  {
2646  const unsigned int first =
2647  data.shape_function_to_row_table[dof * this->n_components() +
2648  this->get_nonzero_components(dof)
2649  .first_selected_component()];
2650
2651  mapping.transform(make_array_view(fe_data.shape_values[dof]),
2653  mapping_internal,
2654  make_array_view(transformed_shape_values));
2655  for (unsigned int q = 0; q < n_q_points; ++q)
2656  {
2657  for (unsigned int d = 0; d < dim; ++d)
2658  {
2659  data.shape_values(first + d, q) =
2660  transformed_shape_values[q][d];
2661  }
2662  }
2663  }
2664  }
2665
2667  {
2668  // Now have all shape_grads stored on the reference cell.
2669  // Must now transform to the physical cell.
2670  std::vector<Tensor<2, dim>> input(n_q_points);
2672  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2673  {
2674  for (unsigned int q = 0; q < n_q_points; ++q)
2675  {
2677  }
2678  mapping.transform(make_array_view(input),
2680  mapping_internal,
2682
2683  const unsigned int first =
2684  data.shape_function_to_row_table[dof * this->n_components() +
2685  this->get_nonzero_components(dof)
2686  .first_selected_component()];
2687
2688  for (unsigned int q = 0; q < n_q_points; ++q)
2689  {
2690  for (unsigned int d1 = 0; d1 < dim; ++d1)
2691  {
2692  for (unsigned int d2 = 0; d2 < dim; ++d2)
2693  {
2695  data.shape_values(first + d2, q) *
2697  }
2698  }
2699  }
2700
2701  for (unsigned int q = 0; q < n_q_points; ++q)
2702  {
2703  for (unsigned int d = 0; d < dim; ++d)
2704  {
2707  }
2708  }
2709  }
2710  }
2711
2712  if (flags & update_hessians)
2713  {
2714  // Now have all shape_grads stored on the reference cell.
2715  // Must now transform to the physical cell.
2716  std::vector<Tensor<3, dim>> input(n_q_points);
2717  std::vector<Tensor<3, dim>> transformed_shape_hessians(n_q_points);
2718  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2719  {
2720  for (unsigned int q = 0; q < n_q_points; ++q)
2721  {
2722  input[q] = fe_data.shape_hessians[dof][q];
2723  }
2724  mapping.transform(make_array_view(input),
2726  mapping_internal,
2727  make_array_view(transformed_shape_hessians));
2728
2729  const unsigned int first =
2730  data.shape_function_to_row_table[dof * this->n_components() +
2731  this->get_nonzero_components(dof)
2732  .first_selected_component()];
2733
2734  for (unsigned int q = 0; q < n_q_points; ++q)
2735  {
2736  for (unsigned int d1 = 0; d1 < dim; ++d1)
2737  {
2738  for (unsigned int d2 = 0; d2 < dim; ++d2)
2739  {
2740  for (unsigned int d3 = 0; d3 < dim; ++d3)
2741  {
2742  for (unsigned int d4 = 0; d4 < dim; ++d4)
2743  {
2744  transformed_shape_hessians[q][d1][d3][d4] -=
2745  (data.shape_values(first + d2, q) *
2746  mapping_data
2747  .jacobian_pushed_forward_2nd_derivatives
2748  [q][d2][d1][d3][d4]) +
2750  mapping_data
2752  [d4]) +
2754  mapping_data
2756  [d4]) +
2758  mapping_data
2760  [d1]);
2761  }
2762  }
2763  }
2764  }
2765  }
2766
2767  for (unsigned int q = 0; q < n_q_points; ++q)
2768  {
2769  for (unsigned int d = 0; d < dim; ++d)
2770  {
2771  data.shape_hessians[first + d][q] =
2772  transformed_shape_hessians[q][d];
2773  }
2774  }
2775  }
2776  }
2777 }
2778
2779
2780
2781 template <int dim, int spacedim>
2782 void
2784  const typename Triangulation<dim, dim>::cell_iterator &cell,
2785  const unsigned int face_no,
2786  const hp::QCollection<dim - 1> & quadrature,
2787  const Mapping<dim, dim> & mapping,
2788  const typename Mapping<dim, dim>::InternalDataBase & mapping_internal,
2789  const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2790  & mapping_data,
2791  const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
2793  &data) const
2794 {
2796
2797  // Note for future improvement:
2798  // We don't have the full quadrature - should use QProjector to create the 2D
2800  //
2801  // For now I am effectively generating all of the shape function vals/grads,
2802  // etc. On all quad points on all faces and then only using them for one face.
2803  // This is obviously inefficient. I should cache the cell number and cache
2804  // all of the shape_values/gradients etc and then reuse them for each face.
2805
2806  // convert data object to internal
2807  // data for this class. fails with
2808  // an exception if that is not
2809  // possible
2810  Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
2811  ExcInternalError());
2812  const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
2813
2814  // Now update the edge-based DoFs, which depend on the cell.
2815  // This will fill in the missing items in the InternalData
2816  // (fe_internal/fe_data) which was not filled in by get_data.
2817  fill_edge_values(cell,
2820  fe_data);
2821  if (dim == 3 && this->degree > 1)
2822  {
2823  fill_face_values(cell,
2826  fe_data);
2827  }
2828
2829  const UpdateFlags flags(fe_data.update_each);
2830  const unsigned int n_q_points = quadrature[0].size();
2831  const auto offset =
2833  face_no,
2834  cell->face_orientation(face_no),
2835  cell->face_flip(face_no),
2836  cell->face_rotation(face_no),
2837  n_q_points);
2838
2839  if (flags & update_values)
2840  {
2841  // Now have all shape_values stored on the reference cell.
2842  // Must now transform to the physical cell.
2843  std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
2844  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2845  {
2846  mapping.transform(make_array_view(fe_data.shape_values[dof],
2847  offset,
2848  n_q_points),
2850  mapping_internal,
2851  make_array_view(transformed_shape_values));
2852
2853  const unsigned int first =
2854  data.shape_function_to_row_table[dof * this->n_components() +
2855  this->get_nonzero_components(dof)
2856  .first_selected_component()];
2857
2858  for (unsigned int q = 0; q < n_q_points; ++q)
2859  {
2860  for (unsigned int d = 0; d < dim; ++d)
2861  {
2862  data.shape_values(first + d, q) =
2863  transformed_shape_values[q][d];
2864  }
2865  }
2866  }
2867  }
2869  {
2870  // Now have all shape_grads stored on the reference cell.
2871  // Must now transform to the physical cell.
2872  std::vector<Tensor<2, dim>> input(n_q_points);
2874  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2875  {
2876  for (unsigned int q = 0; q < n_q_points; ++q)
2877  {
2878  input[q] = fe_data.shape_grads[dof][offset + q];
2879  }
2880  mapping.transform(input,
2882  mapping_internal,
2884
2885  const unsigned int first =
2886  data.shape_function_to_row_table[dof * this->n_components() +
2887  this->get_nonzero_components(dof)
2888  .first_selected_component()];
2889
2890  for (unsigned int q = 0; q < n_q_points; ++q)
2891  {
2892  for (unsigned int d1 = 0; d1 < dim; ++d1)
2893  {
2894  for (unsigned int d2 = 0; d2 < dim; ++d2)
2895  {
2897  data.shape_values(first + d2, q) *
2899  }
2900  }
2901  }
2902
2903  for (unsigned int q = 0; q < n_q_points; ++q)
2904  {
2905  for (unsigned int d = 0; d < dim; ++d)
2906  {
2909  }
2910  }
2911  }
2912  }
2913  if (flags & update_hessians)
2914  {
2915  // Now have all shape_grads stored on the reference cell.
2916  // Must now transform to the physical cell.
2917  std::vector<Tensor<3, dim>> input(n_q_points);
2918  std::vector<Tensor<3, dim>> transformed_shape_hessians(n_q_points);
2919  for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2920  {
2921  for (unsigned int q = 0; q < n_q_points; ++q)
2922  input[q] = fe_data.shape_hessians[dof][offset + q];
2923
2924  mapping.transform(input,
2926  mapping_internal,
2927  make_array_view(transformed_shape_hessians));
2928
2929  const unsigned int first =
2930  data.shape_function_to_row_table[dof * this->n_components() +
2931  this->get_nonzero_components(dof)
2932  .first_selected_component()];
2933
2934  for (unsigned int q = 0; q < n_q_points; ++q)
2935  {
2936  for (unsigned int d1 = 0; d1 < dim; ++d1)
2937  {
2938  for (unsigned int d2 = 0; d2 < dim; ++d2)
2939  {
2940  for (unsigned int d3 = 0; d3 < dim; ++d3)
2941  {
2942  for (unsigned int d4 = 0; d4 < dim; ++d4)
2943  {
2944  transformed_shape_hessians[q][d1][d3][d4] -=
2945  (data.shape_values(first + d2, q) *
2946  mapping_data
2947  .jacobian_pushed_forward_2nd_derivatives
2948  [q][d2][d1][d3][d4]) +
2950  mapping_data
2952  [d4]) +
2954  mapping_data
2956  [d4]) +
2958  mapping_data
2960  [d1]);
2961  }
2962  }
2963  }
2964  }
2965  }
2966
2967  for (unsigned int q = 0; q < n_q_points; ++q)
2968  {
2969  for (unsigned int d = 0; d < dim; ++d)
2970  {
2971  data.shape_hessians[first + d][q] =
2972  transformed_shape_hessians[q][d];
2973  }
2974  }
2975  }
2976  }
2977 }
2978
2979
2980
2981 template <int dim, int spacedim>
2982 void
2984  const typename Triangulation<dim, dim>::cell_iterator & /*cell*/,
2985  const unsigned int /*face_no*/,
2986  const unsigned int /*sub_no*/,
2988  const Mapping<dim, dim> & /*mapping*/,
2989  const typename Mapping<dim, dim>::InternalDataBase & /*mapping_internal*/,
2990  const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2991  & /*mapping_data*/,
2992  const typename FiniteElement<dim, dim>::InternalDataBase & /*fe_internal*/,
2994  & /*data*/) const
2995 {
2996  Assert(false, ExcNotImplemented());
2997 }
2998
2999
3000
3001 template <int dim, int spacedim>
3004  const UpdateFlags flags) const
3005 {
3007
3008  if (flags & update_values)
3010
3012  out |= update_gradients | update_values |
3015
3016  if (flags & update_hessians)
3017  // Assert (false, ExcNotImplemented());
3022
3023  return out;
3024 }
3025
3026
3027
3028 template <int dim, int spacedim>
3029 std::string
3031 {
3032  // note that the FETools::get_fe_by_name function depends on the particular
3033  // format of the string this function returns, so they have to be kept in sync
3034  std::ostringstream namebuf;
3035  namebuf << "FE_NedelecSZ<" << dim << ">(" << this->degree - 1 << ")";
3036
3037  return namebuf.str();
3038 }
3039
3040
3041
3042 template <int dim, int spacedim>
3043 std::unique_ptr<FiniteElement<dim, dim>>
3045 {
3046  return std::make_unique<FE_NedelecSZ<dim, spacedim>>(*this);
3047 }
3048
3049
3050
3051 template <int dim, int spacedim>
3052 std::vector<unsigned int>
3054 {
3055  // internal function to return a vector of "dofs per object"
3056  // where the objects inside the vector refer to:
3057  // 0 = vertex
3058  // 1 = edge
3059  // 2 = face (which is a cell in 2D)
3060  // 3 = cell
3061  std::vector<unsigned int> dpo;
3062
3063  dpo.push_back(0);
3064  dpo.push_back(degree + 1);
3065  if (dim > 1)
3066  dpo.push_back(2 * degree * (degree + 1));
3067  if (dim > 2)
3068  dpo.push_back(3 * degree * degree * (degree + 1));
3069
3070  return dpo;
3071 }
3072
3073
3074
3075 template <int dim, int spacedim>
3076 unsigned int
3077 FE_NedelecSZ<dim, spacedim>::compute_num_dofs(const unsigned int degree) const
3078 {
3079  // Internal function to compute the number of DoFs
3080  // for a given dimension & polynomial order.
3081  switch (dim)
3082  {
3083  case 2:
3084  return 2 * (degree + 1) * (degree + 2);
3085
3086  case 3:
3087  return 3 * (degree + 1) * (degree + 2) * (degree + 2);
3088
3089  default:
3090  {
3091  Assert(false, ExcNotImplemented());
3092  return 0;
3093  }
3094  }
3095 }
3096
3097
3098
3099 template <int dim, int spacedim>
3100 void
3102 {
3103  // fill the 1d polynomials vector:
3104  IntegratedLegendrePolynomials =
3106 }
3107
3108
3109
3110 // explicit instantiations
3111 #include "fe_nedelec_sz.inst"
3112
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:699
std::vector< std::vector< double > > edge_lambda_grads_2d
std::vector< std::vector< Tensor< 1, dim > > > shape_values
std::vector< std::vector< std::vector< double > > > sigma_imj_values
std::vector< std::vector< double > > edge_sigma_grads
std::vector< std::vector< double > > edge_sigma_values
std::vector< std::vector< DerivativeForm< 2, dim, dim > > > shape_hessians
std::vector< std::vector< std::vector< double > > > sigma_imj_grads
std::vector< std::vector< double > > face_lambda_values
std::vector< std::vector< DerivativeForm< 1, dim, dim > > > shape_grads
std::vector< std::vector< std::vector< double > > > edge_lambda_grads_3d
std::vector< std::vector< double > > edge_lambda_values
std::vector< std::vector< double > > face_lambda_grads
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
unsigned int compute_num_dofs(const unsigned int degree) const
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
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
FE_NedelecSZ(const unsigned int order)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
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
void create_polynomials(const unsigned int degree)
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
virtual std::string get_name() const override
virtual Tensor< 2, dim > shape_grad_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual Tensor< 1, dim > shape_grad_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 2, dim > shape_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
void fill_face_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const Quadrature< dim > &quadrature, const InternalData &fedata) const
void fill_edge_values(const typename Triangulation< dim, dim >::cell_iterator &cell, const Quadrature< dim > &quadrature, const InternalData &fedata) const
MappingKind mapping_kind
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
unsigned int n_components() const
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition: fe.h:2562
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
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
Definition: point.h:111
static DataSetDescriptor face(const ReferenceCell &reference_cell, 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:1224
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
unsigned int size() const
Definition: collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
Point< 2 > first
Definition: grid_out.cc:4605
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_covariant_transformation
Covariant transformation.
@ update_default
No update.
Definition: mapping.h:93
@ mapping_covariant
Definition: mapping.h:82
@ mapping_nedelec
Definition: mapping.h:122
@ mapping_covariant_hessian
Definition: mapping.h:143
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:135
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)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)