Reference documentation for deal.II version GIT relicensing-220-g40aaea9664 2024-03-28 18:00:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17#include <deal.II/fe/fe_tools.h>
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 ComponentMask(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 // Compute the face embedding.
51
52 // The implementation assumes that all faces have the same
53 // number of DoFs.
55 const unsigned int face_no = 1;
56 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
57 {
58 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
59 this->n_dofs_per_face(face_no));
60 }
61
62 FETools::compute_face_embedding_matrices<dim, double>(
63 *this, face_embeddings, 0, 0, 1.e-15 * std::exp(std::pow(order, 1.075)));
64
65 switch (dim)
66 {
67 case 1:
68 {
69 this->interface_constraints.reinit(0, 0);
70 break;
71 }
72
73 case 2:
74 {
75 this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no),
76 this->n_dofs_per_face(face_no));
77 for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
78 ++i)
79 {
80 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
81 {
82 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no);
83 ++k)
84 {
86 i * this->n_dofs_per_face(face_no) + j, k) =
87 face_embeddings[i](j, k);
88 }
89 }
90 }
91 break;
92 }
93
94 case 3:
95 {
96 this->interface_constraints.reinit(
97 4 * (this->n_dofs_per_face(face_no) - this->degree),
98 this->n_dofs_per_face(face_no));
99 unsigned int target_row = 0;
100 for (unsigned int i = 0; i < 2; ++i)
101 for (unsigned int j = this->degree; j < 2 * this->degree;
102 ++j, ++target_row)
103 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
104 this->interface_constraints(target_row, k) =
105 face_embeddings[2 * i](j, k);
106 for (unsigned int i = 0; i < 2; ++i)
107 for (unsigned int j = 3 * this->degree;
108 j < GeometryInfo<3>::lines_per_face * this->degree;
109 ++j, ++target_row)
110 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
111 this->interface_constraints(target_row, k) =
112 face_embeddings[i](j, k);
113 for (unsigned int i = 0; i < 2; ++i)
114 for (unsigned int j = 0; j < 2; ++j)
115 for (unsigned int k = i * this->degree;
116 k < (i + 1) * this->degree;
117 ++k, ++target_row)
118 for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
119 ++l)
120 this->interface_constraints(target_row, l) =
121 face_embeddings[i + 2 * j](k, l);
122 for (unsigned int i = 0; i < 2; ++i)
123 for (unsigned int j = 0; j < 2; ++j)
124 for (unsigned int k = (i + 2) * this->degree;
125 k < (i + 3) * this->degree;
126 ++k, ++target_row)
127 for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
128 ++l)
129 this->interface_constraints(target_row, l) =
130 face_embeddings[2 * i + j](k, l);
131 for (unsigned int i = 0; i < GeometryInfo<3>::max_children_per_face;
132 ++i)
133 for (unsigned int j =
134 GeometryInfo<3>::lines_per_face * this->degree;
135 j < this->n_dofs_per_face(face_no);
136 ++j, ++target_row)
137 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
138 this->interface_constraints(target_row, k) =
139 face_embeddings[i](j, k);
140 break;
141 }
142
143 default:
145 }
146}
147
148
149
150// Shape functions:
151template <int dim, int spacedim>
152double
154 const Point<dim> & /*p*/) const
155{
157 return 0.;
158}
159
160
161
162template <int dim, int spacedim>
163double
165 const unsigned int /*i*/,
166 const Point<dim> & /*p*/,
167 const unsigned int /*component*/) const
168{
169 // Not implemented yet:
171 return 0.;
172}
173
174
175
176template <int dim, int spacedim>
179 const Point<dim> & /*p*/) const
180{
182 return Tensor<1, dim>();
183}
184
185
186
187template <int dim, int spacedim>
190 const unsigned int /*i*/,
191 const Point<dim> & /*p*/,
192 const unsigned int /*component*/) const
193{
195 return Tensor<1, dim>();
196}
197
198
199
200template <int dim, int spacedim>
203 const Point<dim> & /*p*/) const
204{
206 return Tensor<2, dim>();
207}
208
209
210
211template <int dim, int spacedim>
214 const unsigned int /*i*/,
215 const Point<dim> & /*p*/,
216 const unsigned int /*component*/) const
217{
219 return Tensor<2, dim>();
220}
221
222
223
224template <int dim, int spacedim>
225std::unique_ptr<typename ::FiniteElement<dim, spacedim>::InternalDataBase>
227 const UpdateFlags update_flags,
228 const Mapping<dim, spacedim> & /*mapping*/,
229 const Quadrature<dim> &quadrature,
231 spacedim>
232 & /*output_data*/) const
233{
234 std::unique_ptr<
235 typename ::FiniteElement<dim, spacedim>::InternalDataBase>
236 data_ptr = std::make_unique<InternalData>();
237 auto &data = dynamic_cast<InternalData &>(*data_ptr);
238 data.update_each = requires_update_flags(update_flags);
239
240 // Useful quantities:
241 const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
242
243 const unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
244 const unsigned int lines_per_cell = GeometryInfo<dim>::lines_per_cell;
245 const unsigned int faces_per_cell = GeometryInfo<dim>::faces_per_cell;
246
247 const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
248
249 // we assume that all quads have the same number of dofs
250 const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
251
252 const UpdateFlags flags(data.update_each);
253 const unsigned int n_q_points = quadrature.size();
254
255 // Resize the internal data storage:
256 data.sigma_imj_values.resize(
257 n_q_points,
258 std::vector<std::vector<double>>(vertices_per_cell,
259 std::vector<double>(vertices_per_cell)));
260
261 data.sigma_imj_grads.resize(vertices_per_cell,
262 std::vector<std::vector<double>>(
263 vertices_per_cell, std::vector<double>(dim)));
264
265 // Resize shape function arrays according to update flags:
266 if (flags & update_values)
267 {
268 data.shape_values.resize(this->n_dofs_per_cell(),
269 std::vector<Tensor<1, dim>>(n_q_points));
270 }
271
272 if (flags & update_gradients)
273 {
274 data.shape_grads.resize(this->n_dofs_per_cell(),
275 std::vector<DerivativeForm<1, dim, dim>>(
276 n_q_points));
277 }
278
279 if (flags & update_hessians)
280 {
281 data.shape_hessians.resize(this->n_dofs_per_cell(),
282 std::vector<DerivativeForm<2, dim, dim>>(
283 n_q_points));
284 }
285
286 std::vector<Point<dim>> p_list(n_q_points);
287 p_list = quadrature.get_points();
288
289
290 switch (dim)
291 {
292 case 2:
293 {
294 // Compute values of sigma & lambda and the sigma differences and
295 // lambda additions.
296 std::vector<std::vector<double>> sigma(
297 n_q_points, std::vector<double>(lines_per_cell));
298 std::vector<std::vector<double>> lambda(
299 n_q_points, std::vector<double>(lines_per_cell));
300
301 for (unsigned int q = 0; q < n_q_points; ++q)
302 {
303 sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
304 sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
305 sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
306 sigma[q][3] = p_list[q][0] + p_list[q][1];
307
308 lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
309 lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
310 lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
311 lambda[q][3] = p_list[q][0] * p_list[q][1];
312 for (unsigned int i = 0; i < vertices_per_cell; ++i)
313 {
314 for (unsigned int j = 0; j < vertices_per_cell; ++j)
315 {
316 data.sigma_imj_values[q][i][j] =
317 sigma[q][i] - sigma[q][j];
318 }
319 }
320 }
321
322 // Calculate the gradient of sigma_imj_values[q][i][j] =
323 // sigma[q][i]-sigma[q][j]
324 // - this depends on the component and the direction of the
325 // corresponding edge.
326 // - the direction of the edge is determined by
327 // sigma_imj_sign[i][j].
328 // Helper arrays:
329 const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2] = {
330 {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
331 int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
332 unsigned int sigma_imj_component[vertices_per_cell]
333 [vertices_per_cell];
334
335 for (unsigned int i = 0; i < vertices_per_cell; ++i)
336 {
337 for (unsigned int j = 0; j < vertices_per_cell; ++j)
338 {
339 // sigma_imj_sign is the sign (+/-) of the coefficient of
340 // x/y/z in sigma_imj_values Due to the numbering of vertices
341 // on the reference element it is easy to find edges in the
342 // positive direction are from smaller to higher local vertex
343 // numbering.
344 sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
345 sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
346
347 // Now store the component which the sigma_i - sigma_j
348 // corresponds to:
349 sigma_imj_component[i][j] = 0;
350 for (unsigned int d = 0; d < dim; ++d)
351 {
352 int temp_imj =
353 sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
354 // Only interested in the first non-zero
355 // as if there is a second, it can not be a valid edge.
356 if (temp_imj != 0)
357 {
358 sigma_imj_component[i][j] = d;
359 break;
360 }
361 }
362 // Can now calculate the gradient, only non-zero in the
363 // component given: Note some i,j combinations will be
364 // incorrect, but only on invalid edges.
365 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
366 2.0 * sigma_imj_sign[i][j];
367 }
368 }
369
370 // Now compute the edge parameterisations for a single element
371 // with global numbering matching that of the reference element:
372
373 // Resize the edge parameterisations
374 data.edge_sigma_values.resize(lines_per_cell);
375 data.edge_sigma_grads.resize(lines_per_cell);
376 for (unsigned int m = 0; m < lines_per_cell; ++m)
377 {
378 data.edge_sigma_values[m].resize(n_q_points);
379
380 // sigma grads are constant in a cell (no need for quad points)
381 data.edge_sigma_grads[m].resize(dim);
382 }
383
384 // Fill the values for edge lambda and edge sigma:
385 const unsigned int
386 edge_sigma_direction[GeometryInfo<2>::lines_per_cell] = {1,
387 1,
388 0,
389 0};
390
391 data.edge_lambda_values.resize(lines_per_cell,
392 std::vector<double>(n_q_points));
393 data.edge_lambda_grads_2d.resize(lines_per_cell,
394 std::vector<double>(dim));
395 for (unsigned int m = 0; m < lines_per_cell; ++m)
396 {
397 // e1=max(reference vertex numbering on this edge)
398 // e2=min(reference vertex numbering on this edge)
399 // Which is guaranteed to be:
400 const unsigned int e1(
402 const unsigned int e2(
404 for (unsigned int q = 0; q < n_q_points; ++q)
405 {
406 data.edge_sigma_values[m][q] =
407 data.sigma_imj_values[q][e2][e1];
408 data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
409 }
410
411 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
412 }
413 data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
414 data.edge_lambda_grads_2d[1] = {1.0, 0.0};
415 data.edge_lambda_grads_2d[2] = {0.0, -1.0};
416 data.edge_lambda_grads_2d[3] = {0.0, 1.0};
417
418 // If the polynomial order is 0, then no more work to do:
419 if (degree < 1)
420 {
421 break;
422 }
423
424 // Otherwise, we can compute the non-cell dependent shape functions.
425 //
426 // Note: the local dof numberings follow the usual order of lines ->
427 // faces -> cells
428 // (we have no vertex-based DoFs in this element).
429 // For a given cell we have:
430 // n_line_dofs = dofs_per_line*lines_per_cell.
431 // n_face_dofs = dofs_per_face*faces_per_cell.
432 // n_cell_dofs = dofs_per_quad (2d)
433 // = dofs_per_hex (3d)
434 //
435 // i.e. For the local dof numbering:
436 // the first line dof is 0,
437 // the first face dof is n_line_dofs,
438 // the first cell dof is n_line_dofs + n_face_dofs.
439 //
440 // On a line, DoFs are ordered first by line_dof and then line_index:
441 // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
442 //
443 // and similarly for faces:
444 // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
445 //
446 // HOWEVER, we have different types of DoFs on a line/face/cell.
447 // On a line we have two types, lowest order and higher order
448 // gradients.
449 // - The numbering is such the lowest order is first, then higher
450 // order.
451 // This is simple enough as there is only 1 lowest order and
452 // degree higher orders DoFs per line.
453 //
454 // On a 2d cell, we have 3 types: Type 1/2/3:
455 // - The ordering done by type:
456 // - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
457 // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index
458 // = ij1.
459 // - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
460 // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index
461 // = degree^2 + ij2
462 // - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
463 // Numbered: ij3 = i3. i.e. cell_dof_index
464 // = 2*(degree^2) + ij3.
465 //
466 // These then fit into the local dof numbering described above:
467 // - local dof numberings are:
468 // line_dofs: local_dof = line_dof_index. 0 <= local_dof <
469 // dofs_per_line*lines_per_cell face_dofs: local_dof =
470 // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
471 // = n_lines_dof + n_face_dofs + cell_dof_index.
472 //
473 // The cell-based shape functions are:
474 //
475 // Type 1 (gradients):
476 // \phi^{C_{1}}_{ij) = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
477 //
478 // 0 <= i,j < degree.
479 //
480 // NOTE: The derivative produced by IntegratedLegendrePolynomials does
481 // not account for the
482 // (2*x-1) or (2*y-1) so we must take this into account when
483 // taking derivatives.
484 const unsigned int cell_type1_offset = n_line_dofs;
485
486 // Type 2:
487 // \phi^{C_{2}}_{ij) = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
488 // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
489 //
490 // 0 <= i,j < degree.
491 const unsigned int cell_type2_offset =
492 cell_type1_offset + degree * degree;
493
494 // Type 3 (two subtypes):
495 // \phi^{C_{3}}_{j) = L_{j+2}(2y-1) \mathbf{e}_{x}
496 //
497 // \phi^{C_{3}}_{j+degree) = L_{j+2}(2x-1) \mathbf{e}_{y}
498 //
499 // 0 <= j < degree
500 const unsigned int cell_type3_offset1 =
501 cell_type2_offset + degree * degree;
502 const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
503
505 {
506 // compute all points we must evaluate the 1d polynomials at:
507 std::vector<Point<dim>> cell_points(n_q_points);
508 for (unsigned int q = 0; q < n_q_points; ++q)
509 {
510 for (unsigned int d = 0; d < dim; ++d)
511 {
512 cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
513 }
514 }
515
516 // Loop through quad points:
517 for (unsigned int q = 0; q < n_q_points; ++q)
518 {
519 // pre-compute values & required derivatives at this quad
520 // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
521 //
522 // for each polyc[d], c=x,y, contains the d-th derivative with
523 // respect to the coordinate c.
524
525 // We only need poly values and 1st derivative for
526 // update_values, but need the 2nd derivative too for
527 // update_gradients. For update_hessians we also need the 3rd
528 // derivatives.
529 const unsigned int poly_length =
530 (flags & update_hessians) ?
531 4 :
532 ((flags & update_gradients) ? 3 : 2);
533
534 std::vector<std::vector<double>> polyx(
535 degree, std::vector<double>(poly_length));
536 std::vector<std::vector<double>> polyy(
537 degree, std::vector<double>(poly_length));
538 for (unsigned int i = 0; i < degree; ++i)
539 {
540 // Compute all required 1d polynomials and their
541 // derivatives, starting at degree 2. e.g. to access
542 // L'_{3}(2x-1) use polyx[1][1].
543 IntegratedLegendrePolynomials[i + 2].value(
544 cell_points[q][0], polyx[i]);
545 IntegratedLegendrePolynomials[i + 2].value(
546 cell_points[q][1], polyy[i]);
547 }
548 // Now use these to compute the shape functions:
549 if (flags & update_values)
550 {
551 for (unsigned int j = 0; j < degree; ++j)
552 {
553 const unsigned int shift_j(j * degree);
554 for (unsigned int i = 0; i < degree; ++i)
555 {
556 const unsigned int shift_ij(i + shift_j);
557
558 // Type 1:
559 const unsigned int dof_index1(cell_type1_offset +
560 shift_ij);
561 data.shape_values[dof_index1][q][0] =
562 2.0 * polyx[i][1] * polyy[j][0];
563 data.shape_values[dof_index1][q][1] =
564 2.0 * polyx[i][0] * polyy[j][1];
565
566 // Type 2:
567 const unsigned int dof_index2(cell_type2_offset +
568 shift_ij);
569 data.shape_values[dof_index2][q][0] =
570 data.shape_values[dof_index1][q][0];
571 data.shape_values[dof_index2][q][1] =
572 -1.0 * data.shape_values[dof_index1][q][1];
573 }
574 // Type 3:
575 const unsigned int dof_index3_1(cell_type3_offset1 +
576 j);
577 data.shape_values[dof_index3_1][q][0] = polyy[j][0];
578 data.shape_values[dof_index3_1][q][1] = 0.0;
579
580 const unsigned int dof_index3_2(cell_type3_offset2 +
581 j);
582 data.shape_values[dof_index3_2][q][0] = 0.0;
583 data.shape_values[dof_index3_2][q][1] = polyx[j][0];
584 }
585 }
586 if (flags & update_gradients)
587 {
588 for (unsigned int j = 0; j < degree; ++j)
589 {
590 const unsigned int shift_j(j * degree);
591 for (unsigned int i = 0; i < degree; ++i)
592 {
593 const unsigned int shift_ij(i + shift_j);
594
595 // Type 1:
596 const unsigned int dof_index1(cell_type1_offset +
597 shift_ij);
598 data.shape_grads[dof_index1][q][0][0] =
599 4.0 * polyx[i][2] * polyy[j][0];
600 data.shape_grads[dof_index1][q][0][1] =
601 4.0 * polyx[i][1] * polyy[j][1];
602 data.shape_grads[dof_index1][q][1][0] =
603 data.shape_grads[dof_index1][q][0][1];
604 data.shape_grads[dof_index1][q][1][1] =
605 4.0 * polyx[i][0] * polyy[j][2];
606
607 // Type 2:
608 const unsigned int dof_index2(cell_type2_offset +
609 shift_ij);
610 data.shape_grads[dof_index2][q][0][0] =
611 data.shape_grads[dof_index1][q][0][0];
612 data.shape_grads[dof_index2][q][0][1] =
613 data.shape_grads[dof_index1][q][0][1];
614 data.shape_grads[dof_index2][q][1][0] =
615 -1.0 * data.shape_grads[dof_index1][q][1][0];
616 data.shape_grads[dof_index2][q][1][1] =
617 -1.0 * data.shape_grads[dof_index1][q][1][1];
618 }
619 // Type 3:
620 const unsigned int dof_index3_1(cell_type3_offset1 +
621 j);
622 data.shape_grads[dof_index3_1][q][0][0] = 0.0;
623 data.shape_grads[dof_index3_1][q][0][1] =
624 2.0 * polyy[j][1];
625 data.shape_grads[dof_index3_1][q][1][0] = 0.0;
626 data.shape_grads[dof_index3_1][q][1][1] = 0.0;
627
628 const unsigned int dof_index3_2(cell_type3_offset2 +
629 j);
630 data.shape_grads[dof_index3_2][q][0][0] = 0.0;
631 data.shape_grads[dof_index3_2][q][0][1] = 0.0;
632 data.shape_grads[dof_index3_2][q][1][0] =
633 2.0 * polyx[j][1];
634 data.shape_grads[dof_index3_2][q][1][1] = 0.0;
635 }
636 }
637 if (flags & update_hessians)
638 {
639 for (unsigned int j = 0; j < degree; ++j)
640 {
641 const unsigned int shift_j(j * degree);
642 for (unsigned int i = 0; i < degree; ++i)
643 {
644 const unsigned int shift_ij(i + shift_j);
645
646 // Type 1:
647 const unsigned int dof_index1(cell_type1_offset +
648 shift_ij);
649 data.shape_hessians[dof_index1][q][0][0][0] =
650 8.0 * polyx[i][3] * polyy[j][0];
651 data.shape_hessians[dof_index1][q][1][0][0] =
652 8.0 * polyx[i][2] * polyy[j][1];
653
654 data.shape_hessians[dof_index1][q][0][1][0] =
655 data.shape_hessians[dof_index1][q][1][0][0];
656 data.shape_hessians[dof_index1][q][1][1][0] =
657 8.0 * polyx[i][1] * polyy[j][2];
658
659 data.shape_hessians[dof_index1][q][0][0][1] =
660 data.shape_hessians[dof_index1][q][1][0][0];
661 data.shape_hessians[dof_index1][q][1][0][1] =
662 data.shape_hessians[dof_index1][q][1][1][0];
663
664 data.shape_hessians[dof_index1][q][0][1][1] =
665 data.shape_hessians[dof_index1][q][1][1][0];
666 data.shape_hessians[dof_index1][q][1][1][1] =
667 8.0 * polyx[i][0] * polyy[j][3];
668
669
670
671 // Type 2:
672 const unsigned int dof_index2(cell_type2_offset +
673 shift_ij);
674 for (unsigned int d = 0; d < dim; ++d)
675 {
676 data.shape_hessians[dof_index2][q][0][0][d] =
677 data.shape_hessians[dof_index1][q][0][0][d];
678 data.shape_hessians[dof_index2][q][0][1][d] =
679 data.shape_hessians[dof_index1][q][0][1][d];
680 data.shape_hessians[dof_index2][q][1][0][d] =
681 -1.0 *
682 data.shape_hessians[dof_index1][q][1][0][d];
683 data.shape_hessians[dof_index2][q][1][1][d] =
684 -1.0 *
685 data.shape_hessians[dof_index1][q][1][1][d];
686 }
687 }
688 // Type 3:
689 const unsigned int dof_index3_1(cell_type3_offset1 +
690 j);
691 data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
692 data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
693 data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
694 data.shape_hessians[dof_index3_1][q][0][1][1] =
695 4.0 * polyy[j][2];
696 data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
697 data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
698 data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
699 data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
700
701 const unsigned int dof_index3_2(cell_type3_offset2 +
702 j);
703 data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
704 data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
705 data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
706 data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
707 data.shape_hessians[dof_index3_2][q][1][0][0] =
708 4.0 * polyx[j][2];
709 data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
710 data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
711 data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
712 }
713 }
714 }
715 }
716 break;
717 }
718 case 3:
719 {
720 // Compute values of sigma & lambda and the sigma differences and
721 // lambda additions.
722 std::vector<std::vector<double>> sigma(
723 n_q_points, std::vector<double>(lines_per_cell));
724 std::vector<std::vector<double>> lambda(
725 n_q_points, std::vector<double>(lines_per_cell));
726 for (unsigned int q = 0; q < n_q_points; ++q)
727 {
728 sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
729 (1 - p_list[q][2]);
730 sigma[q][1] =
731 p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
732 sigma[q][2] =
733 (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
734 sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
735 sigma[q][4] =
736 (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
737 sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
738 sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
739 sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
740
741 lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
742 (1.0 - p_list[q][2]);
743 lambda[q][1] =
744 p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
745 lambda[q][2] =
746 (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
747 lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
748 lambda[q][4] =
749 (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
750 lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
751 lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
752 lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
753
754 // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
755 // and lambda_ipj = \lambda_{i} + \lambda_{j}.
756 for (unsigned int i = 0; i < vertices_per_cell; ++i)
757 {
758 for (unsigned int j = 0; j < vertices_per_cell; ++j)
759 {
760 data.sigma_imj_values[q][i][j] =
761 sigma[q][i] - sigma[q][j];
762 }
763 }
764 }
765
766 // We now want some additional information about
767 // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
768 // calculate values & derivatives of the shape functions we need to
769 // know:
770 // - The component the sigma_imj value corresponds to - this varies
771 // with i & j.
772 // - The gradient of the sigma_imj value
773 // - this depends on the component and the direction of the
774 // corresponding edge.
775 // - the direction of the edge is determined by
776 // sigma_imj_sign[i][j].
777 //
778 // Note that not every i,j combination is a valid edge (there are only
779 // 12 valid edges in 3d), but we compute them all as it simplifies
780 // things.
781
782 // store the sign of each component x, y, z in the sigma list.
783 // can use this to fill in the sigma_imj_component data.
784 const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3] = {
785 {-1, -1, -1},
786 {1, -1, -1},
787 {-1, 1, -1},
788 {1, 1, -1},
789 {-1, -1, 1},
790 {1, -1, 1},
791 {-1, 1, 1},
792 {1, 1, 1}};
793
794 int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
795 unsigned int sigma_imj_component[vertices_per_cell]
796 [vertices_per_cell];
797
798 for (unsigned int i = 0; i < vertices_per_cell; ++i)
799 {
800 for (unsigned int j = 0; j < vertices_per_cell; ++j)
801 {
802 // sigma_imj_sign is the sign (+/-) of the coefficient of
803 // x/y/z in sigma_imj. Due to the numbering of vertices on the
804 // reference element this is easy to work out because edges in
805 // the positive direction go from smaller to higher local
806 // vertex numbering.
807 sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
808 sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
809
810 // Now store the component which the sigma_i - sigma_j
811 // corresponds to:
812 sigma_imj_component[i][j] = 0;
813 for (unsigned int d = 0; d < dim; ++d)
814 {
815 int temp_imj =
816 sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
817 // Only interested in the first non-zero
818 // as if there is a second, it will not be a valid edge.
819 if (temp_imj != 0)
820 {
821 sigma_imj_component[i][j] = d;
822 break;
823 }
824 }
825 // Can now calculate the gradient, only non-zero in the
826 // component given: Note some i,j combinations will be
827 // incorrect, but only on invalid edges.
828 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
829 2.0 * sigma_imj_sign[i][j];
830 }
831 }
832 // Now compute the edge parameterisations for a single element
833 // with global numbering matching that of the reference element:
834
835 // resize the edge parameterisations
836 data.edge_sigma_values.resize(lines_per_cell);
837 data.edge_lambda_values.resize(lines_per_cell);
838 data.edge_sigma_grads.resize(lines_per_cell);
839 data.edge_lambda_grads_3d.resize(lines_per_cell);
840 data.edge_lambda_gradgrads_3d.resize(lines_per_cell);
841 for (unsigned int m = 0; m < lines_per_cell; ++m)
842 {
843 data.edge_sigma_values[m].resize(n_q_points);
844 data.edge_lambda_values[m].resize(n_q_points);
845
846 // sigma grads are constant in a cell (no need for quad points)
847 data.edge_sigma_grads[m].resize(dim);
848
849 data.edge_lambda_grads_3d[m].resize(n_q_points);
850 for (unsigned int q = 0; q < n_q_points; ++q)
851 {
852 data.edge_lambda_grads_3d[m][q].resize(dim);
853 }
854 // lambda_gradgrads are constant in a cell (no need for quad
855 // points)
856 data.edge_lambda_gradgrads_3d[m].resize(dim);
857 for (unsigned int d = 0; d < dim; ++d)
858 {
859 data.edge_lambda_gradgrads_3d[m][d].resize(dim);
860 }
861 }
862
863 // Fill the values:
864 const unsigned int
865 edge_sigma_direction[GeometryInfo<3>::lines_per_cell] = {
866 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
867
868 for (unsigned int m = 0; m < lines_per_cell; ++m)
869 {
870 // e1=max(reference vertex numbering on this edge)
871 // e2=min(reference vertex numbering on this edge)
872 // Which is guaranteed to be:
873 const unsigned int e1(
875 const unsigned int e2(
877
878 for (unsigned int q = 0; q < n_q_points; ++q)
879 {
880 data.edge_sigma_values[m][q] =
881 data.sigma_imj_values[q][e2][e1];
882 data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
883 }
884
885 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
886 }
887 // edge_lambda_grads
888 for (unsigned int q = 0; q < n_q_points; ++q)
889 {
890 double x(p_list[q][0]);
891 double y(p_list[q][1]);
892 double z(p_list[q][2]);
893 data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
894 data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
895 data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
896 data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
897 data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
898 data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
899 data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
900 data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
901 data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
902 data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
903 data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
904 data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
905 }
906 // edge_lambda gradgrads:
907 const int edge_lambda_sign[GeometryInfo<3>::lines_per_cell] = {
908 1,
909 -1,
910 1,
911 -1,
912 -1,
913 1,
914 -1,
915 1,
916 1,
917 -1,
918 -1,
919 1}; // sign of the 2nd derivative for each edge.
920 const unsigned int
921 edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2] = {
922 {0, 2},
923 {0, 2},
924 {1, 2},
925 {1, 2},
926 {0, 2},
927 {0, 2},
928 {1, 2},
929 {1, 2},
930 {0, 1},
931 {0, 1},
932 {0, 1},
933 {0, 1}}; // component which edge_lambda[m] depends on.
934 for (unsigned int m = 0; m < lines_per_cell; ++m)
935 {
936 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
937 [edge_lambda_directions[m][1]] =
938 edge_lambda_sign[m];
939 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
940 [edge_lambda_directions[m][0]] =
941 edge_lambda_sign[m];
942 }
943 // Precomputation for higher order shape functions,
944 // and the face parameterisation.
945 if (degree > 0)
946 {
947 // resize required data:
948 data.face_lambda_values.resize(faces_per_cell);
949 data.face_lambda_grads.resize(faces_per_cell);
950 // for face-based shape functions:
951 for (unsigned int m = 0; m < faces_per_cell; ++m)
952 {
953 data.face_lambda_values[m].resize(n_q_points);
954 data.face_lambda_grads[m].resize(3);
955 }
956 // Fill in the values (these don't change between cells).
957 for (unsigned int q = 0; q < n_q_points; ++q)
958 {
959 double x(p_list[q][0]);
960 double y(p_list[q][1]);
961 double z(p_list[q][2]);
962 data.face_lambda_values[0][q] = 1.0 - x;
963 data.face_lambda_values[1][q] = x;
964 data.face_lambda_values[2][q] = 1.0 - y;
965 data.face_lambda_values[3][q] = y;
966 data.face_lambda_values[4][q] = 1.0 - z;
967 data.face_lambda_values[5][q] = z;
968 }
969 // gradients are constant:
970 data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
971 data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
972 data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
973 data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
974 data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
975 data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
976
977 // for cell-based shape functions:
978 // these don't depend on the cell, so can precompute all here:
980 {
981 // Cell-based shape functions:
982 //
983 // Type-1 (gradients):
984 // \phi^{C_{1}}_{ijk} = grad(
985 // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
986 //
987 // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
988 const unsigned int cell_type1_offset(n_line_dofs +
989 n_face_dofs);
990 // Type-2:
991 //
992 // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
993 // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
994 // -1)\phi^{C_{1}}_{ijk}
995 //
996 // 0 <= i,j,k < degree. (subtypes in groups of
997 // degree*degree*degree)
998 //
999 // here we order so that all of subtype 1 comes first, then
1000 // subtype 2.
1001 const unsigned int cell_type2_offset1(
1002 cell_type1_offset + degree * degree * degree);
1003 const unsigned int cell_type2_offset2(
1004 cell_type2_offset1 + degree * degree * degree);
1005 // Type-3
1006 // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
1007 // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
1008 // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
1009 //
1010 // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
1011 //
1012 // again we order so we compute all of subtype 1 first, then
1013 // subtype 2, etc.
1014 const unsigned int cell_type3_offset1(
1015 cell_type2_offset2 + degree * degree * degree);
1016 const unsigned int cell_type3_offset2(cell_type3_offset1 +
1017 degree * degree);
1018 const unsigned int cell_type3_offset3(cell_type3_offset2 +
1019 degree * degree);
1020
1021 // compute all points we must evaluate the 1d polynomials at:
1022 std::vector<Point<dim>> cell_points(n_q_points);
1023 for (unsigned int q = 0; q < n_q_points; ++q)
1024 {
1025 for (unsigned int d = 0; d < dim; ++d)
1026 {
1027 cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
1028 }
1029 }
1030
1031 // We only need poly values and 1st derivative for
1032 // update_values, but need the 2nd derivative too for
1033 // update_gradients. For update_hessians we also need 3rd
1034 // derivative.
1035 const unsigned int poly_length =
1036 (flags & update_hessians) ?
1037 4 :
1038 ((flags & update_gradients) ? 3 : 2);
1039
1040 // Loop through quad points:
1041 for (unsigned int q = 0; q < n_q_points; ++q)
1042 {
1043 // pre-compute values & required derivatives at this quad
1044 // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
1045 // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
1046 //
1047 // for each polyc[d], c=x,y,z, contains the d-th
1048 // derivative with respect to the coordinate c.
1049 std::vector<std::vector<double>> polyx(
1050 degree, std::vector<double>(poly_length));
1051 std::vector<std::vector<double>> polyy(
1052 degree, std::vector<double>(poly_length));
1053 std::vector<std::vector<double>> polyz(
1054 degree, std::vector<double>(poly_length));
1055 for (unsigned int i = 0; i < degree; ++i)
1056 {
1057 // compute all required 1d polynomials for i
1058 IntegratedLegendrePolynomials[i + 2].value(
1059 cell_points[q][0], polyx[i]);
1060 IntegratedLegendrePolynomials[i + 2].value(
1061 cell_points[q][1], polyy[i]);
1062 IntegratedLegendrePolynomials[i + 2].value(
1063 cell_points[q][2], polyz[i]);
1064 }
1065 // Now use these to compute the shape functions:
1066 if (flags & update_values)
1067 {
1068 for (unsigned int k = 0; k < degree; ++k)
1069 {
1070 const unsigned int shift_k(k * degree * degree);
1071 const unsigned int shift_j(
1072 k * degree); // Used below when subbing k for j
1073 // (type 3)
1074 for (unsigned int j = 0; j < degree; ++j)
1075 {
1076 const unsigned int shift_jk(j * degree +
1077 shift_k);
1078 for (unsigned int i = 0; i < degree; ++i)
1079 {
1080 const unsigned int shift_ijk(shift_jk +
1081 i);
1082
1083 // Type 1:
1084 const unsigned int dof_index1(
1085 cell_type1_offset + shift_ijk);
1086
1087 data.shape_values[dof_index1][q][0] =
1088 2.0 * polyx[i][1] * polyy[j][0] *
1089 polyz[k][0];
1090 data.shape_values[dof_index1][q][1] =
1091 2.0 * polyx[i][0] * polyy[j][1] *
1092 polyz[k][0];
1093 data.shape_values[dof_index1][q][2] =
1094 2.0 * polyx[i][0] * polyy[j][0] *
1095 polyz[k][1];
1096
1097 // Type 2:
1098 const unsigned int dof_index2_1(
1099 cell_type2_offset1 + shift_ijk);
1100 const unsigned int dof_index2_2(
1101 cell_type2_offset2 + shift_ijk);
1102
1103 data.shape_values[dof_index2_1][q][0] =
1104 data.shape_values[dof_index1][q][0];
1105 data.shape_values[dof_index2_1][q][1] =
1106 -1.0 *
1107 data.shape_values[dof_index1][q][1];
1108 data.shape_values[dof_index2_1][q][2] =
1109 data.shape_values[dof_index1][q][2];
1110
1111 data.shape_values[dof_index2_2][q][0] =
1112 data.shape_values[dof_index1][q][0];
1113 data.shape_values[dof_index2_2][q][1] =
1114 -1.0 *
1115 data.shape_values[dof_index1][q][1];
1116 data.shape_values[dof_index2_2][q][2] =
1117 -1.0 *
1118 data.shape_values[dof_index1][q][2];
1119 }
1120 // Type 3: (note we re-use k and j for
1121 // convenience):
1122 const unsigned int shift_ij(
1123 j + shift_j); // here we've subbed j for i,
1124 // k for j.
1125 const unsigned int dof_index3_1(
1126 cell_type3_offset1 + shift_ij);
1127 const unsigned int dof_index3_2(
1128 cell_type3_offset2 + shift_ij);
1129 const unsigned int dof_index3_3(
1130 cell_type3_offset3 + shift_ij);
1131
1132 data.shape_values[dof_index3_1][q][0] =
1133 polyy[j][0] * polyz[k][0];
1134 data.shape_values[dof_index3_1][q][1] = 0.0;
1135 data.shape_values[dof_index3_1][q][2] = 0.0;
1136
1137 data.shape_values[dof_index3_2][q][0] = 0.0;
1138 data.shape_values[dof_index3_2][q][1] =
1139 polyx[j][0] * polyz[k][0];
1140 data.shape_values[dof_index3_2][q][2] = 0.0;
1141
1142 data.shape_values[dof_index3_3][q][0] = 0.0;
1143 data.shape_values[dof_index3_3][q][1] = 0.0;
1144 data.shape_values[dof_index3_3][q][2] =
1145 polyx[j][0] * polyy[k][0];
1146 }
1147 }
1148 }
1149 if (flags & update_gradients)
1150 {
1151 for (unsigned int k = 0; k < degree; ++k)
1152 {
1153 const unsigned int shift_k(k * degree * degree);
1154 const unsigned int shift_j(
1155 k * degree); // Used below when subbing k for j
1156 // (type 3)
1157 for (unsigned int j = 0; j < degree; ++j)
1158 {
1159 const unsigned int shift_jk(j * degree +
1160 shift_k);
1161 for (unsigned int i = 0; i < degree; ++i)
1162 {
1163 const unsigned int shift_ijk(shift_jk +
1164 i);
1165
1166 // Type 1:
1167 const unsigned int dof_index1(
1168 cell_type1_offset + shift_ijk);
1169
1170 data.shape_grads[dof_index1][q][0][0] =
1171 4.0 * polyx[i][2] * polyy[j][0] *
1172 polyz[k][0];
1173 data.shape_grads[dof_index1][q][0][1] =
1174 4.0 * polyx[i][1] * polyy[j][1] *
1175 polyz[k][0];
1176 data.shape_grads[dof_index1][q][0][2] =
1177 4.0 * polyx[i][1] * polyy[j][0] *
1178 polyz[k][1];
1179
1180 data.shape_grads[dof_index1][q][1][0] =
1181 data.shape_grads[dof_index1][q][0][1];
1182 data.shape_grads[dof_index1][q][1][1] =
1183 4.0 * polyx[i][0] * polyy[j][2] *
1184 polyz[k][0];
1185 data.shape_grads[dof_index1][q][1][2] =
1186 4.0 * polyx[i][0] * polyy[j][1] *
1187 polyz[k][1];
1188
1189 data.shape_grads[dof_index1][q][2][0] =
1190 data.shape_grads[dof_index1][q][0][2];
1191 data.shape_grads[dof_index1][q][2][1] =
1192 data.shape_grads[dof_index1][q][1][2];
1193 data.shape_grads[dof_index1][q][2][2] =
1194 4.0 * polyx[i][0] * polyy[j][0] *
1195 polyz[k][2];
1196
1197 // Type 2:
1198 const unsigned int dof_index2_1(
1199 cell_type2_offset1 + shift_ijk);
1200 const unsigned int dof_index2_2(
1201 cell_type2_offset2 + shift_ijk);
1202
1203 for (unsigned int d = 0; d < dim; ++d)
1204 {
1205 data.shape_grads[dof_index2_1][q][0]
1206 [d] =
1207 data
1208 .shape_grads[dof_index1][q][0][d];
1209 data.shape_grads[dof_index2_1][q][1]
1210 [d] =
1211 -1.0 *
1212 data
1213 .shape_grads[dof_index1][q][1][d];
1214 data.shape_grads[dof_index2_1][q][2]
1215 [d] =
1216 data
1217 .shape_grads[dof_index1][q][2][d];
1218
1219 data.shape_grads[dof_index2_2][q][0]
1220 [d] =
1221 data
1222 .shape_grads[dof_index1][q][0][d];
1223 data.shape_grads[dof_index2_2][q][1]
1224 [d] =
1225 -1.0 *
1226 data
1227 .shape_grads[dof_index1][q][1][d];
1228 data.shape_grads[dof_index2_2][q][2]
1229 [d] =
1230 -1.0 *
1231 data
1232 .shape_grads[dof_index1][q][2][d];
1233 }
1234 }
1235 // Type 3: (note we re-use k and j for
1236 // convenience):
1237 const unsigned int shift_ij(
1238 j + shift_j); // here we've subbed j for i,
1239 // k for j.
1240 const unsigned int dof_index3_1(
1241 cell_type3_offset1 + shift_ij);
1242 const unsigned int dof_index3_2(
1243 cell_type3_offset2 + shift_ij);
1244 const unsigned int dof_index3_3(
1245 cell_type3_offset3 + shift_ij);
1246 for (unsigned int d1 = 0; d1 < dim; ++d1)
1247 {
1248 for (unsigned int d2 = 0; d2 < dim; ++d2)
1249 {
1250 data.shape_grads[dof_index3_1][q][d1]
1251 [d2] = 0.0;
1252 data.shape_grads[dof_index3_2][q][d1]
1253 [d2] = 0.0;
1254 data.shape_grads[dof_index3_3][q][d1]
1255 [d2] = 0.0;
1256 }
1257 }
1258 data.shape_grads[dof_index3_1][q][0][1] =
1259 2.0 * polyy[j][1] * polyz[k][0];
1260 data.shape_grads[dof_index3_1][q][0][2] =
1261 2.0 * polyy[j][0] * polyz[k][1];
1262
1263 data.shape_grads[dof_index3_2][q][1][0] =
1264 2.0 * polyx[j][1] * polyz[k][0];
1265 data.shape_grads[dof_index3_2][q][1][2] =
1266 2.0 * polyx[j][0] * polyz[k][1];
1267
1268 data.shape_grads[dof_index3_3][q][2][0] =
1269 2.0 * polyx[j][1] * polyy[k][0];
1270 data.shape_grads[dof_index3_3][q][2][1] =
1271 2.0 * polyx[j][0] * polyy[k][1];
1272 }
1273 }
1274 }
1275 if (flags & update_hessians)
1276 {
1277 for (unsigned int k = 0; k < degree; ++k)
1278 {
1279 const unsigned int shift_k(k * degree * degree);
1280 const unsigned int shift_j(
1281 k * degree); // Used below when subbing k for j
1282 // type 3
1283
1284 for (unsigned int j = 0; j < degree; ++j)
1285 {
1286 const unsigned int shift_jk(j * degree +
1287 shift_k);
1288 for (unsigned int i = 0; i < degree; ++i)
1289 {
1290 const unsigned int shift_ijk(shift_jk +
1291 i);
1292
1293 // Type 1:
1294 const unsigned int dof_index1(
1295 cell_type1_offset + shift_ijk);
1296
1297 data.shape_hessians[dof_index1][q][0][0]
1298 [0] =
1299 8.0 * polyx[i][3] * polyy[j][0] *
1300 polyz[k][0];
1301 data.shape_hessians[dof_index1][q][1][0]
1302 [0] =
1303 8.0 * polyx[i][2] * polyy[j][1] *
1304 polyz[k][0];
1305 data.shape_hessians[dof_index1][q][2][0]
1306 [0] =
1307 8.0 * polyx[i][2] * polyy[j][0] *
1308 polyz[k][1];
1309
1310 data.shape_hessians[dof_index1][q][0][1]
1311 [0] =
1312 data.shape_hessians[dof_index1][q][1][0]
1313 [0];
1314 data.shape_hessians[dof_index1][q][1][1]
1315 [0] =
1316 8.0 * polyx[i][1] * polyy[j][2] *
1317 polyz[k][0];
1318 data.shape_hessians[dof_index1][q][2][1]
1319 [0] =
1320 8.0 * polyx[i][1] * polyy[j][1] *
1321 polyz[k][1];
1322
1323 data.shape_hessians[dof_index1][q][0][2]
1324 [0] =
1325 data.shape_hessians[dof_index1][q][2][0]
1326 [0];
1327 data.shape_hessians[dof_index1][q][1][2]
1328 [0] =
1329 data.shape_hessians[dof_index1][q][2][1]
1330 [0];
1331 data.shape_hessians[dof_index1][q][2][2]
1332 [0] =
1333 8.0 * polyx[i][1] * polyy[j][0] *
1334 polyz[k][2];
1335
1336
1337 data.shape_hessians[dof_index1][q][0][0]
1338 [1] =
1339 data.shape_hessians[dof_index1][q][1][0]
1340 [0];
1341 data.shape_hessians[dof_index1][q][1][0]
1342 [1] =
1343 data.shape_hessians[dof_index1][q][1][1]
1344 [0];
1345 data.shape_hessians[dof_index1][q][2][0]
1346 [1] =
1347 data.shape_hessians[dof_index1][q][2][1]
1348 [0];
1349
1350 data.shape_hessians[dof_index1][q][0][1]
1351 [1] =
1352 data.shape_hessians[dof_index1][q][1][1]
1353 [0];
1354 data.shape_hessians[dof_index1][q][1][1]
1355 [1] =
1356 8.0 * polyx[i][0] * polyy[j][3] *
1357 polyz[k][0];
1358 data.shape_hessians[dof_index1][q][2][1]
1359 [1] =
1360 8.0 * polyx[i][0] * polyy[j][2] *
1361 polyz[k][1];
1362
1363 data.shape_hessians[dof_index1][q][0][2]
1364 [1] =
1365 data.shape_hessians[dof_index1][q][2][1]
1366 [0];
1367 data.shape_hessians[dof_index1][q][1][2]
1368 [1] =
1369 data.shape_hessians[dof_index1][q][2][1]
1370 [1];
1371 data.shape_hessians[dof_index1][q][2][2]
1372 [1] =
1373 8.0 * polyx[i][0] * polyy[j][1] *
1374 polyz[k][2];
1375
1376
1377 data.shape_hessians[dof_index1][q][0][0]
1378 [2] =
1379 data.shape_hessians[dof_index1][q][2][0]
1380 [0];
1381 data.shape_hessians[dof_index1][q][1][0]
1382 [2] =
1383 data.shape_hessians[dof_index1][q][2][1]
1384 [0];
1385 data.shape_hessians[dof_index1][q][2][0]
1386 [2] =
1387 data.shape_hessians[dof_index1][q][2][2]
1388 [0];
1389
1390 data.shape_hessians[dof_index1][q][0][1]
1391 [2] =
1392 data.shape_hessians[dof_index1][q][2][1]
1393 [0];
1394 data.shape_hessians[dof_index1][q][1][1]
1395 [2] =
1396 data.shape_hessians[dof_index1][q][2][1]
1397 [1];
1398 data.shape_hessians[dof_index1][q][2][1]
1399 [2] =
1400 data.shape_hessians[dof_index1][q][2][2]
1401 [1];
1402
1403 data.shape_hessians[dof_index1][q][0][2]
1404 [2] =
1405 data.shape_hessians[dof_index1][q][2][2]
1406 [0];
1407 data.shape_hessians[dof_index1][q][1][2]
1408 [2] =
1409 data.shape_hessians[dof_index1][q][2][2]
1410 [1];
1411 data.shape_hessians[dof_index1][q][2][2]
1412 [2] =
1413 8.0 * polyx[i][0] * polyy[j][0] *
1414 polyz[k][3];
1415
1416
1417 // Type 2:
1418 const unsigned int dof_index2_1(
1419 cell_type2_offset1 + shift_ijk);
1420 const unsigned int dof_index2_2(
1421 cell_type2_offset2 + shift_ijk);
1422
1423 for (unsigned int d1 = 0; d1 < dim; ++d1)
1424 {
1425 for (unsigned int d2 = 0; d2 < dim;
1426 ++d2)
1427 {
1428 data
1429 .shape_hessians[dof_index2_1][q]
1430 [0][d1][d2] =
1431 data
1432 .shape_hessians[dof_index1][q]
1433 [0][d1][d2];
1434 data
1435 .shape_hessians[dof_index2_1][q]
1436 [1][d1][d2] =
1437 -1.0 *
1438 data
1439 .shape_hessians[dof_index1][q]
1440 [1][d1][d2];
1441 data
1442 .shape_hessians[dof_index2_1][q]
1443 [2][d1][d2] =
1444 data
1445 .shape_hessians[dof_index1][q]
1446 [2][d1][d2];
1447
1448 data
1449 .shape_hessians[dof_index2_2][q]
1450 [0][d1][d2] =
1451 data
1452 .shape_hessians[dof_index1][q]
1453 [0][d1][d2];
1454 data
1455 .shape_hessians[dof_index2_2][q]
1456 [1][d1][d2] =
1457 -1.0 *
1458 data
1459 .shape_hessians[dof_index1][q]
1460 [1][d1][d2];
1461 data
1462 .shape_hessians[dof_index2_2][q]
1463 [2][d1][d2] =
1464 -1.0 *
1465 data
1466 .shape_hessians[dof_index1][q]
1467 [2][d1][d2];
1468 }
1469 }
1470 }
1471 // Type 3: (note we re-use k and j for
1472 // convenience):
1473 const unsigned int shift_ij(
1474 j + shift_j); // here we've subbed j for i,
1475 // k for j.
1476 const unsigned int dof_index3_1(
1477 cell_type3_offset1 + shift_ij);
1478 const unsigned int dof_index3_2(
1479 cell_type3_offset2 + shift_ij);
1480 const unsigned int dof_index3_3(
1481 cell_type3_offset3 + shift_ij);
1482 for (unsigned int d1 = 0; d1 < dim; ++d1)
1483 {
1484 for (unsigned int d2 = 0; d2 < dim; ++d2)
1485 {
1486 for (unsigned int d3 = 0; d3 < dim;
1487 ++d3)
1488 {
1489 data
1490 .shape_hessians[dof_index3_1][q]
1491 [d1][d2][d3] =
1492 0.0;
1493 data
1494 .shape_hessians[dof_index3_2][q]
1495 [d1][d2][d3] =
1496 0.0;
1497 data
1498 .shape_hessians[dof_index3_3][q]
1499 [d1][d2][d3] =
1500 0.0;
1501 }
1502 }
1503 }
1504 data
1505 .shape_hessians[dof_index3_1][q][0][1][1] =
1506 4.0 * polyy[j][2] * polyz[k][0];
1507 data
1508 .shape_hessians[dof_index3_1][q][0][1][2] =
1509 4.0 * polyy[j][1] * polyz[k][1];
1510
1511 data
1512 .shape_hessians[dof_index3_1][q][0][2][1] =
1513 data
1514 .shape_hessians[dof_index3_1][q][0][1][2];
1515 data
1516 .shape_hessians[dof_index3_1][q][0][2][2] =
1517 4.0 * polyy[j][0] * polyz[k][2];
1518
1519
1520 data
1521 .shape_hessians[dof_index3_2][q][1][0][0] =
1522 4.0 * polyx[j][2] * polyz[k][0];
1523 data
1524 .shape_hessians[dof_index3_2][q][1][0][2] =
1525 4.0 * polyx[j][1] * polyz[k][1];
1526
1527 data
1528 .shape_hessians[dof_index3_2][q][1][2][0] =
1529 data
1530 .shape_hessians[dof_index3_2][q][1][0][2];
1531 data
1532 .shape_hessians[dof_index3_2][q][1][2][2] =
1533 4.0 * polyx[j][0] * polyz[k][2];
1534
1535
1536 data
1537 .shape_hessians[dof_index3_3][q][2][0][0] =
1538 4.0 * polyx[j][2] * polyy[k][0];
1539 data
1540 .shape_hessians[dof_index3_3][q][2][0][1] =
1541 4.0 * polyx[j][1] * polyy[k][1];
1542
1543 data
1544 .shape_hessians[dof_index3_3][q][2][1][0] =
1545 data
1546 .shape_hessians[dof_index3_3][q][2][0][1];
1547 data
1548 .shape_hessians[dof_index3_3][q][2][1][1] =
1549 4.0 * polyx[j][0] * polyy[k][2];
1550 }
1551 }
1552 }
1553 }
1554 }
1555 }
1556 break;
1557 }
1558 default:
1559 {
1561 }
1562 }
1563 return data_ptr;
1564}
1565
1566
1567
1568template <int dim, int spacedim>
1569void
1571 const typename Triangulation<dim, dim>::cell_iterator &cell,
1572 const Quadrature<dim> &quadrature,
1573 const InternalData &fe_data) const
1574{
1575 // This function handles the cell-dependent construction of the EDGE-based
1576 // shape functions.
1577 //
1578 // Note it will handle both 2d and 3d, in 2d, the edges are faces, but we
1579 // handle them here.
1580 //
1581 // It will fill in the missing parts of fe_data which were not possible to
1582 // fill in the get_data routine, with respect to the edge-based shape
1583 // functions.
1584 //
1585 // It should be called by the fill_fe_*_values routines in order to complete
1586 // the basis set at quadrature points on the current cell for each edge.
1587
1588 const UpdateFlags flags(fe_data.update_each);
1589 const unsigned int n_q_points = quadrature.size();
1590
1591 Assert(!(flags & update_values) ||
1592 fe_data.shape_values.size() == this->n_dofs_per_cell(),
1593 ExcDimensionMismatch(fe_data.shape_values.size(),
1594 this->n_dofs_per_cell()));
1595 Assert(!(flags & update_values) ||
1596 fe_data.shape_values[0].size() == n_q_points,
1597 ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
1598
1599 // Useful constants:
1600 const unsigned int degree(
1601 this->degree -
1602 1); // Note: constructor takes input degree + 1, so need to knock 1 off.
1603
1604 // Useful geometry info:
1605 const unsigned int vertices_per_line(2);
1606 const unsigned int lines_per_cell(GeometryInfo<dim>::lines_per_cell);
1607
1608 // Calculate edge orderings:
1609 std::vector<std::vector<unsigned int>> edge_order(
1610 lines_per_cell, std::vector<unsigned int>(vertices_per_line));
1611
1612
1613 switch (dim)
1614 {
1615 case 2:
1616 {
1618 {
1619 // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1},
1620 // e^{m}_{2}] e1 = higher global numbering of the two local
1621 // vertices e2 = lower global numbering of the two local vertices
1622 std::vector<int> edge_sign(lines_per_cell);
1623 for (unsigned int m = 0; m < lines_per_cell; ++m)
1624 {
1625 unsigned int v0_loc =
1627 unsigned int v1_loc =
1629 unsigned int v0_glob = cell->vertex_index(v0_loc);
1630 unsigned int v1_glob = cell->vertex_index(v1_loc);
1631
1632 // Check for hanging edges on the current face. If we
1633 // encounter a hanging edge, we use the vertex indices
1634 // from the parent.
1635 if (cell->face(m)->at_boundary() == false)
1636 if (cell->neighbor_is_coarser(m))
1637 {
1638 v0_glob = cell->parent()->vertex_index(v0_loc);
1639 v1_glob = cell->parent()->vertex_index(v1_loc);
1640 }
1641
1642 if (v0_glob > v1_glob)
1643 {
1644 // Opposite to global numbering on our reference element
1645 edge_sign[m] = -1.0;
1646 }
1647 else
1648 {
1649 // Aligns with global numbering on our reference element.
1650 edge_sign[m] = 1.0;
1651 }
1652 }
1653
1654 // Define \sigma_{m} = sigma_{e^{m}_{2}} - sigma_{e^{m}_{1}}
1655 // \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
1656 //
1657 // To help things, in fe_data, we have precomputed (sigma_{i} -
1658 // sigma_{j}) and (lambda_{i} + lambda_{j}) for 0<= i,j <
1659 // lines_per_cell.
1660 //
1661 // There are two types:
1662 // - lower order (1 per edge, m):
1663 // \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
1664 //
1665 // - higher order (degree per edge, m):
1666 // \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
1667 //
1668 // NOTE: sigma_{m} and lambda_{m} are either a function of x OR
1669 // y
1670 // and if sigma is of x, then lambda is of y, and vice
1671 // versa. This means that grad(\sigma) requires
1672 // multiplication by d(sigma)/dx_{i} for the i^th comp of
1673 // grad(sigma) and similarly when taking derivatives of
1674 // lambda.
1675 //
1676 // First handle the lowest order edges (dofs 0 to 3)
1677 // 0 and 1 are the edges in the y dir. (sigma is function of y,
1678 // lambda is function of x). 2 and 3 are the edges in the x dir.
1679 // (sigma is function of x, lambda is function of y).
1680 //
1681 // More more info: see GeometryInfo for picture of the standard
1682 // element.
1683 //
1684 // Fill edge-based points:
1685 // std::vector<std::vector< Point<dim> > >
1686 // edge_points(lines_per_cell, std::vector<Point<dim>>
1687 // (n_q_points));
1688
1689 std::vector<std::vector<double>> edge_sigma_values(
1690 fe_data.edge_sigma_values);
1691 std::vector<std::vector<double>> edge_sigma_grads(
1692 fe_data.edge_sigma_grads);
1693
1694 std::vector<std::vector<double>> edge_lambda_values(
1695 fe_data.edge_lambda_values);
1696 std::vector<std::vector<double>> edge_lambda_grads(
1697 fe_data.edge_lambda_grads_2d);
1698
1699 // Adjust the edge_sigma_* for the current cell:
1700 for (unsigned int m = 0; m < lines_per_cell; ++m)
1701 {
1702 std::transform(edge_sigma_values[m].begin(),
1703 edge_sigma_values[m].end(),
1704 edge_sigma_values[m].begin(),
1705 [&](const double edge_sigma_value) {
1706 return edge_sign[m] * edge_sigma_value;
1707 });
1708
1709 std::transform(edge_sigma_grads[m].begin(),
1710 edge_sigma_grads[m].end(),
1711 edge_sigma_grads[m].begin(),
1712 [&](const double edge_sigma_grad) {
1713 return edge_sign[m] * edge_sigma_grad;
1714 });
1715 }
1716
1717 // If we want to generate shape gradients then we need second
1718 // derivatives of the 1d polynomials, but only first derivatives
1719 // for the shape values.
1720 const unsigned int poly_length =
1721 (flags & update_hessians) ?
1722 4 :
1723 ((flags & update_gradients) ? 3 : 2);
1724
1725
1726 for (unsigned int m = 0; m < lines_per_cell; ++m)
1727 {
1728 const unsigned int shift_m(m * this->n_dofs_per_line());
1729 for (unsigned int q = 0; q < n_q_points; ++q)
1730 {
1731 // Only compute 1d polynomials if degree>0.
1732 std::vector<std::vector<double>> poly(
1733 degree, std::vector<double>(poly_length));
1734 for (unsigned int i = 1; i < degree + 1; ++i)
1735 {
1736 // Compute all required 1d polynomials and their
1737 // derivatives, starting at degree 2. e.g. to access
1738 // L'_{i+2}(edge_sigma) use polyx[i][1].
1739 IntegratedLegendrePolynomials[i + 1].value(
1740 edge_sigma_values[m][q], poly[i - 1]);
1741 }
1742 if (flags & update_values)
1743 {
1744 // Lowest order edge shape functions:
1745 for (unsigned int d = 0; d < dim; ++d)
1746 {
1747 fe_data.shape_values[shift_m][q][d] =
1748 0.5 * edge_sigma_grads[m][d] *
1749 edge_lambda_values[m][q];
1750 }
1751 // Higher order edge shape functions:
1752 for (unsigned int i = 1; i < degree + 1; ++i)
1753 {
1754 const unsigned int poly_index = i - 1;
1755 const unsigned int dof_index(i + shift_m);
1756 for (unsigned int d = 0; d < dim; ++d)
1757 {
1758 fe_data.shape_values[dof_index][q][d] =
1759 edge_sigma_grads[m][d] *
1760 poly[poly_index][1] *
1761 edge_lambda_values[m][q] +
1762 poly[poly_index][0] *
1763 edge_lambda_grads[m][d];
1764 }
1765 }
1766 }
1767 if (flags & update_gradients)
1768 {
1769 // Lowest order edge shape functions:
1770 for (unsigned int d1 = 0; d1 < dim; ++d1)
1771 {
1772 for (unsigned int d2 = 0; d2 < dim; ++d2)
1773 {
1774 // Note: gradient is constant for a given
1775 // edge.
1776 fe_data.shape_grads[shift_m][q][d1][d2] =
1777 0.5 * edge_sigma_grads[m][d1] *
1778 edge_lambda_grads[m][d2];
1779 }
1780 }
1781 // Higher order edge shape functions:
1782 for (unsigned int i = 1; i < degree + 1; ++i)
1783 {
1784 const unsigned int poly_index = i - 1;
1785 const unsigned int dof_index(i + shift_m);
1786
1787 fe_data.shape_grads[dof_index][q][0][0] =
1788 edge_sigma_grads[m][0] *
1789 edge_sigma_grads[m][0] *
1790 edge_lambda_values[m][q] * poly[poly_index][2];
1791
1792 fe_data.shape_grads[dof_index][q][0][1] =
1793 (edge_sigma_grads[m][0] *
1794 edge_lambda_grads[m][1] +
1795 edge_sigma_grads[m][1] *
1796 edge_lambda_grads[m][0]) *
1797 poly[poly_index][1];
1798
1799 fe_data.shape_grads[dof_index][q][1][0] =
1800 fe_data.shape_grads[dof_index][q][0][1];
1801
1802 fe_data.shape_grads[dof_index][q][1][1] =
1803 edge_sigma_grads[m][1] *
1804 edge_sigma_grads[m][1] *
1805 edge_lambda_values[m][q] * poly[poly_index][2];
1806 }
1807 }
1808 if (flags & update_hessians)
1809 {
1810 // Lowest order edge shape function
1811 for (unsigned int d1 = 0; d1 < dim; ++d1)
1812 {
1813 for (unsigned int d2 = 0; d2 < dim; ++d2)
1814 {
1815 for (unsigned int d3 = 0; d3 < dim; ++d3)
1816 {
1817 fe_data.shape_hessians[shift_m][q][d1][d2]
1818 [d3] = 0;
1819 }
1820 }
1821 }
1822
1823 // Higher order edge shape function
1824 for (unsigned int i = 0; i < degree; ++i)
1825 {
1826 const unsigned int dof_index(i + 1 + shift_m);
1827
1828 for (unsigned int d1 = 0; d1 < dim; ++d1)
1829 {
1830 for (unsigned int d2 = 0; d2 < dim; ++d2)
1831 {
1832 for (unsigned int d3 = 0; d3 < dim; ++d3)
1833 {
1834 fe_data.shape_hessians[dof_index][q]
1835 [d1][d2][d3] =
1836 edge_sigma_grads[m][d1] *
1837 edge_sigma_grads[m][d2] *
1838 edge_sigma_grads[m][d3] *
1839 poly[i][3] *
1840 edge_lambda_values[m][q] +
1841 poly[i][2] *
1842 (edge_sigma_grads[m][d1] *
1843 edge_sigma_grads[m][d2] *
1844 edge_lambda_grads[m][d3] +
1845 edge_sigma_grads[m][d3] *
1846 edge_sigma_grads[m][d1] *
1847 edge_lambda_grads[m][d2] +
1848 edge_sigma_grads[m][d3] *
1849 edge_sigma_grads[m][d2] *
1850 edge_lambda_grads[m][d1]);
1851 }
1852 }
1853 }
1854 }
1855 }
1856 }
1857 }
1858 }
1859 break;
1860 }
1861 case 3:
1862 {
1864 {
1865 // Define an edge numbering so that each edge, E_{m} = [e^{m}_{1},
1866 // e^{m}_{2}] e1 = higher global numbering of the two local
1867 // vertices e2 = lower global numbering of the two local vertices
1868 std::vector<int> edge_sign(lines_per_cell);
1869 for (unsigned int m = 0; m < lines_per_cell; ++m)
1870 {
1871 unsigned int v0_loc =
1873 unsigned int v1_loc =
1875 unsigned int v0_glob = cell->vertex_index(v0_loc);
1876 unsigned int v1_glob = cell->vertex_index(v1_loc);
1877
1878 if (v0_glob > v1_glob)
1879 {
1880 // Opposite to global numbering on our reference element
1881 edge_sign[m] = -1.0;
1882 }
1883 else
1884 {
1885 // Aligns with global numbering on our reference element.
1886 edge_sign[m] = 1.0;
1887 }
1888 }
1889
1890 // Check for hanging faces. If we encounter hanging faces,
1891 // we use the vertex indices from the parent.
1892 for (unsigned int f = 0; f < 6; ++f)
1893 if (cell->face(f)->at_boundary() == false)
1894 if (cell->neighbor_is_coarser(f))
1895 for (unsigned int m = 0; m < 4; ++m)
1896 {
1897 unsigned int parent_m =
1899
1900 unsigned int v0_loc =
1902 unsigned int v1_loc =
1904
1905 unsigned int v0_glob =
1906 cell->parent()->vertex_index(v0_loc);
1907 unsigned int v1_glob =
1908 cell->parent()->vertex_index(v1_loc);
1909
1910 if (v0_glob > v1_glob)
1911 {
1912 // Opposite to global numbering on our reference
1913 // element
1914 edge_sign[parent_m] = -1.0;
1915 }
1916 else
1917 {
1918 // Aligns with global numbering on our reference
1919 // element.
1920 edge_sign[parent_m] = 1.0;
1921 }
1922 }
1923
1924 // Next, we cover the case where we encounter a hanging edge but
1925 // not a hanging face. This is always the case if the cell
1926 // currently considered shares all faces with cells of the same
1927 // refinement level but shares one edge with a coarser cell,
1928 // e.g.:
1929 // *----*----*---------*
1930 // / / / / |
1931 // *----*----* / |
1932 // / / / / |
1933 // *----*----*---------* *
1934 // | | | | / |
1935 // *----*----* | * *
1936 // | | | | / | / |
1937 // *----*----*----*----* * *
1938 // | | | | | / | /
1939 // *----*----*----*----* *
1940 // | | | | | /
1941 // *----*----*----*----*
1942 // where the cell at the left bottom is the currently
1943 // considered cell.
1944 //
1945 // In that case, we determine the direction of the hanging edge
1946 // based on the vertex indices from the parent cell.
1947 //
1948 // Note: We assume here that at most eight cells are adjacent to
1949 // a single vertex.
1950 std::vector<unsigned int> adjacent_faces = {2, 2, 4, 4, 0, 0};
1951 for (unsigned int f = 0; f < 6; ++f)
1952 {
1953 if (!cell->face(f)->at_boundary() &&
1954 !cell->face(adjacent_faces[f])->at_boundary())
1955 if (!cell->neighbor_is_coarser(f) &&
1956 !cell->neighbor_is_coarser(adjacent_faces[f]))
1957 if (!cell->neighbor(f)
1958 ->face(adjacent_faces[f])
1959 ->at_boundary())
1960 if (cell->neighbor(f)->neighbor_is_coarser(
1961 adjacent_faces[f]))
1962 {
1963 unsigned int parent_m =
1965 unsigned int v0_loc =
1967 0);
1968 unsigned int v1_loc =
1970 1);
1971 unsigned int v0_glob =
1972 cell->parent()->vertex_index(v0_loc);
1973 unsigned int v1_glob =
1974 cell->parent()->vertex_index(v1_loc);
1975 if (v0_glob > v1_glob)
1976 {
1977 // Opposite to global numbering on our reference
1978 // element
1979 edge_sign[parent_m] = -1.0;
1980 }
1981 else
1982 {
1983 // Aligns with global numbering on our reference
1984 // element.
1985 edge_sign[parent_m] = 1.0;
1986 }
1987 }
1988
1989 if (!cell->face(f)->at_boundary() &&
1990 !cell->face(adjacent_faces[f] + 1)->at_boundary())
1991 if (!cell->neighbor_is_coarser(f) &&
1992 !cell->neighbor_is_coarser(adjacent_faces[f] + 1))
1993 if (!cell->neighbor(f)
1994 ->face(adjacent_faces[f] + 1)
1995 ->at_boundary())
1996 if (cell->neighbor(f)->neighbor_is_coarser(
1997 adjacent_faces[f] + 1))
1998 {
1999 unsigned int parent_m =
2001 unsigned int v0_loc =
2003 0);
2004 unsigned int v1_loc =
2006 1);
2007 unsigned int v0_glob =
2008 cell->parent()->vertex_index(v0_loc);
2009 unsigned int v1_glob =
2010 cell->parent()->vertex_index(v1_loc);
2011 if (v0_glob > v1_glob)
2012 {
2013 // Opposite to global numbering on our reference
2014 // element
2015 edge_sign[parent_m] = -1.0;
2016 }
2017 else
2018 {
2019 // Aligns with global numbering on our reference
2020 // element.
2021 edge_sign[parent_m] = 1.0;
2022 }
2023 }
2024 }
2025
2026 // Define \sigma_{m} = sigma_{e^{m}_{1}} - sigma_{e^{m}_{2}}
2027 // \lambda_{m} = \lambda_{e^{m}_{1}} + \lambda_{e^{m}_{2}}
2028 //
2029 // To help things, in fe_data, we have precomputed (sigma_{i} -
2030 // sigma_{j}) and (lambda_{i} + lambda_{j}) for 0<= i,j <
2031 // lines_per_cell.
2032 //
2033 // There are two types:
2034 // - lower order (1 per edge, m):
2035 // \phi_{m}^{\mathcal{N}_{0}} = 1/2 grad(\sigma_{m})\lambda_{m}
2036 //
2037 // - higher order (degree per edge, m):
2038 // \phi_{i}^{E_{m}} = grad( L_{i+2}(\sigma_{m}) (\lambda_{m}) ).
2039 //
2040 // NOTE: In the ref cell, sigma_{m} is a function of x OR y OR Z
2041 // and lambda_{m} a function of the remaining co-ords.
2042 // for example, if sigma is of x, then lambda is of y AND
2043 // z, and so on. This means that grad(\sigma) requires
2044 // multiplication by d(sigma)/dx_{i} for the i^th comp of
2045 // grad(sigma) and similarly when taking derivatives of
2046 // lambda.
2047 //
2048 // First handle the lowest order edges (dofs 0 to 11)
2049 // 0 and 1 are the edges in the y dir at z=0. (sigma is a fn of y,
2050 // lambda is a fn of x & z). 2 and 3 are the edges in the x dir at
2051 // z=0. (sigma is a fn of x, lambda is a fn of y & z). 4 and 5 are
2052 // the edges in the y dir at z=1. (sigma is a fn of y, lambda is a
2053 // fn of x & z). 6 and 7 are the edges in the x dir at z=1. (sigma
2054 // is a fn of x, lambda is a fn of y & z). 8 and 9 are the edges
2055 // in the z dir at y=0. (sigma is a fn of z, lambda is a fn of x &
2056 // y). 10 and 11 are the edges in the z dir at y=1. (sigma is a fn
2057 // of z, lambda is a fn of x & y).
2058 //
2059 // For more info: see GeometryInfo for picture of the standard
2060 // element.
2061
2062 // Copy over required edge-based data:
2063 std::vector<std::vector<double>> edge_sigma_values(
2064 fe_data.edge_sigma_values);
2065 std::vector<std::vector<double>> edge_lambda_values(
2066 fe_data.edge_lambda_values);
2067 std::vector<std::vector<double>> edge_sigma_grads(
2068 fe_data.edge_sigma_grads);
2069 std::vector<std::vector<std::vector<double>>> edge_lambda_grads(
2070 fe_data.edge_lambda_grads_3d);
2071 std::vector<std::vector<std::vector<double>>>
2072 edge_lambda_gradgrads_3d(fe_data.edge_lambda_gradgrads_3d);
2073
2074 // Adjust the edge_sigma_* for the current cell:
2075 for (unsigned int m = 0; m < lines_per_cell; ++m)
2076 {
2077 std::transform(edge_sigma_values[m].begin(),
2078 edge_sigma_values[m].end(),
2079 edge_sigma_values[m].begin(),
2080 [&](const double edge_sigma_value) {
2081 return edge_sign[m] * edge_sigma_value;
2082 });
2083 std::transform(edge_sigma_grads[m].begin(),
2084 edge_sigma_grads[m].end(),
2085 edge_sigma_grads[m].begin(),
2086 [&](const double edge_sigma_grad) {
2087 return edge_sign[m] * edge_sigma_grad;
2088 });
2089 }
2090
2091 // Now calculate the edge-based shape functions:
2092 // If we want to generate shape gradients then we need second
2093 // derivatives of the 1d polynomials, but only first derivatives
2094 // for the shape values.
2095 const unsigned int poly_length =
2096 (flags & update_hessians) ?
2097 4 :
2098 ((flags & update_gradients) ? 3 : 2);
2099
2100 std::vector<std::vector<double>> poly(
2101 degree, std::vector<double>(poly_length));
2102 for (unsigned int m = 0; m < lines_per_cell; ++m)
2103 {
2104 const unsigned int shift_m(m * this->n_dofs_per_line());
2105 for (unsigned int q = 0; q < n_q_points; ++q)
2106 {
2107 // precompute values of all 1d polynomials required:
2108 // for example poly[i][1] = L'_{i+2}(edge_sigma_values)
2109 if (degree > 0)
2110 {
2111 for (unsigned int i = 0; i < degree; ++i)
2112 {
2113 IntegratedLegendrePolynomials[i + 2].value(
2114 edge_sigma_values[m][q], poly[i]);
2115 }
2116 }
2117 if (flags & update_values)
2118 {
2119 // Lowest order edge shape functions:
2120 for (unsigned int d = 0; d < dim; ++d)
2121 {
2122 fe_data.shape_values[shift_m][q][d] =
2123 0.5 * edge_sigma_grads[m][d] *
2124 edge_lambda_values[m][q];
2125 }
2126 // Higher order edge shape functions
2127 if (degree > 0)
2128 {
2129 for (unsigned int i = 0; i < degree; ++i)
2130 {
2131 const unsigned int dof_index(i + 1 + shift_m);
2132 for (unsigned int d = 0; d < dim; ++d)
2133 {
2134 fe_data.shape_values[dof_index][q][d] =
2135 edge_sigma_grads[m][d] * poly[i][1] *
2136 edge_lambda_values[m][q] +
2137 poly[i][0] * edge_lambda_grads[m][q][d];
2138 }
2139 }
2140 }
2141 }
2142 if (flags & update_gradients)
2143 {
2144 // Lowest order edge shape functions:
2145 for (unsigned int d1 = 0; d1 < dim; ++d1)
2146 {
2147 for (unsigned int d2 = 0; d2 < dim; ++d2)
2148 {
2149 fe_data.shape_grads[shift_m][q][d1][d2] =
2150 0.5 * edge_sigma_grads[m][d1] *
2151 edge_lambda_grads[m][q][d2];
2152 }
2153 }
2154 // Higher order edge shape functions
2155 if (degree > 0)
2156 {
2157 for (unsigned int i = 0; i < degree; ++i)
2158 {
2159 const unsigned int dof_index(i + 1 + shift_m);
2160
2161 for (unsigned int d1 = 0; d1 < dim; ++d1)
2162 {
2163 for (unsigned int d2 = 0; d2 < dim; ++d2)
2164 {
2165 fe_data
2166 .shape_grads[dof_index][q][d1][d2] =
2167 edge_sigma_grads[m][d1] *
2168 edge_sigma_grads[m][d2] *
2169 poly[i][2] *
2170 edge_lambda_values[m][q] +
2171 (edge_sigma_grads[m][d1] *
2172 edge_lambda_grads[m][q][d2] +
2173 edge_sigma_grads[m][d2] *
2174 edge_lambda_grads[m][q][d1]) *
2175 poly[i][1] +
2176 edge_lambda_gradgrads_3d[m][d1]
2177 [d2] *
2178 poly[i][0];
2179 }
2180 }
2181 }
2182 }
2183 }
2184 if (flags & update_hessians)
2185 {
2186 // Lowest order edge shape functions:
2187 for (unsigned int d1 = 0; d1 < dim; ++d1)
2188 {
2189 for (unsigned int d2 = 0; d2 < dim; ++d2)
2190 {
2191 for (unsigned int d3 = 0; d3 < dim; ++d3)
2192 {
2193 fe_data.shape_hessians[shift_m][q][d1][d2]
2194 [d3] =
2195 0.5 * edge_sigma_grads[m][d1] *
2196 edge_lambda_gradgrads_3d[m][d3][d2];
2197 }
2198 }
2199 }
2200
2201 // Higher order edge shape functions
2202 if (degree > 0)
2203 {
2204 for (unsigned int i = 0; i < degree; ++i)
2205 {
2206 const unsigned int dof_index(i + 1 + shift_m);
2207
2208 for (unsigned int d1 = 0; d1 < dim; ++d1)
2209 {
2210 for (unsigned int d2 = 0; d2 < dim; ++d2)
2211 {
2212 for (unsigned int d3 = 0; d3 < dim;
2213 ++d3)
2214 {
2215 fe_data
2216 .shape_hessians[dof_index][q]
2217 [d1][d2][d3] =
2218 edge_sigma_grads[m][d1] *
2219 edge_sigma_grads[m][d2] *
2220 edge_sigma_grads[m][d3] *
2221 poly[i][3] *
2222 edge_lambda_values[m][q] +
2223 poly[i][2] *
2224 (edge_sigma_grads[m][d1] *
2225 edge_sigma_grads[m][d2] *
2226 edge_lambda_grads[m][q]
2227 [d3] +
2228 edge_sigma_grads[m][d3] *
2229 edge_sigma_grads[m][d1] *
2230 edge_lambda_grads[m][q]
2231 [d2] +
2232 edge_sigma_grads[m][d3] *
2233 edge_sigma_grads[m][d2] *
2234 edge_lambda_grads[m][q]
2235 [d1]) +
2236 poly[i][1] *
2237 (edge_sigma_grads[m][d1] *
2238 edge_lambda_gradgrads_3d
2239 [m][d3][d2] +
2240 edge_sigma_grads[m][d2] *
2241 edge_lambda_gradgrads_3d
2242 [m][d3][d1] +
2243 edge_sigma_grads[m][d3] *
2244 edge_lambda_gradgrads_3d
2245 [m][d2][d1]);
2246 }
2247 }
2248 }
2249 }
2250 }
2251 }
2252 }
2253 }
2254 }
2255 break;
2256 }
2257 default:
2258 {
2260 }
2261 }
2262}
2263
2264
2265
2266template <int dim, int spacedim>
2267void
2269 const typename Triangulation<dim, dim>::cell_iterator &cell,
2270 const Quadrature<dim> &quadrature,
2271 const InternalData &fe_data) const
2272{
2273 // This function handles the cell-dependent construction of the FACE-based
2274 // shape functions.
2275 //
2276 // Note that it should only be called in 3d.
2277 Assert(dim == 3, ExcDimensionMismatch(dim, 3));
2278 //
2279 // It will fill in the missing parts of fe_data which were not possible to
2280 // fill in the get_data routine, with respect to face-based shape functions.
2281 //
2282 // It should be called by the fill_fe_*_values routines in order to complete
2283 // the basis set at quadrature points on the current cell for each face.
2284
2285 // Useful constants:
2286 const unsigned int degree(
2287 this->degree -
2288 1); // Note: constructor takes input degree + 1, so need to knock 1 off.
2289
2290 // Do nothing if FE degree is 0.
2291 if (degree > 0)
2292 {
2293 const UpdateFlags flags(fe_data.update_each);
2294
2296 {
2297 const unsigned int n_q_points = quadrature.size();
2298
2299 Assert(!(flags & update_values) ||
2300 fe_data.shape_values.size() == this->n_dofs_per_cell(),
2301 ExcDimensionMismatch(fe_data.shape_values.size(),
2302 this->n_dofs_per_cell()));
2303 Assert(!(flags & update_values) ||
2304 fe_data.shape_values[0].size() == n_q_points,
2305 ExcDimensionMismatch(fe_data.shape_values[0].size(),
2306 n_q_points));
2307
2308 // Useful geometry info:
2309 const unsigned int vertices_per_face(
2311 const unsigned int faces_per_cell(GeometryInfo<dim>::faces_per_cell);
2312
2313 // DoF info:
2314 const unsigned int n_line_dofs =
2315 this->n_dofs_per_line() * GeometryInfo<dim>::lines_per_cell;
2316
2317 // First we find the global face orientations on the current cell.
2318 std::vector<std::vector<unsigned int>> face_orientation(
2319 faces_per_cell, std::vector<unsigned int>(vertices_per_face));
2320
2321 const unsigned int
2322 vertex_opposite_on_face[GeometryInfo<3>::vertices_per_face] = {3,
2323 2,
2324 1,
2325 0};
2326
2327 const unsigned int
2328 vertices_adjacent_on_face[GeometryInfo<3>::vertices_per_face][2] = {
2329 {1, 2}, {0, 3}, {3, 0}, {2, 1}};
2330
2331 for (unsigned int m = 0; m < faces_per_cell; ++m)
2332 {
2333 // Check, if we are on a hanging face.
2334 bool cell_has_coarser_neighbor = false;
2335 if (cell->face(m)->at_boundary() == false)
2336 if (cell->neighbor_is_coarser(m))
2337 cell_has_coarser_neighbor = true;
2338
2339 // Find the local vertex on this face with the highest global
2340 // numbering. This is f^m_0.
2341 unsigned int current_max = 0;
2342
2343 // We start with the hanging face case, where the face
2344 // orientation is determined based on the vertex indices
2345 // of the parent cell.
2346 if (cell_has_coarser_neighbor)
2347 {
2348 unsigned int current_glob = cell->parent()->vertex_index(
2350 for (unsigned int v = 1; v < vertices_per_face; ++v)
2351 {
2352 if (current_glob <
2353 cell->parent()->vertex_index(
2355 {
2356 current_max = v;
2357 current_glob = cell->parent()->vertex_index(
2359 }
2360 }
2361 }
2362 // Otherwise, the face orientation is based on its own
2363 // vertex indices.
2364 else
2365 {
2366 unsigned int current_glob = cell->vertex_index(
2368 for (unsigned int v = 1; v < vertices_per_face; ++v)
2369 {
2370 if (current_glob <
2371 cell->vertex_index(
2373 {
2374 current_max = v;
2375 current_glob = cell->vertex_index(
2377 }
2378 }
2379 }
2380
2381 face_orientation[m][0] =
2383
2384 // f^m_2 is the vertex opposite f^m_0.
2385 face_orientation[m][2] = GeometryInfo<dim>::face_to_cell_vertices(
2386 m, vertex_opposite_on_face[current_max]);
2387
2388 // Finally, f^m_1 is the vertex with the greater global numbering
2389 // of the remaining two local vertices. Then, f^m_3 is the other.
2390 // Again, we need to distinguish between the hanging face and the
2391 // non-hanging face cases. In the case of hanging faces, we
2392 // consider the vertex indices from the parent. Otherwise, we
2393 // consider the vertex indices of the face itself.
2394 if (cell_has_coarser_neighbor)
2395 {
2396 if (cell->parent()->vertex_index(
2398 m, vertices_adjacent_on_face[current_max][0])) >
2399 cell->parent()->vertex_index(
2401 m, vertices_adjacent_on_face[current_max][1])))
2402 {
2403 face_orientation[m][1] =
2405 m, vertices_adjacent_on_face[current_max][0]);
2406 face_orientation[m][3] =
2408 m, vertices_adjacent_on_face[current_max][1]);
2409 }
2410 else
2411 {
2412 face_orientation[m][1] =
2414 m, vertices_adjacent_on_face[current_max][1]);
2415 face_orientation[m][3] =
2417 m, vertices_adjacent_on_face[current_max][0]);
2418 }
2419 }
2420 else
2421 {
2422 if (cell->vertex_index(
2424 m, vertices_adjacent_on_face[current_max][0])) >
2425 cell->vertex_index(
2427 m, vertices_adjacent_on_face[current_max][1])))
2428 {
2429 face_orientation[m][1] =
2431 m, vertices_adjacent_on_face[current_max][0]);
2432 face_orientation[m][3] =
2434 m, vertices_adjacent_on_face[current_max][1]);
2435 }
2436 else
2437 {
2438 face_orientation[m][1] =
2440 m, vertices_adjacent_on_face[current_max][1]);
2441 face_orientation[m][3] =
2443 m, vertices_adjacent_on_face[current_max][0]);
2444 }
2445 }
2446 }
2447 // Now we know the face orientation on the current cell, we can
2448 // generate the parameterisation:
2449 std::vector<std::vector<double>> face_xi_values(
2450 faces_per_cell, std::vector<double>(n_q_points));
2451 std::vector<std::vector<double>> face_xi_grads(
2452 faces_per_cell, std::vector<double>(dim));
2453 std::vector<std::vector<double>> face_eta_values(
2454 faces_per_cell, std::vector<double>(n_q_points));
2455 std::vector<std::vector<double>> face_eta_grads(
2456 faces_per_cell, std::vector<double>(dim));
2457
2458 std::vector<std::vector<double>> face_lambda_values(
2459 faces_per_cell, std::vector<double>(n_q_points));
2460 std::vector<std::vector<double>> face_lambda_grads(
2461 faces_per_cell, std::vector<double>(dim));
2462 for (unsigned int m = 0; m < faces_per_cell; ++m)
2463 {
2464 for (unsigned int q = 0; q < n_q_points; ++q)
2465 {
2466 face_xi_values[m][q] =
2467 fe_data.sigma_imj_values[q][face_orientation[m][0]]
2468 [face_orientation[m][1]];
2469 face_eta_values[m][q] =
2470 fe_data.sigma_imj_values[q][face_orientation[m][0]]
2471 [face_orientation[m][3]];
2472 face_lambda_values[m][q] = fe_data.face_lambda_values[m][q];
2473 }
2474 for (unsigned int d = 0; d < dim; ++d)
2475 {
2476 face_xi_grads[m][d] =
2477 fe_data.sigma_imj_grads[face_orientation[m][0]]
2478 [face_orientation[m][1]][d];
2479 face_eta_grads[m][d] =
2480 fe_data.sigma_imj_grads[face_orientation[m][0]]
2481 [face_orientation[m][3]][d];
2482
2483 face_lambda_grads[m][d] = fe_data.face_lambda_grads[m][d];
2484 }
2485 }
2486 // Now can generate the basis
2487 const unsigned int poly_length =
2488 (flags & update_hessians) ? 4 :
2489 ((flags & update_gradients) ? 3 : 2);
2490
2491
2492 std::vector<std::vector<double>> polyxi(
2493 degree, std::vector<double>(poly_length));
2494 std::vector<std::vector<double>> polyeta(
2495 degree, std::vector<double>(poly_length));
2496
2497 // Loop through quad points:
2498 for (unsigned int m = 0; m < faces_per_cell; ++m)
2499 {
2500 // we assume that all quads have the same number of dofs
2501 const unsigned int shift_m(m * this->n_dofs_per_quad(0));
2502 // Calculate the offsets for each face-based shape function:
2503 //
2504 // Type-1 (gradients)
2505 // \phi^{F_m,1}_{ij} = \nabla( L_{i+2}(\xi_{F_{m}})
2506 // L_{j+2}(\eta_{F_{m}}) \lambda_{F_{m}} )
2507 //
2508 // 0 <= i,j < degree (in a group of degree*degree)
2509 const unsigned int face_type1_offset(n_line_dofs + shift_m);
2510 // Type-2:
2511 //
2512 // \phi^{F_m,2}_{ij} = ( L'_{i+2}(\xi_{F_{m}})
2513 // L_{j+2}(\eta_{F_{m}}) \nabla\xi_{F_{m}}
2514 // - L_{i+2}(\xi_{F_{m}})
2515 // L'_{j+2}(\eta_{F_{m}}) \nabla\eta_{F_{m}}
2516 // ) \lambda_{F_{m}}
2517 //
2518 // 0 <= i,j < degree (in a group of degree*degree)
2519 const unsigned int face_type2_offset(face_type1_offset +
2520 degree * degree);
2521 // Type-3:
2522 //
2523 // \phi^{F_m,3}_{i} = L_{i+2}(\eta_{F_{m}}) \lambda_{F_{m}}
2524 // \nabla\xi_{F_{m}}
2525 // \phi^{F_m,3}_{i+p} = L_{i+2}(\xi_{F_{m}})
2526 // \lambda_{F_{m}} \nabla\eta_{F_{m}}
2527 //
2528 // 0 <= i < degree.
2529 //
2530 // here we order so that all of subtype 1 comes first, then
2531 // subtype 2.
2532 const unsigned int face_type3_offset1(face_type2_offset +
2533 degree * degree);
2534 const unsigned int face_type3_offset2(face_type3_offset1 +
2535 degree);
2536
2537 // Loop over all faces:
2538 for (unsigned int q = 0; q < n_q_points; ++q)
2539 {
2540 // pre-compute values & required derivatives at this quad
2541 // point: polyxi = L_{i+2}(\xi_{F_{m}}), polyeta =
2542 // L_{j+2}(\eta_{F_{m}}),
2543 //
2544 // each polypoint[k][d], contains the dth derivative of
2545 // L_{k+2} at the point \xi or \eta. Note that this doesn't
2546 // include the derivative of xi/eta via the chain rule.
2547 for (unsigned int i = 0; i < degree; ++i)
2548 {
2549 // compute all required 1d polynomials:
2550 IntegratedLegendrePolynomials[i + 2].value(
2551 face_xi_values[m][q], polyxi[i]);
2552 IntegratedLegendrePolynomials[i + 2].value(
2553 face_eta_values[m][q], polyeta[i]);
2554 }
2555 // Now use these to compute the shape functions:
2556 if (flags & update_values)
2557 {
2558 for (unsigned int j = 0; j < degree; ++j)
2559 {
2560 const unsigned int shift_j(j * degree);
2561 for (unsigned int i = 0; i < degree; ++i)
2562 {
2563 const unsigned int shift_ij(shift_j + i);
2564 // Type 1:
2565 const unsigned int dof_index1(face_type1_offset +
2566 shift_ij);
2567 for (unsigned int d = 0; d < dim; ++d)
2568 {
2569 fe_data.shape_values[dof_index1][q][d] =
2570 (face_xi_grads[m][d] * polyxi[i][1] *
2571 polyeta[j][0] +
2572 face_eta_grads[m][d] * polyxi[i][0] *
2573 polyeta[j][1]) *
2574 face_lambda_values[m][q] +
2575 face_lambda_grads[m][d] * polyxi[i][0] *
2576 polyeta[j][0];
2577 }
2578 // Type 2:
2579 const unsigned int dof_index2(face_type2_offset +
2580 shift_ij);
2581 for (unsigned int d = 0; d < dim; ++d)
2582 {
2583 fe_data.shape_values[dof_index2][q][d] =
2584 (face_xi_grads[m][d] * polyxi[i][1] *
2585 polyeta[j][0] -
2586 face_eta_grads[m][d] * polyxi[i][0] *
2587 polyeta[j][1]) *
2588 face_lambda_values[m][q];
2589 }
2590 }
2591 // Type 3:
2592 const unsigned int dof_index3_1(face_type3_offset1 +
2593 j);
2594 const unsigned int dof_index3_2(face_type3_offset2 +
2595 j);
2596 for (unsigned int d = 0; d < dim; ++d)
2597 {
2598 fe_data.shape_values[dof_index3_1][q][d] =
2599 face_xi_grads[m][d] * polyeta[j][0] *
2600 face_lambda_values[m][q];
2601
2602 fe_data.shape_values[dof_index3_2][q][d] =
2603 face_eta_grads[m][d] * polyxi[j][0] *
2604 face_lambda_values[m][q];
2605 }
2606 }
2607 }
2608 if (flags & update_gradients)
2609 {
2610 for (unsigned int j = 0; j < degree; ++j)
2611 {
2612 const unsigned int shift_j(j * degree);
2613 for (unsigned int i = 0; i < degree; ++i)
2614 {
2615 const unsigned int shift_ij(shift_j + i);
2616 // Type 1:
2617 const unsigned int dof_index1(face_type1_offset +
2618 shift_ij);
2619 for (unsigned int d1 = 0; d1 < dim; ++d1)
2620 {
2621 for (unsigned int d2 = 0; d2 < dim; ++d2)
2622 {
2623 fe_data
2624 .shape_grads[dof_index1][q][d1][d2] =
2625 (face_xi_grads[m][d1] *
2626 face_xi_grads[m][d2] * polyxi[i][2] *
2627 polyeta[j][0] +
2628 (face_xi_grads[m][d1] *
2629 face_eta_grads[m][d2] +
2630 face_xi_grads[m][d2] *
2631 face_eta_grads[m][d1]) *
2632 polyxi[i][1] * polyeta[j][1] +
2633 face_eta_grads[m][d1] *
2634 face_eta_grads[m][d2] *
2635 polyxi[i][0] * polyeta[j][2]) *
2636 face_lambda_values[m][q] +
2637 (face_xi_grads[m][d2] * polyxi[i][1] *
2638 polyeta[j][0] +
2639 face_eta_grads[m][d2] * polyxi[i][0] *
2640 polyeta[j][1]) *
2641 face_lambda_grads[m][d1] +
2642 (face_xi_grads[m][d1] * polyxi[i][1] *
2643 polyeta[j][0] +
2644 face_eta_grads[m][d1] * polyxi[i][0] *
2645 polyeta[j][1]) *
2646 face_lambda_grads[m][d2];
2647 }
2648 }
2649 // Type 2:
2650 const unsigned int dof_index2(face_type2_offset +
2651 shift_ij);
2652 for (unsigned int d1 = 0; d1 < dim; ++d1)
2653 {
2654 for (unsigned int d2 = 0; d2 < dim; ++d2)
2655 {
2656 fe_data
2657 .shape_grads[dof_index2][q][d1][d2] =
2658 (face_xi_grads[m][d1] *
2659 face_xi_grads[m][d2] * polyxi[i][2] *
2660 polyeta[j][0] +
2661 (face_xi_grads[m][d1] *
2662 face_eta_grads[m][d2] -
2663 face_xi_grads[m][d2] *
2664 face_eta_grads[m][d1]) *
2665 polyxi[i][1] * polyeta[j][1] -
2666 face_eta_grads[m][d1] *
2667 face_eta_grads[m][d2] *
2668 polyxi[i][0] * polyeta[j][2]) *
2669 face_lambda_values[m][q] +
2670 (face_xi_grads[m][d1] * polyxi[i][1] *
2671 polyeta[j][0] -
2672 face_eta_grads[m][d1] * polyxi[i][0] *
2673 polyeta[j][1]) *
2674 face_lambda_grads[m][d2];
2675 }
2676 }
2677 }
2678 // Type 3:
2679 const unsigned int dof_index3_1(face_type3_offset1 +
2680 j);
2681 const unsigned int dof_index3_2(face_type3_offset2 +
2682 j);
2683 for (unsigned int d1 = 0; d1 < dim; ++d1)
2684 {
2685 for (unsigned int d2 = 0; d2 < dim; ++d2)
2686 {
2687 fe_data.shape_grads[dof_index3_1][q][d1][d2] =
2688 face_xi_grads[m][d1] *
2689 (face_eta_grads[m][d2] * polyeta[j][1] *
2690 face_lambda_values[m][q] +
2691 face_lambda_grads[m][d2] * polyeta[j][0]);
2692
2693 fe_data.shape_grads[dof_index3_2][q][d1][d2] =
2694 face_eta_grads[m][d1] *
2695 (face_xi_grads[m][d2] * polyxi[j][1] *
2696 face_lambda_values[m][q] +
2697 face_lambda_grads[m][d2] * polyxi[j][0]);
2698 }
2699 }
2700 }
2701 }
2702 if (flags & update_hessians)
2703 {
2704 for (unsigned int j = 0; j < degree; ++j)
2705 {
2706 const unsigned int shift_j(j * degree);
2707 for (unsigned int i = 0; i < degree; ++i)
2708 {
2709 const unsigned int shift_ij(shift_j + i);
2710
2711 // Type 1:
2712 const unsigned int dof_index1(face_type1_offset +
2713 shift_ij);
2714 for (unsigned int d1 = 0; d1 < dim; ++d1)
2715 {
2716 for (unsigned int d2 = 0; d2 < dim; ++d2)
2717 {
2718 for (unsigned int d3 = 0; d3 < dim; ++d3)
2719 {
2720 fe_data.shape_hessians[dof_index1][q]
2721 [d1][d2][d3] =
2722 polyxi[i][1] *
2723 face_xi_grads[m][d3] *
2724 (face_eta_grads[m][d1] *
2725 (polyeta[j][2] *
2726 face_eta_grads[m][d2] *
2727 face_lambda_values[m][q] +
2728 polyeta[j][1] *
2729 face_lambda_grads[m][d2]) +
2730 polyeta[j][1] *
2731 face_eta_grads[m][d2] *
2732 face_lambda_grads[m][d1]) +
2733 polyxi[i][0] *
2734 (polyeta[j][3] *
2735 face_eta_grads[m][d1] *
2736 face_eta_grads[m][d2] *
2737 face_eta_grads[m][d3] *
2738 face_lambda_values[m][q] +
2739 polyeta[j][2] *
2740 (face_eta_grads[m][d1] *
2741 face_eta_grads[m][d2] *
2742 face_lambda_grads[m][d3] +
2743 face_eta_grads[m][d3] *
2744 (face_eta_grads[m][d1] *
2745 face_lambda_grads[m]
2746 [d2] +
2747 face_eta_grads[m][d2] *
2748 face_lambda_grads
2749 [m][d1]))) +
2750 (polyxi[i][1] * polyeta[j][1] *
2751 face_eta_grads[m][d3] +
2752 polyxi[i][2] * polyeta[j][0] *
2753 face_xi_grads[m][d3]) *
2754 (face_xi_grads[m][d1] *
2755 face_lambda_grads[m][d2] +
2756 face_xi_grads[m][d2] *
2757 face_lambda_grads[m][d1]) +
2758 face_lambda_grads[m][d3] *
2759 (polyxi[i][2] * polyeta[j][0] *
2760 face_xi_grads[m][d1] *
2761 face_xi_grads[m][d2] +
2762 polyxi[i][1] * polyeta[j][1] *
2763 (face_xi_grads[m][d1] *
2764 face_eta_grads[m][d2] +
2765 face_xi_grads[m][d2] *
2766 face_eta_grads[m][d1])) +
2767 face_lambda_values[m][q] *
2768 (polyxi[i][3] * polyeta[j][0] *
2769 face_xi_grads[m][d1] *
2770 face_xi_grads[m][d2] *
2771 face_xi_grads[m][d3] +
2772 polyxi[i][1] * polyeta[j][2] *
2773 face_eta_grads[m][d3] *
2774 (face_xi_grads[m][d1] *
2775 face_eta_grads[m][d2] +
2776 face_xi_grads[m][d2] *
2777 face_eta_grads[m][d1]) +
2778 polyxi[i][2] * polyeta[j][1] *
2779 (face_xi_grads[m][d3] *
2780 face_xi_grads[m][d2] *
2781 face_eta_grads[m][d1] +
2782 face_xi_grads[m][d1] *
2783 (face_xi_grads[m][d2] *
2784 face_eta_grads[m][d3] +
2785 face_xi_grads[m][d3] *
2786 face_eta_grads[m][d2])));
2787 }
2788 }
2789 }
2790
2791 // Type 2:
2792 const unsigned int dof_index2(face_type2_offset +
2793 shift_ij);
2794 for (unsigned int d1 = 0; d1 < dim; ++d1)
2795 {
2796 for (unsigned int d2 = 0; d2 < dim; ++d2)
2797 {
2798 for (unsigned int d3 = 0; d3 < dim; ++d3)
2799 {
2800 fe_data.shape_hessians[dof_index2][q]
2801 [d1][d2][d3] =
2802 face_xi_grads[m][d1] *
2803 (polyxi[i][1] * polyeta[j][1] *
2804 (face_eta_grads[m][d2] *
2805 face_lambda_grads[m][d3] +
2806 face_eta_grads[m][d3] *
2807 face_lambda_grads[m][d2]) +
2808 polyxi[i][2] * polyeta[j][0] *
2809 (face_xi_grads[m][d2] *
2810 face_lambda_grads[m][d3] +
2811 face_xi_grads[m][d3] *
2812 face_lambda_grads[m][d2]) +
2813 face_lambda_values[m][q] *
2814 (face_eta_grads[m][d2] *
2815 (polyxi[i][1] *
2816 polyeta[j][2] *
2817 face_eta_grads[m][d3] +
2818 polyxi[i][2] *
2819 polyeta[j][1] *
2820 face_xi_grads[m][d3]) +
2821 face_xi_grads[m][d2] *
2822 (polyxi[i][2] *
2823 polyeta[j][1] *
2824 face_eta_grads[m][d3] +
2825 polyxi[i][3] *
2826 polyeta[j][0] *
2827 face_xi_grads[m][d3]))) -
2828 polyxi[i][0] *
2829 face_eta_grads[m][d1] *
2830 (face_eta_grads[m][d2] *
2831 (polyeta[j][3] *
2832 face_eta_grads[m][d3] *
2833 face_lambda_values[m][q] +
2834 polyeta[j][2] *
2835 face_lambda_grads[m][d3]) +
2836 polyeta[j][2] *
2837 face_eta_grads[m][d3] *
2838 face_lambda_grads[m][d2]) -
2839 face_eta_grads[m][d1] *
2840 (polyxi[i][1] *
2841 face_xi_grads[m][d3] *
2842 (polyeta[j][2] *
2843 face_eta_grads[m][d2] *
2844 face_lambda_values[m][q] +
2845 polyeta[j][1] *
2846 face_lambda_grads[m][d2]) +
2847 face_xi_grads[m][d2] *
2848 (polyxi[i][1] *
2849 (polyeta[j][2] *
2850 face_eta_grads[m][d3] *
2851 face_lambda_values[m]
2852 [q] +
2853 polyeta[j][1] *
2854 face_lambda_grads[m]
2855 [d3]) +
2856 polyxi[i][2] * polyeta[j][1] *
2857 face_xi_grads[m][d3] *
2858 face_lambda_values[m][q]));
2859 }
2860 }
2861 }
2862 }
2863 // Type 3:
2864 const unsigned int dof_index3_1(face_type3_offset1 +
2865 j);
2866 const unsigned int dof_index3_2(face_type3_offset2 +
2867 j);
2868 for (unsigned int d1 = 0; d1 < dim; ++d1)
2869 {
2870 for (unsigned int d2 = 0; d2 < dim; ++d2)
2871 {
2872 for (unsigned int d3 = 0; d3 < dim; ++d3)
2873 {
2874 fe_data.shape_hessians[dof_index3_1][q]
2875 [d1][d2][d3] =
2876 face_xi_grads[m][d1] *
2877 (face_eta_grads[m][d2] *
2878 (polyeta[j][2] *
2879 face_eta_grads[m][d3] *
2880 face_lambda_values[m][q] +
2881 polyeta[j][1] *
2882 face_lambda_grads[m][d3]) +
2883 face_lambda_grads[m][d2] *
2884 polyeta[j][1] *
2885 face_eta_grads[m][d3]);
2886
2887 fe_data.shape_hessians[dof_index3_2][q]
2888 [d1][d2][d3] =
2889 face_eta_grads[m][d1] *
2890 (face_xi_grads[m][d2] *
2891 (polyxi[j][2] *
2892 face_xi_grads[m][d3] *
2893 face_lambda_values[m][q] +
2894 polyxi[j][1] *
2895 face_lambda_grads[m][d3]) +
2896 face_lambda_grads[m][d2] *
2897 polyxi[j][1] * face_xi_grads[m][d3]);
2898 }
2899 }
2900 }
2901 }
2902 }
2903 }
2904 }
2905 }
2906 }
2907}
2908
2909
2910
2911template <int dim, int spacedim>
2912void
2914 const typename Triangulation<dim, dim>::cell_iterator &cell,
2915 const CellSimilarity::Similarity /*cell_similarity*/,
2916 const Quadrature<dim> &quadrature,
2917 const Mapping<dim, dim> &mapping,
2918 const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
2920 &mapping_data,
2921 const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
2923 &data) const
2924{
2925 // Convert to the correct internal data class for this FE class.
2926 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
2928 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
2929
2930 // Now update the edge-based DoFs, which depend on the cell.
2931 // This will fill in the missing items in the InternalData
2932 // (fe_internal/fe_data) which was not filled in by get_data.
2933 fill_edge_values(cell, quadrature, fe_data);
2934 if (dim == 3 && this->degree > 1)
2935 {
2936 fill_face_values(cell, quadrature, fe_data);
2937 }
2938
2939 const UpdateFlags flags(fe_data.update_each);
2940 const unsigned int n_q_points = quadrature.size();
2941
2942 Assert(!(flags & update_values) ||
2943 fe_data.shape_values.size() == this->n_dofs_per_cell(),
2944 ExcDimensionMismatch(fe_data.shape_values.size(),
2945 this->n_dofs_per_cell()));
2946 Assert(!(flags & update_values) ||
2947 fe_data.shape_values[0].size() == n_q_points,
2948 ExcDimensionMismatch(fe_data.shape_values[0].size(), n_q_points));
2949
2950 if (flags & update_values)
2951 {
2952 // Now have all shape_values stored on the reference cell.
2953 // Must now transform to the physical cell.
2954 std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
2955 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2956 {
2957 const unsigned int first =
2958 data.shape_function_to_row_table[dof * this->n_components() +
2959 this->get_nonzero_components(dof)
2960 .first_selected_component()];
2961
2962 mapping.transform(make_array_view(fe_data.shape_values[dof]),
2964 mapping_internal,
2965 make_array_view(transformed_shape_values));
2966 for (unsigned int q = 0; q < n_q_points; ++q)
2967 {
2968 for (unsigned int d = 0; d < dim; ++d)
2969 {
2970 data.shape_values(first + d, q) =
2971 transformed_shape_values[q][d];
2972 }
2973 }
2974 }
2975 }
2976
2977 if (flags & update_gradients)
2978 {
2979 // Now have all shape_grads stored on the reference cell.
2980 // Must now transform to the physical cell.
2981 std::vector<Tensor<2, dim>> input(n_q_points);
2982 std::vector<Tensor<2, dim>> transformed_shape_grads(n_q_points);
2983 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2984 {
2985 for (unsigned int q = 0; q < n_q_points; ++q)
2986 {
2987 input[q] = fe_data.shape_grads[dof][q];
2988 }
2989 mapping.transform(make_array_view(input),
2991 mapping_internal,
2992 make_array_view(transformed_shape_grads));
2993
2994 const unsigned int first =
2995 data.shape_function_to_row_table[dof * this->n_components() +
2996 this->get_nonzero_components(dof)
2997 .first_selected_component()];
2998
2999 for (unsigned int q = 0; q < n_q_points; ++q)
3000 {
3001 for (unsigned int d1 = 0; d1 < dim; ++d1)
3002 {
3003 for (unsigned int d2 = 0; d2 < dim; ++d2)
3004 {
3005 transformed_shape_grads[q][d1] -=
3006 data.shape_values(first + d2, q) *
3007 mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
3008 }
3009 }
3010 }
3011
3012 for (unsigned int q = 0; q < n_q_points; ++q)
3013 {
3014 for (unsigned int d = 0; d < dim; ++d)
3015 {
3016 data.shape_gradients[first + d][q] =
3017 transformed_shape_grads[q][d];
3018 }
3019 }
3020 }
3021 }
3022
3023 if (flags & update_hessians)
3024 {
3025 // Now have all shape_grads stored on the reference cell.
3026 // Must now transform to the physical cell.
3027 std::vector<Tensor<3, dim>> input(n_q_points);
3028 std::vector<Tensor<3, dim>> transformed_shape_hessians(n_q_points);
3029 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
3030 {
3031 for (unsigned int q = 0; q < n_q_points; ++q)
3032 {
3033 input[q] = fe_data.shape_hessians[dof][q];
3034 }
3035 mapping.transform(make_array_view(input),
3037 mapping_internal,
3038 make_array_view(transformed_shape_hessians));
3039
3040 const unsigned int first =
3041 data.shape_function_to_row_table[dof * this->n_components() +
3042 this->get_nonzero_components(dof)
3043 .first_selected_component()];
3044
3045 for (unsigned int q = 0; q < n_q_points; ++q)
3046 {
3047 for (unsigned int d1 = 0; d1 < dim; ++d1)
3048 {
3049 for (unsigned int d2 = 0; d2 < dim; ++d2)
3050 {
3051 for (unsigned int d3 = 0; d3 < dim; ++d3)
3052 {
3053 for (unsigned int d4 = 0; d4 < dim; ++d4)
3054 {
3055 transformed_shape_hessians[q][d1][d3][d4] -=
3056 (data.shape_values(first + d2, q) *
3057 mapping_data
3059 [q][d2][d1][d3][d4]) +
3060 (data.shape_gradients[first + d1][q][d2] *
3061 mapping_data
3063 [d4]) +
3064 (data.shape_gradients[first + d2][q][d3] *
3065 mapping_data
3067 [d4]) +
3068 (data.shape_gradients[first + d2][q][d4] *
3069 mapping_data
3071 [d1]);
3072 }
3073 }
3074 }
3075 }
3076 }
3077
3078 for (unsigned int q = 0; q < n_q_points; ++q)
3079 {
3080 for (unsigned int d = 0; d < dim; ++d)
3081 {
3082 data.shape_hessians[first + d][q] =
3083 transformed_shape_hessians[q][d];
3084 }
3085 }
3086 }
3087 }
3088}
3089
3090
3091
3092template <int dim, int spacedim>
3093void
3095 const typename Triangulation<dim, dim>::cell_iterator &cell,
3096 const unsigned int face_no,
3097 const hp::QCollection<dim - 1> &quadrature,
3098 const Mapping<dim, dim> &mapping,
3099 const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
3101 &mapping_data,
3102 const typename FiniteElement<dim, dim>::InternalDataBase &fe_internal,
3104 &data) const
3105{
3106 AssertDimension(quadrature.size(), 1);
3107
3108 // Note for future improvement:
3109 // We don't have the full quadrature - should use QProjector to create the 2d
3110 // quadrature.
3111 //
3112 // For now I am effectively generating all of the shape function vals/grads,
3113 // etc. On all quad points on all faces and then only using them for one face.
3114 // This is obviously inefficient. I should cache the cell number and cache
3115 // all of the shape_values/gradients etc and then reuse them for each face.
3116
3117 // convert data object to internal
3118 // data for this class. fails with
3119 // an exception if that is not
3120 // possible
3121 Assert(dynamic_cast<const InternalData *>(&fe_internal) != nullptr,
3123 const InternalData &fe_data = static_cast<const InternalData &>(fe_internal);
3124
3125 // Now update the edge-based DoFs, which depend on the cell.
3126 // This will fill in the missing items in the InternalData
3127 // (fe_internal/fe_data) which was not filled in by get_data.
3128 fill_edge_values(cell,
3129 QProjector<dim>::project_to_all_faces(this->reference_cell(),
3130 quadrature[0]),
3131 fe_data);
3132 if (dim == 3 && this->degree > 1)
3133 {
3134 fill_face_values(cell,
3136 this->reference_cell(), quadrature[0]),
3137 fe_data);
3138 }
3139
3140 const UpdateFlags flags(fe_data.update_each);
3141 const unsigned int n_q_points = quadrature[0].size();
3142 const auto offset =
3143 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
3144 face_no,
3145 cell->face_orientation(face_no),
3146 cell->face_flip(face_no),
3147 cell->face_rotation(face_no),
3148 n_q_points);
3149
3150 if (flags & update_values)
3151 {
3152 // Now have all shape_values stored on the reference cell.
3153 // Must now transform to the physical cell.
3154 std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
3155 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
3156 {
3157 mapping.transform(make_array_view(fe_data.shape_values[dof],
3158 offset,
3159 n_q_points),
3161 mapping_internal,
3162 make_array_view(transformed_shape_values));
3163
3164 const unsigned int first =
3165 data.shape_function_to_row_table[dof * this->n_components() +
3166 this->get_nonzero_components(dof)
3167 .first_selected_component()];
3168
3169 for (unsigned int q = 0; q < n_q_points; ++q)
3170 {
3171 for (unsigned int d = 0; d < dim; ++d)
3172 {
3173 data.shape_values(first + d, q) =
3174 transformed_shape_values[q][d];
3175 }
3176 }
3177 }
3178 }
3179 if (flags & update_gradients)
3180 {
3181 // Now have all shape_grads stored on the reference cell.
3182 // Must now transform to the physical cell.
3183 std::vector<Tensor<2, dim>> input(n_q_points);
3184 std::vector<Tensor<2, dim>> transformed_shape_grads(n_q_points);
3185 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
3186 {
3187 for (unsigned int q = 0; q < n_q_points; ++q)
3188 {
3189 input[q] = fe_data.shape_grads[dof][offset + q];
3190 }
3191 mapping.transform(input,
3193 mapping_internal,
3194 make_array_view(transformed_shape_grads));
3195
3196 const unsigned int first =
3197 data.shape_function_to_row_table[dof * this->n_components() +
3198 this->get_nonzero_components(dof)
3199 .first_selected_component()];
3200
3201 for (unsigned int q = 0; q < n_q_points; ++q)
3202 {
3203 for (unsigned int d1 = 0; d1 < dim; ++d1)
3204 {
3205 for (unsigned int d2 = 0; d2 < dim; ++d2)
3206 {
3207 transformed_shape_grads[q][d1] -=
3208 data.shape_values(first + d2, q) *
3209 mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
3210 }
3211 }
3212 }
3213
3214 for (unsigned int q = 0; q < n_q_points; ++q)
3215 {
3216 for (unsigned int d = 0; d < dim; ++d)
3217 {
3218 data.shape_gradients[first + d][q] =
3219 transformed_shape_grads[q][d];
3220 }
3221 }
3222 }
3223 }
3224 if (flags & update_hessians)
3225 {
3226 // Now have all shape_grads stored on the reference cell.
3227 // Must now transform to the physical cell.
3228 std::vector<Tensor<3, dim>> input(n_q_points);
3229 std::vector<Tensor<3, dim>> transformed_shape_hessians(n_q_points);
3230 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
3231 {
3232 for (unsigned int q = 0; q < n_q_points; ++q)
3233 input[q] = fe_data.shape_hessians[dof][offset + q];
3234
3235 mapping.transform(input,
3237 mapping_internal,
3238 make_array_view(transformed_shape_hessians));
3239
3240 const unsigned int first =
3241 data.shape_function_to_row_table[dof * this->n_components() +
3242 this->get_nonzero_components(dof)
3243 .first_selected_component()];
3244
3245 for (unsigned int q = 0; q < n_q_points; ++q)
3246 {
3247 for (unsigned int d1 = 0; d1 < dim; ++d1)
3248 {
3249 for (unsigned int d2 = 0; d2 < dim; ++d2)
3250 {
3251 for (unsigned int d3 = 0; d3 < dim; ++d3)
3252 {
3253 for (unsigned int d4 = 0; d4 < dim; ++d4)
3254 {
3255 transformed_shape_hessians[q][d1][d3][d4] -=
3256 (data.shape_values(first + d2, q) *
3257 mapping_data
3259 [q][d2][d1][d3][d4]) +
3260 (data.shape_gradients[first + d1][q][d2] *
3261 mapping_data
3263 [d4]) +
3264 (data.shape_gradients[first + d2][q][d3] *
3265 mapping_data
3267 [d4]) +
3268 (data.shape_gradients[first + d2][q][d4] *
3269 mapping_data
3271 [d1]);
3272 }
3273 }
3274 }
3275 }
3276 }
3277
3278 for (unsigned int q = 0; q < n_q_points; ++q)
3279 {
3280 for (unsigned int d = 0; d < dim; ++d)
3281 {
3282 data.shape_hessians[first + d][q] =
3283 transformed_shape_hessians[q][d];
3284 }
3285 }
3286 }
3287 }
3288}
3289
3290
3291
3292template <int dim, int spacedim>
3293void
3295 const typename Triangulation<dim, dim>::cell_iterator & /*cell*/,
3296 const unsigned int /*face_no*/,
3297 const unsigned int /*sub_no*/,
3298 const Quadrature<dim - 1> & /*quadrature*/,
3299 const Mapping<dim, dim> & /*mapping*/,
3300 const typename Mapping<dim, dim>::InternalDataBase & /*mapping_internal*/,
3302 & /*mapping_data*/,
3303 const typename FiniteElement<dim, dim>::InternalDataBase & /*fe_internal*/,
3305 & /*data*/) const
3306{
3308}
3309
3310
3311
3312template <int dim, int spacedim>
3335
3336
3337
3338template <int dim, int spacedim>
3339std::string
3341{
3342 // note that the FETools::get_fe_by_name function depends on the particular
3343 // format of the string this function returns, so they have to be kept in sync
3344 std::ostringstream namebuf;
3345 namebuf << "FE_NedelecSZ<" << dim << ">(" << this->degree - 1 << ")";
3346
3347 return namebuf.str();
3348}
3349
3350
3351
3352template <int dim, int spacedim>
3353std::unique_ptr<FiniteElement<dim, dim>>
3355{
3356 return std::make_unique<FE_NedelecSZ<dim, spacedim>>(*this);
3357}
3358
3359
3360
3361template <int dim, int spacedim>
3362std::vector<unsigned int>
3364{
3365 // internal function to return a vector of "dofs per object"
3366 // where the objects inside the vector refer to:
3367 // 0 = vertex
3368 // 1 = edge
3369 // 2 = face (which is a cell in 2d)
3370 // 3 = cell
3371 std::vector<unsigned int> dpo;
3372
3373 dpo.push_back(0);
3374 dpo.push_back(degree + 1);
3375 if (dim > 1)
3376 dpo.push_back(2 * degree * (degree + 1));
3377 if (dim > 2)
3378 dpo.push_back(3 * degree * degree * (degree + 1));
3379
3380 return dpo;
3381}
3382
3383
3384
3385template <int dim, int spacedim>
3386unsigned int
3388{
3389 // Internal function to compute the number of DoFs
3390 // for a given dimension & polynomial order.
3391 switch (dim)
3392 {
3393 case 2:
3394 return 2 * (degree + 1) * (degree + 2);
3395
3396 case 3:
3397 return 3 * (degree + 1) * (degree + 2) * (degree + 2);
3398
3399 default:
3400 {
3402 return 0;
3403 }
3404 }
3405}
3406
3407
3408
3409template <int dim, int spacedim>
3410void
3412{
3413 // fill the 1d polynomials vector:
3414 IntegratedLegendrePolynomials =
3416}
3417
3418
3419
3420// explicit instantiations
3421#include "fe_nedelec_sz.inst"
3422
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:949
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
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_components() const
unsigned int n_unique_faces() const
std::vector< std::pair< std::pair< unsigned int, unsigned int >, unsigned int > > component_to_base_table
Definition fe.h:2558
FullMatrix< double > interface_constraints
Definition fe.h:2423
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
Abstract base class for mapping classes.
Definition mapping.h:316
Definition point.h:111
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
unsigned int size() const
Definition collection.h:264
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > first
Definition grid_out.cc:4613
#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:98
@ mapping_covariant
Definition mapping.h:87
@ mapping_nedelec
Definition mapping.h:127
@ mapping_covariant_hessian
Definition mapping.h:148
#define DEAL_II_NOT_IMPLEMENTED()
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
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)
static unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)