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