Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_nedelec_sz.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2015 - 2023 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
18
19#include <memory>
20
22
23// Constructor:
24template <int dim, int spacedim>
26 : FiniteElement<dim, dim>(
27 FiniteElementData<dim>(get_dpo_vector(order),
28 dim,
29 order + 1,
30 FiniteElementData<dim>::Hcurl),
31 std::vector<bool>(compute_num_dofs(order), true),
32 std::vector<ComponentMask>(compute_num_dofs(order),
33 std::vector<bool>(dim, true)))
34{
35 Assert(dim >= 2, ExcImpossibleInDim(dim));
36
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:
53template <int dim, int spacedim>
54double
56 const Point<dim> & /*p*/) const
57{
59 return 0.;
60}
61
62
63
64template <int dim, int spacedim>
65double
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
78template <int dim, int spacedim>
81 const Point<dim> & /*p*/) const
82{
84 return Tensor<1, dim>();
85}
86
87
88
89template <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
102template <int dim, int spacedim>
105 const Point<dim> & /*p*/) const
106{
108 return Tensor<2, dim>();
109}
110
111
112
113template <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
126template <int dim, int spacedim>
127std::unique_ptr<typename ::FiniteElement<dim, spacedim>::InternalDataBase>
129 const UpdateFlags update_flags,
130 const Mapping<dim, spacedim> & /*mapping*/,
131 const Quadrature<dim> &quadrature,
133 spacedim>
134 & /*output_data*/) const
135{
136 std::unique_ptr<
137 typename ::FiniteElement<dim, spacedim>::InternalDataBase>
138 data_ptr = std::make_unique<InternalData>();
139 auto &data = dynamic_cast<InternalData &>(*data_ptr);
140 data.update_each = requires_update_flags(update_flags);
141
142 // Useful quantities:
143 const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
144
145 const unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
146 const unsigned int lines_per_cell = GeometryInfo<dim>::lines_per_cell;
147 const unsigned int faces_per_cell = GeometryInfo<dim>::faces_per_cell;
148
149 const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
150
151 // we assume that all quads have the same number of dofs
152 const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
153
154 const UpdateFlags flags(data.update_each);
155 const unsigned int n_q_points = quadrature.size();
156
157 // Resize the internal data storage:
158 data.sigma_imj_values.resize(
159 n_q_points,
160 std::vector<std::vector<double>>(vertices_per_cell,
161 std::vector<double>(vertices_per_cell)));
162
163 data.sigma_imj_grads.resize(vertices_per_cell,
164 std::vector<std::vector<double>>(
165 vertices_per_cell, std::vector<double>(dim)));
166
167 // Resize shape function arrays according to update flags:
168 if (flags & update_values)
169 {
170 data.shape_values.resize(this->n_dofs_per_cell(),
171 std::vector<Tensor<1, dim>>(n_q_points));
172 }
173
174 if (flags & update_gradients)
175 {
176 data.shape_grads.resize(this->n_dofs_per_cell(),
177 std::vector<DerivativeForm<1, dim, dim>>(
178 n_q_points));
179 }
180
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);
189 p_list = quadrature.get_points();
190
191
192 switch (dim)
193 {
194 case 2:
195 {
196 // Compute values of sigma & lambda and the sigma differences and
197 // lambda additions.
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.
267 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
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);
277 data.edge_sigma_grads.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)
283 data.edge_sigma_grads[m].resize(dim);
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));
295 data.edge_lambda_grads_2d.resize(lines_per_cell,
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
313 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
314 }
315 data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
316 data.edge_lambda_grads_2d[1] = {1.0, 0.0};
317 data.edge_lambda_grads_2d[2] = {0.0, -1.0};
318 data.edge_lambda_grads_2d[3] = {0.0, 1.0};
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
350 // gradients.
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 //
377 // Type 1 (gradients):
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 }
488 if (flags & update_gradients)
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);
500 data.shape_grads[dof_index1][q][0][0] =
501 4.0 * polyx[i][2] * polyy[j][0];
502 data.shape_grads[dof_index1][q][0][1] =
503 4.0 * polyx[i][1] * polyy[j][1];
504 data.shape_grads[dof_index1][q][1][0] =
505 data.shape_grads[dof_index1][q][0][1];
506 data.shape_grads[dof_index1][q][1][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);
512 data.shape_grads[dof_index2][q][0][0] =
513 data.shape_grads[dof_index1][q][0][0];
514 data.shape_grads[dof_index2][q][0][1] =
515 data.shape_grads[dof_index1][q][0][1];
516 data.shape_grads[dof_index2][q][1][0] =
517 -1.0 * data.shape_grads[dof_index1][q][1][0];
518 data.shape_grads[dof_index2][q][1][1] =
519 -1.0 * data.shape_grads[dof_index1][q][1][1];
520 }
521 // Type 3:
522 const unsigned int dof_index3_1(cell_type3_offset1 +
523 j);
524 data.shape_grads[dof_index3_1][q][0][0] = 0.0;
525 data.shape_grads[dof_index3_1][q][0][1] =
526 2.0 * polyy[j][1];
527 data.shape_grads[dof_index3_1][q][1][0] = 0.0;
528 data.shape_grads[dof_index3_1][q][1][1] = 0.0;
529
530 const unsigned int dof_index3_2(cell_type3_offset2 +
531 j);
532 data.shape_grads[dof_index3_2][q][0][0] = 0.0;
533 data.shape_grads[dof_index3_2][q][0][1] = 0.0;
534 data.shape_grads[dof_index3_2][q][1][0] =
535 2.0 * polyx[j][1];
536 data.shape_grads[dof_index3_2][q][1][1] = 0.0;
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
623 // lambda additions.
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
668 // We now want some additional information about
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.
730 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
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);
740 data.edge_sigma_grads.resize(lines_per_cell);
741 data.edge_lambda_grads_3d.resize(lines_per_cell);
742 data.edge_lambda_gradgrads_3d.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)
749 data.edge_sigma_grads[m].resize(dim);
750
751 data.edge_lambda_grads_3d[m].resize(n_q_points);
752 for (unsigned int q = 0; q < n_q_points; ++q)
753 {
754 data.edge_lambda_grads_3d[m][q].resize(dim);
755 }
756 // lambda_gradgrads are constant in a cell (no need for quad
757 // points)
758 data.edge_lambda_gradgrads_3d[m].resize(dim);
759 for (unsigned int d = 0; d < dim; ++d)
760 {
761 data.edge_lambda_gradgrads_3d[m][d].resize(dim);
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
787 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
788 }
789 // edge_lambda_grads
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 }
808 // edge_lambda gradgrads:
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 {
838 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
839 [edge_lambda_directions[m][1]] =
840 edge_lambda_sign[m];
841 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
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);
851 data.face_lambda_grads.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);
856 data.face_lambda_grads[m].resize(3);
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 }
871 // gradients are constant:
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 //
885 // Type-1 (gradients):
886 // \phi^{C_{1}}_{ijk} = grad(
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 }
1051 if (flags & update_gradients)
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
1072 data.shape_grads[dof_index1][q][0][0] =
1073 4.0 * polyx[i][2] * polyy[j][0] *
1074 polyz[k][0];
1075 data.shape_grads[dof_index1][q][0][1] =
1076 4.0 * polyx[i][1] * polyy[j][1] *
1077 polyz[k][0];
1078 data.shape_grads[dof_index1][q][0][2] =
1079 4.0 * polyx[i][1] * polyy[j][0] *
1080 polyz[k][1];
1081
1082 data.shape_grads[dof_index1][q][1][0] =
1083 data.shape_grads[dof_index1][q][0][1];
1084 data.shape_grads[dof_index1][q][1][1] =
1085 4.0 * polyx[i][0] * polyy[j][2] *
1086 polyz[k][0];
1087 data.shape_grads[dof_index1][q][1][2] =
1088 4.0 * polyx[i][0] * polyy[j][1] *
1089 polyz[k][1];
1090
1091 data.shape_grads[dof_index1][q][2][0] =
1092 data.shape_grads[dof_index1][q][0][2];
1093 data.shape_grads[dof_index1][q][2][1] =
1094 data.shape_grads[dof_index1][q][1][2];
1095 data.shape_grads[dof_index1][q][2][2] =
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 {
1107 data.shape_grads[dof_index2_1][q][0]
1108 [d] =
1109 data
1110 .shape_grads[dof_index1][q][0][d];
1111 data.shape_grads[dof_index2_1][q][1]
1112 [d] =
1113 -1.0 *
1114 data
1115 .shape_grads[dof_index1][q][1][d];
1116 data.shape_grads[dof_index2_1][q][2]
1117 [d] =
1118 data
1119 .shape_grads[dof_index1][q][2][d];
1120
1121 data.shape_grads[dof_index2_2][q][0]
1122 [d] =
1123 data
1124 .shape_grads[dof_index1][q][0][d];
1125 data.shape_grads[dof_index2_2][q][1]
1126 [d] =
1127 -1.0 *
1128 data
1129 .shape_grads[dof_index1][q][1][d];
1130 data.shape_grads[dof_index2_2][q][2]
1131 [d] =
1132 -1.0 *
1133 data
1134 .shape_grads[dof_index1][q][2][d];
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 {
1152 data.shape_grads[dof_index3_1][q][d1]
1153 [d2] = 0.0;
1154 data.shape_grads[dof_index3_2][q][d1]
1155 [d2] = 0.0;
1156 data.shape_grads[dof_index3_3][q][d1]
1157 [d2] = 0.0;
1158 }
1159 }
1160 data.shape_grads[dof_index3_1][q][0][1] =
1161 2.0 * polyy[j][1] * polyz[k][0];
1162 data.shape_grads[dof_index3_1][q][0][2] =
1163 2.0 * polyy[j][0] * polyz[k][1];
1164
1165 data.shape_grads[dof_index3_2][q][1][0] =
1166 2.0 * polyx[j][1] * polyz[k][0];
1167 data.shape_grads[dof_index3_2][q][1][2] =
1168 2.0 * polyx[j][0] * polyz[k][1];
1169
1170 data.shape_grads[dof_index3_3][q][2][0] =
1171 2.0 * polyx[j][1] * polyy[k][0];
1172 data.shape_grads[dof_index3_3][q][2][1] =
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
1470template <int dim, int spacedim>
1471void
1473 const typename Triangulation<dim, dim>::cell_iterator &cell,
1474 const Quadrature<dim> & quadrature,
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);
1583 std::vector<std::vector<double>> edge_sigma_grads(
1584 fe_data.edge_sigma_grads);
1585
1586 std::vector<std::vector<double>> edge_lambda_values(
1587 fe_data.edge_lambda_values);
1588 std::vector<std::vector<double>> edge_lambda_grads(
1589 fe_data.edge_lambda_grads_2d);
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
1601 std::transform(edge_sigma_grads[m].begin(),
1602 edge_sigma_grads[m].end(),
1603 edge_sigma_grads[m].begin(),
1604 [&](const double edge_sigma_grad) {
1605 return edge_sign[m] * edge_sigma_grad;
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] =
1640 0.5 * edge_sigma_grads[m][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] =
1651 edge_sigma_grads[m][d] *
1652 poly[poly_index][1] *
1653 edge_lambda_values[m][q] +
1654 poly[poly_index][0] *
1655 edge_lambda_grads[m][d];
1656 }
1657 }
1658 }
1659 if (flags & update_gradients)
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.
1668 fe_data.shape_grads[shift_m][q][d1][d2] =
1669 0.5 * edge_sigma_grads[m][d1] *
1670 edge_lambda_grads[m][d2];
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
1679 fe_data.shape_grads[dof_index][q][0][0] =
1680 edge_sigma_grads[m][0] *
1681 edge_sigma_grads[m][0] *
1682 edge_lambda_values[m][q] * poly[poly_index][2];
1683
1684 fe_data.shape_grads[dof_index][q][0][1] =
1685 (edge_sigma_grads[m][0] *
1686 edge_lambda_grads[m][1] +
1687 edge_sigma_grads[m][1] *
1688 edge_lambda_grads[m][0]) *
1689 poly[poly_index][1];
1690
1691 fe_data.shape_grads[dof_index][q][1][0] =
1692 fe_data.shape_grads[dof_index][q][0][1];
1693
1694 fe_data.shape_grads[dof_index][q][1][1] =
1695 edge_sigma_grads[m][1] *
1696 edge_sigma_grads[m][1] *
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] =
1728 edge_sigma_grads[m][d1] *
1729 edge_sigma_grads[m][d2] *
1730 edge_sigma_grads[m][d3] *
1731 poly[i][3] *
1732 edge_lambda_values[m][q] +
1733 poly[i][2] *
1734 (edge_sigma_grads[m][d1] *
1735 edge_sigma_grads[m][d2] *
1736 edge_lambda_grads[m][d3] +
1737 edge_sigma_grads[m][d3] *
1738 edge_sigma_grads[m][d1] *
1739 edge_lambda_grads[m][d2] +
1740 edge_sigma_grads[m][d3] *
1741 edge_sigma_grads[m][d2] *
1742 edge_lambda_grads[m][d1]);
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);
1823 std::vector<std::vector<double>> edge_sigma_grads(
1824 fe_data.edge_sigma_grads);
1825 std::vector<std::vector<std::vector<double>>> edge_lambda_grads(
1826 fe_data.edge_lambda_grads_3d);
1827 std::vector<std::vector<std::vector<double>>>
1828 edge_lambda_gradgrads_3d(fe_data.edge_lambda_gradgrads_3d);
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 });
1839 std::transform(edge_sigma_grads[m].begin(),
1840 edge_sigma_grads[m].end(),
1841 edge_sigma_grads[m].begin(),
1842 [&](const double edge_sigma_grad) {
1843 return edge_sign[m] * edge_sigma_grad;
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] =
1879 0.5 * edge_sigma_grads[m][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] =
1891 edge_sigma_grads[m][d] * poly[i][1] *
1892 edge_lambda_values[m][q] +
1893 poly[i][0] * edge_lambda_grads[m][q][d];
1894 }
1895 }
1896 }
1897 }
1898 if (flags & update_gradients)
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 {
1905 fe_data.shape_grads[shift_m][q][d1][d2] =
1906 0.5 * edge_sigma_grads[m][d1] *
1907 edge_lambda_grads[m][q][d2];
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
1922 .shape_grads[dof_index][q][d1][d2] =
1923 edge_sigma_grads[m][d1] *
1924 edge_sigma_grads[m][d2] *
1925 poly[i][2] *
1926 edge_lambda_values[m][q] +
1927 (edge_sigma_grads[m][d1] *
1928 edge_lambda_grads[m][q][d2] +
1929 edge_sigma_grads[m][d2] *
1930 edge_lambda_grads[m][q][d1]) *
1931 poly[i][1] +
1932 edge_lambda_gradgrads_3d[m][d1]
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] =
1951 0.5 * edge_sigma_grads[m][d1] *
1952 edge_lambda_gradgrads_3d[m][d3][d2];
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] =
1974 edge_sigma_grads[m][d1] *
1975 edge_sigma_grads[m][d2] *
1976 edge_sigma_grads[m][d3] *
1977 poly[i][3] *
1978 edge_lambda_values[m][q] +
1979 poly[i][2] *
1980 (edge_sigma_grads[m][d1] *
1981 edge_sigma_grads[m][d2] *
1982 edge_lambda_grads[m][q]
1983 [d3] +
1984 edge_sigma_grads[m][d3] *
1985 edge_sigma_grads[m][d1] *
1986 edge_lambda_grads[m][q]
1987 [d2] +
1988 edge_sigma_grads[m][d3] *
1989 edge_sigma_grads[m][d2] *
1990 edge_lambda_grads[m][q]
1991 [d1]) +
1992 poly[i][1] *
1993 (edge_sigma_grads[m][d1] *
1994 edge_lambda_gradgrads_3d
1995 [m][d3][d2] +
1996 edge_sigma_grads[m][d2] *
1997 edge_lambda_gradgrads_3d
1998 [m][d3][d1] +
1999 edge_sigma_grads[m][d3] *
2000 edge_lambda_gradgrads_3d
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
2022template <int dim, int spacedim>
2023void
2025 const typename Triangulation<dim, dim>::cell_iterator &cell,
2026 const Quadrature<dim> & quadrature,
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
2084 vertices_adjacent_on_face[GeometryInfo<3>::vertices_per_face][2] = {
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(
2115 m, vertices_adjacent_on_face[current_max][0])) >
2117 m, vertices_adjacent_on_face[current_max][1])))
2118 {
2119 face_orientation[m][1] =
2121 m, vertices_adjacent_on_face[current_max][0]);
2122 face_orientation[m][3] =
2124 m, vertices_adjacent_on_face[current_max][1]);
2125 }
2126 else
2127 {
2128 face_orientation[m][1] =
2130 m, vertices_adjacent_on_face[current_max][1]);
2131 face_orientation[m][3] =
2133 m, vertices_adjacent_on_face[current_max][0]);
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));
2141 std::vector<std::vector<double>> face_xi_grads(
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));
2145 std::vector<std::vector<double>> face_eta_grads(
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));
2150 std::vector<std::vector<double>> face_lambda_grads(
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 {
2166 face_xi_grads[m][d] =
2167 fe_data.sigma_imj_grads[face_orientation[m][0]]
2168 [face_orientation[m][1]][d];
2169 face_eta_grads[m][d] =
2170 fe_data.sigma_imj_grads[face_orientation[m][0]]
2171 [face_orientation[m][3]][d];
2172
2173 face_lambda_grads[m][d] = fe_data.face_lambda_grads[m][d];
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 //
2194 // Type-1 (gradients)
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] =
2259 (face_xi_grads[m][d] * polyxi[i][1] *
2260 polyeta[j][0] +
2261 face_eta_grads[m][d] * polyxi[i][0] *
2262 polyeta[j][1]) *
2263 face_lambda_values[m][q] +
2264 face_lambda_grads[m][d] * polyxi[i][0] *
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] =
2273 (face_xi_grads[m][d] * polyxi[i][1] *
2274 polyeta[j][0] -
2275 face_eta_grads[m][d] * polyxi[i][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] =
2288 face_xi_grads[m][d] * polyeta[j][0] *
2289 face_lambda_values[m][q];
2290
2291 fe_data.shape_values[dof_index3_2][q][d] =
2292 face_eta_grads[m][d] * polyxi[j][0] *
2293 face_lambda_values[m][q];
2294 }
2295 }
2296 }
2297 if (flags & update_gradients)
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
2313 .shape_grads[dof_index1][q][d1][d2] =
2314 (face_xi_grads[m][d1] *
2315 face_xi_grads[m][d2] * polyxi[i][2] *
2316 polyeta[j][0] +
2317 (face_xi_grads[m][d1] *
2318 face_eta_grads[m][d2] +
2319 face_xi_grads[m][d2] *
2320 face_eta_grads[m][d1]) *
2321 polyxi[i][1] * polyeta[j][1] +
2322 face_eta_grads[m][d1] *
2323 face_eta_grads[m][d2] *
2324 polyxi[i][0] * polyeta[j][2]) *
2325 face_lambda_values[m][q] +
2326 (face_xi_grads[m][d2] * polyxi[i][1] *
2327 polyeta[j][0] +
2328 face_eta_grads[m][d2] * polyxi[i][0] *
2329 polyeta[j][1]) *
2330 face_lambda_grads[m][d1] +
2331 (face_xi_grads[m][d1] * polyxi[i][1] *
2332 polyeta[j][0] +
2333 face_eta_grads[m][d1] * polyxi[i][0] *
2334 polyeta[j][1]) *
2335 face_lambda_grads[m][d2];
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
2346 .shape_grads[dof_index2][q][d1][d2] =
2347 (face_xi_grads[m][d1] *
2348 face_xi_grads[m][d2] * polyxi[i][2] *
2349 polyeta[j][0] +
2350 (face_xi_grads[m][d1] *
2351 face_eta_grads[m][d2] -
2352 face_xi_grads[m][d2] *
2353 face_eta_grads[m][d1]) *
2354 polyxi[i][1] * polyeta[j][1] -
2355 face_eta_grads[m][d1] *
2356 face_eta_grads[m][d2] *
2357 polyxi[i][0] * polyeta[j][2]) *
2358 face_lambda_values[m][q] +
2359 (face_xi_grads[m][d1] * polyxi[i][1] *
2360 polyeta[j][0] -
2361 face_eta_grads[m][d1] * polyxi[i][0] *
2362 polyeta[j][1]) *
2363 face_lambda_grads[m][d2];
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 {
2376 fe_data.shape_grads[dof_index3_1][q][d1][d2] =
2377 face_xi_grads[m][d1] *
2378 (face_eta_grads[m][d2] * polyeta[j][1] *
2379 face_lambda_values[m][q] +
2380 face_lambda_grads[m][d2] * polyeta[j][0]);
2381
2382 fe_data.shape_grads[dof_index3_2][q][d1][d2] =
2383 face_eta_grads[m][d1] *
2384 (face_xi_grads[m][d2] * polyxi[j][1] *
2385 face_lambda_values[m][q] +
2386 face_lambda_grads[m][d2] * polyxi[j][0]);
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] *
2412 face_xi_grads[m][d3] *
2413 (face_eta_grads[m][d1] *
2414 (polyeta[j][2] *
2415 face_eta_grads[m][d2] *
2416 face_lambda_values[m][q] +
2417 polyeta[j][1] *
2418 face_lambda_grads[m][d2]) +
2419 polyeta[j][1] *
2420 face_eta_grads[m][d2] *
2421 face_lambda_grads[m][d1]) +
2422 polyxi[i][0] *
2423 (polyeta[j][3] *
2424 face_eta_grads[m][d1] *
2425 face_eta_grads[m][d2] *
2426 face_eta_grads[m][d3] *
2427 face_lambda_values[m][q] +
2428 polyeta[j][2] *
2429 (face_eta_grads[m][d1] *
2430 face_eta_grads[m][d2] *
2431 face_lambda_grads[m][d3] +
2432 face_eta_grads[m][d3] *
2433 (face_eta_grads[m][d1] *
2434 face_lambda_grads[m]
2435 [d2] +
2436 face_eta_grads[m][d2] *
2437 face_lambda_grads
2438 [m][d1]))) +
2439 (polyxi[i][1] * polyeta[j][1] *
2440 face_eta_grads[m][d3] +
2441 polyxi[i][2] * polyeta[j][0] *
2442 face_xi_grads[m][d3]) *
2443 (face_xi_grads[m][d1] *
2444 face_lambda_grads[m][d2] +
2445 face_xi_grads[m][d2] *
2446 face_lambda_grads[m][d1]) +
2447 face_lambda_grads[m][d3] *
2448 (polyxi[i][2] * polyeta[j][0] *
2449 face_xi_grads[m][d1] *
2450 face_xi_grads[m][d2] +
2451 polyxi[i][1] * polyeta[j][1] *
2452 (face_xi_grads[m][d1] *
2453 face_eta_grads[m][d2] +
2454 face_xi_grads[m][d2] *
2455 face_eta_grads[m][d1])) +
2456 face_lambda_values[m][q] *
2457 (polyxi[i][3] * polyeta[j][0] *
2458 face_xi_grads[m][d1] *
2459 face_xi_grads[m][d2] *
2460 face_xi_grads[m][d3] +
2461 polyxi[i][1] * polyeta[j][2] *
2462 face_eta_grads[m][d3] *
2463 (face_xi_grads[m][d1] *
2464 face_eta_grads[m][d2] +
2465 face_xi_grads[m][d2] *
2466 face_eta_grads[m][d1]) +
2467 polyxi[i][2] * polyeta[j][1] *
2468 (face_xi_grads[m][d3] *
2469 face_xi_grads[m][d2] *
2470 face_eta_grads[m][d1] +
2471 face_xi_grads[m][d1] *
2472 (face_xi_grads[m][d2] *
2473 face_eta_grads[m][d3] +
2474 face_xi_grads[m][d3] *
2475 face_eta_grads[m][d2])));
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] =
2491 face_xi_grads[m][d1] *
2492 (polyxi[i][1] * polyeta[j][1] *
2493 (face_eta_grads[m][d2] *
2494 face_lambda_grads[m][d3] +
2495 face_eta_grads[m][d3] *
2496 face_lambda_grads[m][d2]) +
2497 polyxi[i][2] * polyeta[j][0] *
2498 (face_xi_grads[m][d2] *
2499 face_lambda_grads[m][d3] +
2500 face_xi_grads[m][d3] *
2501 face_lambda_grads[m][d2]) +
2502 face_lambda_values[m][q] *
2503 (face_eta_grads[m][d2] *
2504 (polyxi[i][1] *
2505 polyeta[j][2] *
2506 face_eta_grads[m][d3] +
2507 polyxi[i][2] *
2508 polyeta[j][1] *
2509 face_xi_grads[m][d3]) +
2510 face_xi_grads[m][d2] *
2511 (polyxi[i][2] *
2512 polyeta[j][1] *
2513 face_eta_grads[m][d3] +
2514 polyxi[i][3] *
2515 polyeta[j][0] *
2516 face_xi_grads[m][d3]))) -
2517 polyxi[i][0] *
2518 face_eta_grads[m][d1] *
2519 (face_eta_grads[m][d2] *
2520 (polyeta[j][3] *
2521 face_eta_grads[m][d3] *
2522 face_lambda_values[m][q] +
2523 polyeta[j][2] *
2524 face_lambda_grads[m][d3]) +
2525 polyeta[j][2] *
2526 face_eta_grads[m][d3] *
2527 face_lambda_grads[m][d2]) -
2528 face_eta_grads[m][d1] *
2529 (polyxi[i][1] *
2530 face_xi_grads[m][d3] *
2531 (polyeta[j][2] *
2532 face_eta_grads[m][d2] *
2533 face_lambda_values[m][q] +
2534 polyeta[j][1] *
2535 face_lambda_grads[m][d2]) +
2536 face_xi_grads[m][d2] *
2537 (polyxi[i][1] *
2538 (polyeta[j][2] *
2539 face_eta_grads[m][d3] *
2540 face_lambda_values[m]
2541 [q] +
2542 polyeta[j][1] *
2543 face_lambda_grads[m]
2544 [d3]) +
2545 polyxi[i][2] * polyeta[j][1] *
2546 face_xi_grads[m][d3] *
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] =
2565 face_xi_grads[m][d1] *
2566 (face_eta_grads[m][d2] *
2567 (polyeta[j][2] *
2568 face_eta_grads[m][d3] *
2569 face_lambda_values[m][q] +
2570 polyeta[j][1] *
2571 face_lambda_grads[m][d3]) +
2572 face_lambda_grads[m][d2] *
2573 polyeta[j][1] *
2574 face_eta_grads[m][d3]);
2575
2576 fe_data.shape_hessians[dof_index3_2][q]
2577 [d1][d2][d3] =
2578 face_eta_grads[m][d1] *
2579 (face_xi_grads[m][d2] *
2580 (polyxi[j][2] *
2581 face_xi_grads[m][d3] *
2582 face_lambda_values[m][q] +
2583 polyxi[j][1] *
2584 face_lambda_grads[m][d3]) +
2585 face_lambda_grads[m][d2] *
2586 polyxi[j][1] * face_xi_grads[m][d3]);
2587 }
2588 }
2589 }
2590 }
2591 }
2592 }
2593 }
2594 }
2595 }
2596}
2597
2598
2599
2600template <int dim, int spacedim>
2601void
2603 const typename Triangulation<dim, dim>::cell_iterator &cell,
2604 const CellSimilarity::Similarity /*cell_similarity*/,
2605 const Quadrature<dim> & quadrature,
2606 const Mapping<dim, dim> & mapping,
2607 const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
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,
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.
2622 fill_edge_values(cell, quadrature, fe_data);
2623 if (dim == 3 && this->degree > 1)
2624 {
2625 fill_face_values(cell, quadrature, fe_data);
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
2666 if (flags & update_gradients)
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);
2671 std::vector<Tensor<2, dim>> transformed_shape_grads(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 {
2676 input[q] = fe_data.shape_grads[dof][q];
2677 }
2678 mapping.transform(make_array_view(input),
2680 mapping_internal,
2681 make_array_view(transformed_shape_grads));
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 {
2694 transformed_shape_grads[q][d1] -=
2695 data.shape_values(first + d2, q) *
2696 mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
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 {
2705 data.shape_gradients[first + d][q] =
2706 transformed_shape_grads[q][d];
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
2748 [q][d2][d1][d3][d4]) +
2749 (data.shape_gradients[first + d1][q][d2] *
2750 mapping_data
2752 [d4]) +
2753 (data.shape_gradients[first + d2][q][d3] *
2754 mapping_data
2756 [d4]) +
2757 (data.shape_gradients[first + d2][q][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
2781template <int dim, int spacedim>
2782void
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,
2790 & mapping_data,
2791 const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
2793 &data) const
2794{
2795 AssertDimension(quadrature.size(), 1);
2796
2797 // Note for future improvement:
2798 // We don't have the full quadrature - should use QProjector to create the 2d
2799 // quadrature.
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,
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,
2818 QProjector<dim>::project_to_all_faces(this->reference_cell(),
2819 quadrature[0]),
2820 fe_data);
2821 if (dim == 3 && this->degree > 1)
2822 {
2823 fill_face_values(cell,
2825 this->reference_cell(), quadrature[0]),
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 =
2832 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
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 }
2868 if (flags & update_gradients)
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);
2873 std::vector<Tensor<2, dim>> transformed_shape_grads(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,
2883 make_array_view(transformed_shape_grads));
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 {
2896 transformed_shape_grads[q][d1] -=
2897 data.shape_values(first + d2, q) *
2898 mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
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 {
2907 data.shape_gradients[first + d][q] =
2908 transformed_shape_grads[q][d];
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
2948 [q][d2][d1][d3][d4]) +
2949 (data.shape_gradients[first + d1][q][d2] *
2950 mapping_data
2952 [d4]) +
2953 (data.shape_gradients[first + d2][q][d3] *
2954 mapping_data
2956 [d4]) +
2957 (data.shape_gradients[first + d2][q][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
2981template <int dim, int spacedim>
2982void
2984 const typename Triangulation<dim, dim>::cell_iterator & /*cell*/,
2985 const unsigned int /*face_no*/,
2986 const unsigned int /*sub_no*/,
2987 const Quadrature<dim - 1> & /*quadrature*/,
2988 const Mapping<dim, dim> & /*mapping*/,
2989 const typename Mapping<dim, dim>::InternalDataBase & /*mapping_internal*/,
2991 & /*mapping_data*/,
2992 const typename FiniteElement<dim, dim>::InternalDataBase & /*fe_internal*/,
2994 & /*data*/) const
2995{
2996 Assert(false, ExcNotImplemented());
2997}
2998
2999
3000
3001template <int dim, int spacedim>
3004 const UpdateFlags flags) const
3005{
3007
3008 if (flags & update_values)
3010
3011 if (flags & update_gradients)
3015
3016 if (flags & update_hessians)
3021
3022 return out;
3023}
3024
3025
3026
3027template <int dim, int spacedim>
3028std::string
3030{
3031 // note that the FETools::get_fe_by_name function depends on the particular
3032 // format of the string this function returns, so they have to be kept in sync
3033 std::ostringstream namebuf;
3034 namebuf << "FE_NedelecSZ<" << dim << ">(" << this->degree - 1 << ")";
3035
3036 return namebuf.str();
3037}
3038
3039
3040
3041template <int dim, int spacedim>
3042std::unique_ptr<FiniteElement<dim, dim>>
3044{
3045 return std::make_unique<FE_NedelecSZ<dim, spacedim>>(*this);
3046}
3047
3048
3049
3050template <int dim, int spacedim>
3051std::vector<unsigned int>
3053{
3054 // internal function to return a vector of "dofs per object"
3055 // where the objects inside the vector refer to:
3056 // 0 = vertex
3057 // 1 = edge
3058 // 2 = face (which is a cell in 2d)
3059 // 3 = cell
3060 std::vector<unsigned int> dpo;
3061
3062 dpo.push_back(0);
3063 dpo.push_back(degree + 1);
3064 if (dim > 1)
3065 dpo.push_back(2 * degree * (degree + 1));
3066 if (dim > 2)
3067 dpo.push_back(3 * degree * degree * (degree + 1));
3068
3069 return dpo;
3070}
3071
3072
3073
3074template <int dim, int spacedim>
3075unsigned int
3077{
3078 // Internal function to compute the number of DoFs
3079 // for a given dimension & polynomial order.
3080 switch (dim)
3081 {
3082 case 2:
3083 return 2 * (degree + 1) * (degree + 2);
3084
3085 case 3:
3086 return 3 * (degree + 1) * (degree + 2) * (degree + 2);
3087
3088 default:
3089 {
3090 Assert(false, ExcNotImplemented());
3091 return 0;
3092 }
3093 }
3094}
3095
3096
3097
3098template <int dim, int spacedim>
3099void
3101{
3102 // fill the 1d polynomials vector:
3103 IntegratedLegendrePolynomials =
3105}
3106
3107
3108
3109// explicit instantiations
3110#include "fe_nedelec_sz.inst"
3111
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:704
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< std::vector< double > > > edge_lambda_gradgrads_3d
std::vector< std::vector< double > > face_lambda_grads
unsigned int compute_num_dofs(const unsigned int degree) const
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
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
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
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)
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 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
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:2563
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
Abstract base class for mapping classes.
Definition mapping.h:317
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:112
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)
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
unsigned int size() const
Definition collection.h:265
std::vector< Tensor< 4, spacedim > > jacobian_pushed_forward_2nd_derivatives
std::vector< Tensor< 3, spacedim > > jacobian_pushed_forward_grads
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > first
Definition grid_out.cc:4615
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
UpdateFlags
@ update_jacobian_pushed_forward_2nd_derivatives
@ update_jacobian_pushed_forward_grads
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_covariant_transformation
Covariant transformation.
@ update_gradients
Shape function gradients.
@ update_default
No update.
@ mapping_covariant_gradient
Definition mapping.h:99
@ mapping_covariant
Definition mapping.h:88
@ mapping_nedelec
Definition mapping.h:128
@ mapping_covariant_hessian
Definition mapping.h:149
STL namespace.
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)