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