Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3087-ga35b476a78 2025-04-19 20:30:01+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
fe_nedelec.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2013 - 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
20
22
25#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/mapping.h>
28
29#include <deal.II/grid/tria.h>
31
33#include <deal.II/lac/vector.h>
34
35#include <iostream>
36#include <memory>
37#include <sstream>
38
40
41// #define DEBUG_NEDELEC
42
43namespace internal
44{
45 namespace FE_Nedelec
46 {
47 namespace
48 {
49 double
50 get_embedding_computation_tolerance(const unsigned int p)
51 {
52 // This heuristic was computed by monitoring the worst residual
53 // resulting from the least squares computation when computing
54 // the face embedding matrices in the FE_Nedelec constructor.
55 // The residual growth is exponential, but is bounded by this
56 // function up to degree 12.
57 return 1.e-15 * std::exp(std::pow(p, 1.075));
58 }
59 } // namespace
60 } // namespace FE_Nedelec
61} // namespace internal
62
63
64// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
65// adjust_line_dof_index_for_line_orientation_table fields, and write tests
66// similar to bits/face_orientation_and_fe_q_*
67
68template <int dim>
69FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
70 : FE_PolyTensor<dim>(
71 PolynomialsNedelec<dim>(order),
72 FiniteElementData<dim>(get_dpo_vector(order),
73 dim,
74 order + 1,
75 FiniteElementData<dim>::Hcurl),
76 std::vector<bool>(PolynomialsNedelec<dim>::n_polynomials(order), true),
77 std::vector<ComponentMask>(PolynomialsNedelec<dim>::n_polynomials(order),
78 ComponentMask(std::vector<bool>(dim, true))))
79{
80#ifdef DEBUG_NEDELEC
81 deallog << get_name() << std::endl;
82#endif
83
84 Assert(dim >= 2, ExcImpossibleInDim(dim));
85
87 // First, initialize the
88 // generalized support points and
89 // quadrature weights, since they
90 // are required for interpolation.
92
93 // We already use the correct basis, so no basis transformation is required
94 // from the polynomial space we have described above to the one that is dual
95 // to the node functionals. As documented in the base class, this is
96 // expressed by setting the inverse node matrix to the empty matrix.
97 this->inverse_node_matrix.clear();
98
99 // do not initialize embedding and restriction here. these matrices are
100 // initialized on demand in get_restriction_matrix and
101 // get_prolongation_matrix
102
103#ifdef DEBUG_NEDELEC
104 deallog << "Face Embedding" << std::endl;
105#endif
107
108 // TODO: the implementation makes the assumption that all faces have the
109 // same number of dofs
111 const unsigned int face_no = 0;
112
113 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
114 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
115 this->n_dofs_per_face(face_no));
116
117 FETools::compute_face_embedding_matrices<dim, double>(
118 *this,
119 face_embeddings,
120 0,
121 0,
122 internal::FE_Nedelec::get_embedding_computation_tolerance(order));
123
124 switch (dim)
125 {
126 case 1:
127 {
128 this->interface_constraints.reinit(0, 0);
129 break;
130 }
131
132 case 2:
133 {
134 this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no),
135 this->n_dofs_per_face(face_no));
136
137 for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
138 ++i)
139 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
140 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
141 this->interface_constraints(i * this->n_dofs_per_face(face_no) +
142 j,
143 k) = face_embeddings[i](j, k);
144
145 break;
146 }
147
148 case 3:
149 {
150 this->interface_constraints.reinit(
151 4 * (this->n_dofs_per_face(face_no) - this->degree),
152 this->n_dofs_per_face(face_no));
153
154 unsigned int target_row = 0;
155
156 for (unsigned int i = 0; i < 2; ++i)
157 for (unsigned int j = this->degree; j < 2 * this->degree;
158 ++j, ++target_row)
159 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
160 this->interface_constraints(target_row, k) =
161 face_embeddings[2 * i](j, k);
162
163 for (unsigned int i = 0; i < 2; ++i)
164 for (unsigned int j = 3 * this->degree;
165 j < GeometryInfo<3>::lines_per_face * this->degree;
166 ++j, ++target_row)
167 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
168 this->interface_constraints(target_row, k) =
169 face_embeddings[i](j, k);
170
171 for (unsigned int i = 0; i < 2; ++i)
172 for (unsigned int j = 0; j < 2; ++j)
173 for (unsigned int k = i * this->degree;
174 k < (i + 1) * this->degree;
175 ++k, ++target_row)
176 for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
177 ++l)
178 this->interface_constraints(target_row, l) =
179 face_embeddings[i + 2 * j](k, l);
180
181 for (unsigned int i = 0; i < 2; ++i)
182 for (unsigned int j = 0; j < 2; ++j)
183 for (unsigned int k = (i + 2) * this->degree;
184 k < (i + 3) * this->degree;
185 ++k, ++target_row)
186 for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
187 ++l)
188 this->interface_constraints(target_row, l) =
189 face_embeddings[2 * i + j](k, l);
190
191 for (unsigned int i = 0; i < GeometryInfo<3>::max_children_per_face;
192 ++i)
193 for (unsigned int j =
194 GeometryInfo<3>::lines_per_face * this->degree;
195 j < this->n_dofs_per_face(face_no);
196 ++j, ++target_row)
197 for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
198 this->interface_constraints(target_row, k) =
199 face_embeddings[i](j, k);
200
201 break;
202 }
203
204 default:
206 }
207
208 // We need to initialize the dof permutation table and the one for the sign
209 // change.
211}
212
213
214template <int dim>
215void
217{
218 // The order of the Nedelec elements equals the tensor degree minus one,
219 // k = n - 1. In the three-dimensional space the Nedelec elements of the
220 // lowermost order, k = 0, have only 12 line (edge) dofs. The Nedelec
221 // elements of the higher orders, k > 0, have 3*(k+1)*(k+2)^2 dofs in
222 // total if dim=3. The dofs in a cell are distributed between lines
223 // (edges), quads (faces), and the hex (the interior of the cell) as the
224 // following:
225 //
226 // 12*(k+1) line dofs; (k+1) dofs per line.
227 // 2*6*k*(k+1) quad dofs; 2*k*(k+1) dofs per quad.
228 // 3*(k+2)^2*(k+1) hex dofs;
229 //
230 // The dofs are indexed in the following order: first all line dofs,
231 // then all quad dofs, and then all hex dofs.
232 //
233 // The line dofs need only sign adjustments. No permutation of line
234 // dofs is needed. The line dofs are treated by
235 // internal::FE_PolyTensor::get_dof_sign_change_nedelec(...)
236 // in fe_poly_tensor.cc.
237 //
238 // The hex dofs need no adjustments: they are not shared between
239 // neighbouring mesh cells.
240 //
241 // The two-dimensional Nedelec finite elements share no quad dofs between
242 // neighbouring mesh cells. The zero-order three-dimensional Nedelec
243 // finite elements have no quad dofs. Consequently, here we treat only
244 // quad dofs of the three-dimensional Nedelec finite elements of the
245 // higher orders, k>0. The questions how the curl looks like in the
246 // higher-dimensional spaces and what does it mean to be curl-conforming
247 // if dim>3 we leave unanswered.
248 //
249 // In this function we need to change some entries in the following two
250 // vectors of tables:
251 // adjust_quad_dof_index_for_face_orientation_table
252 // and
253 // adjust_quad_dof_sign_for_face_orientation_table.
254 // These tables specify the permutations and sign adjustments of the quad
255 // dofs only. The tables are already filled with zeros meaning no
256 // permutations or sign change are required. We need to change some
257 // entries of the tables such that the shape functions that correspond to
258 // the quad dofs and are shared between neighbouring cells have consistent
259 // orientations.
260 //
261 // The swap tables below describe the dof permutations and sign changes
262 // that need to be done. The function
263 // FE_Nedelec<dim>::initialize_quad_dof_index_permutation_and_sign_change()
264 // simply reads the information in the swap tables below and puts it into
265 // tables
266 // adjust_quad_dof_index_for_face_orientation_table
267 // and
268 // adjust_quad_dof_sign_for_face_orientation_table.
269 // A good question is: why don't we put the information into the tables of
270 // deal.II right away? The answer is the following. The information on the
271 // necessary dof permutations and sign changes is derived by plotting the
272 // shape functions and observing them on faces of different orientations.
273 // It is convenient to put the observations first in the format of the
274 // swap tables below and then convert the swap tables into the format used
275 // by deal.II.
276 //
277 // The dofs on a quad are indexed as the following:
278 //
279 // | x0, x1, x2, x3, ..., xk | y0, y1, y2, y3 ..., yk |
280 // | | |
281 // |<------ k*(k+1) --------->|<------ k*(k+1) -------->|
282 // | |
283 // |<------------------- 2*k*(k+1) -------------------->|
284 //
285 // Only one type of dof permutation is needed: swap between two dofs; one
286 // dof being xi, another yj. That is, if x4 is replaced with y7,
287 // then y7 must be replaced with x4. Such swaps can be ordered as
288 // illustrated by the following example:
289 //
290 // *
291 // y0, y9, y1, y2, ..., yk
292 // --------------------------- (swap)
293 // x0, x1, x2, x3, ..., xk
294 // *
295 //
296 // An x-dof below the line is swapped with the corresponding y-dof above
297 // the line. A dof marked by the asterisk must change its sign before the
298 // swap.
299 //
300 // The x-dofs are assumed to have the normal order. There is no need to
301 // encode it. Therefore, the swap tables need to encode the following
302 // information: indices of the y-dofs, the sign change of the x-dofs, and
303 // sign change of the y-dofs. The swap above is encoded as the following:
304 //
305 // swap = { 0, 9, 1, 2, ...., yk, // indices of the y-dofs
306 // 1, 0, 0, 0, ...., 0, // sign change of the x-dofs,
307 // 0, 1, 0, 0, ...., 0}; // sign change of the y-dofs.
308 //
309 // If no swap is needed, -1 is placed instead of the y-dof index.
310 //
311 // Such swaps are assembled into the swap table:
312 //
313 // swap_table = {swap_0, swap_1, ... swap_7};
314 //
315 // Each swap table contains eight swaps - one swap for each possible quad
316 // orientation. The deal.II encodes the orientation of a quad using
317 // three boolean parameters:
318 // face_orientation - true if face is in standard orientation
319 // and false otherwise;
320 // face_rotation - rotation by 90 deg counterclockwise if true;
321 // face_flip - rotation by 180 deg counterclockwise if true.
322 // See the documentation of GeometryInfo<dim>.
323 //
324 // The combined face orientation is computes as
325 // orientation_no = face_flip*4 + face_rotation*2 + face_orientation*1;
326 // See tria_orientation.h.
327 //
328 // The parameter orientation_no (0...7) indexes the swaps in a swap table.
329 //
330 // Nedelec elements of order k have their own swap table, swap_table_k.
331 // Recall, the swap_table_0 is empty as the Nedelec finite elements of the
332 // lowermost order have no quad dofs.
333
334 static const int c_swap_table_0 = 0;
335
336 static const int c_swap_table_1[8][3][2] = { // 0 1
337 {{-1, -1}, // 0
338 {0, 0},
339 {0, 0}},
340 {{0, 1}, // 1
341 {0, 0},
342 {0, 0}},
343 {{0, 1}, // 2
344 {1, 0},
345 {0, 0}},
346 {{-1, -1}, // 3
347 {0, 0},
348 {1, 0}},
349 {{-1, -1}, // 4
350 {1, 0},
351 {1, 0}},
352 {{0, 1}, // 5
353 {1, 0},
354 {1, 0}},
355 {{0, 1}, // 6
356 {0, 0},
357 {1, 0}},
358 {{-1, -1}, // 7
359 {1, 0},
360 {0, 0}}};
361
362 static const int c_swap_table_2[8][3][6] = {// 0 1 2 3 4 5
363 {{-1, -1, -1, -1, -1, -1}, // 0
364 {0, 0, 0, 0, 0, 0},
365 {0, 0, 0, 0, 0, 0}},
366 {{0, 3, 1, 4, 2, 5}, // 1
367 {0, 0, 0, 0, 0, 0},
368 {0, 0, 0, 0, 0, 0}},
369 {{0, 3, 1, 4, 2, 5}, // 2
370 {1, 1, 0, 0, 1, 1},
371 {0, 0, 0, 1, 1, 1}},
372 {{-1, -1, -1, -1, -1, -1}, // 3
373 {0, 1, 0, 1, 0, 1},
374 {1, 0, 1, 1, 0, 1}},
375 {{-1, -1, -1, -1, -1, -1}, // 4
376 {1, 0, 0, 1, 1, 0},
377 {1, 0, 1, 0, 1, 0}},
378 {{0, 3, 1, 4, 2, 5}, // 5
379 {1, 0, 0, 1, 1, 0},
380 {1, 0, 1, 0, 1, 0}},
381 {{0, 3, 1, 4, 2, 5}, // 6
382 {0, 1, 0, 1, 0, 1},
383 {1, 0, 1, 1, 0, 1}},
384 {{-1, -1, -1, -1, -1, -1}, // 7
385 {1, 1, 0, 0, 1, 1},
386 {0, 0, 0, 1, 1, 1}}};
387
388 static const int c_swap_table_3[8][3][12] = {
389 // 0 1 2 3 4 5 6 7 8 9 10 11
390 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 0
391 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
392 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
393 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 1
394 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
395 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
396 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 2
397 {1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0},
398 {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}},
399 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 3
400 {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
401 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}},
402 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 4
403 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
404 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
405 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 5
406 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
407 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
408 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11}, // 6
409 {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
410 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}},
411 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 7
412 {1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0},
413 {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}}};
414
415 static const int c_swap_table_4[8][3][20] = {
416 // Swap sign_X and sign_Y rows if k=4. Why?...
417 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
418 // 19
419 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
420 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 0
421 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
422 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
423 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
424 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 1
425 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
426 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
427 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
428 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 2
429 {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
430 {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}},
431 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
432 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 3
433 {1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1},
434 {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}},
435 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
436 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 4
437 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
438 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
439 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
440 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 5
441 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
442 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
443 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
444 12, 17, 3, 8, 13, 18, 4, 9, 14, 19}, // 6
445 {1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1},
446 {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}},
447 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
448 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // 7
449 {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
450 {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}}};
451
452 static const int *swap_table_array[5] = {&c_swap_table_0,
453 &c_swap_table_1[0][0][0],
454 &c_swap_table_2[0][0][0],
455 &c_swap_table_3[0][0][0],
456 &c_swap_table_4[0][0][0]};
457
458 static const int row_length[5] = {0, 2, 6, 12, 20};
459 static const int table_size[5] = {
460 0, 8 * 3 * 2, 8 * 3 * 6, 8 * 3 * 12, 8 * 3 * 20};
461
462 // Only three-dimensional Nedelec finite elements are treated. The
463 // two-dimensional Nedelec finite elements only need sign adjustments of the
464 // line dofs. These adjustments are done by
465 // internal::FE_PolyTensor::get_dof_sign_change_nedelec(...)
466 // in fe_poly_tensor.cc. The notions of curl and curl-conforming finite
467 // elements in higher-dimensional spaces, dim >3, are somewhat unclear as
468 // curl, strictly peaking, exists only in the three-dimensional space.
469 if (dim != 3)
470 return;
471
472 const unsigned int k = this->tensor_degree() - 1;
473
474 // The Nedelec finite elements of the lowermost order have no quad dofs.
475 if (k == 0)
476 return;
477
478 // The finite element orders > 4 are not implemented.
480
481 // TODO: the implementation makes the assumption that all quads have the
482 // same number of dofs
483 AssertDimension(this->n_unique_faces(), 1);
484 const unsigned int face_no = 0;
485
486 Assert(
487 this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
488 this->reference_cell().n_face_orientations(face_no) *
489 this->n_dofs_per_quad(face_no),
491
492 Assert(
493 this->adjust_quad_dof_sign_for_face_orientation_table[0].n_elements() ==
494 this->reference_cell().n_face_orientations(face_no) *
495 this->n_dofs_per_quad(face_no),
497
498 // The 3D Nedelec finite elements have 2*k*(k+1) dofs per each quad.
499 Assert(2 * k * (k + 1) == this->n_dofs_per_quad(face_no), ExcInternalError());
500
501 const int *swap_table = swap_table_array[k];
502
503 const unsigned int half_dofs = k * (k + 1); // see below;
504
505 const int rl = row_length[k];
506
507 int offset = table_size[0];
508 // The assignment above is to prevent the compiler from complaining about the
509 // unused table_size.
510
511 int value = 0;
512
513 for (const bool face_orientation : {false, true})
514 for (const bool face_rotation : {false, true})
515 for (const bool face_flip : {false, true})
516 {
517 const auto case_no =
519 face_rotation,
520 face_flip);
521
522 // The dofs on a quad are indexed as the following:
523 //
524 // | x0, x1, x2, x3, ..., xk | y0, y1, y2, y3 ..., yk |
525 // | | |
526 // |-- half_ dofs = k*(k+1) --|-- half_dofs = k*(k+1) --|
527 // | |
528 // |-------------------- 2*k*(k+1) ---------------------|
529
530 for (unsigned int indx_x = 0; indx_x < half_dofs; indx_x++)
531 {
532 offset = 3 * rl * case_no + 0 * rl + indx_x;
533 Assert(offset < table_size[k], ExcInternalError());
534 value = *(swap_table + offset);
535 Assert(value < table_size[k], ExcInternalError());
536 Assert(value > -2, ExcInternalError());
537
538 if (value != -1)
539 {
540 const unsigned int indx_y =
541 half_dofs + static_cast<unsigned int>(value);
542
543 // dofs swap
544 this
545 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
546 indx_x, case_no) = indx_y - indx_x;
547
548 this
549 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
550 indx_y, case_no) = indx_x - indx_y;
551 }
552
553 // dof sign change
554 offset = 3 * rl * case_no + 1 * rl + indx_x;
555 Assert(offset < table_size[k], ExcInternalError());
556 value = *(swap_table + offset);
557 Assert((value == 0) || (value == 1), ExcInternalError());
558
559 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
560 indx_x, case_no) = static_cast<bool>(value);
561
562
563 offset = 3 * rl * case_no + 2 * rl + indx_x;
564 Assert(offset < table_size[k], ExcInternalError());
565 value = *(swap_table + offset);
566 Assert((value == 0) || (value == 1), ExcInternalError());
567
568 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
569 indx_x + half_dofs, case_no) = static_cast<bool>(value);
570 }
571 }
572
573 return;
574}
575
576
577template <int dim>
578std::string
580{
581 // note that the
582 // FETools::get_fe_by_name
583 // function depends on the
584 // particular format of the string
585 // this function returns, so they
586 // have to be kept in synch
587
588 std::ostringstream namebuf;
589 namebuf << "FE_Nedelec<" << dim << ">(" << this->degree - 1 << ")";
590
591 return namebuf.str();
592}
593
594
595template <int dim>
596std::unique_ptr<FiniteElement<dim, dim>>
598{
599 return std::make_unique<FE_Nedelec<dim>>(*this);
600}
601
602//---------------------------------------------------------------------------
603// Auxiliary and internal functions
604//---------------------------------------------------------------------------
605
606
607
608// Set the generalized support
609// points and precompute the
610// parts of the projection-based
611// interpolation, which does
612// not depend on the interpolated
613// function.
614template <>
615void
620
621
622
623template <>
624void
626{
627 const int dim = 2;
628
629 // TODO: the implementation makes the assumption that all faces have the
630 // same number of dofs
631 AssertDimension(this->n_unique_faces(), 1);
632 const unsigned int face_no = 0;
633
634 // Create polynomial basis.
635 const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
637 std::vector<Polynomials::Polynomial<double>> lobatto_polynomials_grad(order +
638 1);
639
640 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
641 lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative();
642
643 // Initialize quadratures to obtain
644 // quadrature points later on.
645 const QGauss<dim - 1> reference_edge_quadrature(order + 1);
646 const unsigned int n_edge_points = reference_edge_quadrature.size();
647 const unsigned int n_boundary_points =
648 GeometryInfo<dim>::lines_per_cell * n_edge_points;
649 const Quadrature<dim> edge_quadrature =
650 QProjector<dim>::project_to_all_faces(this->reference_cell(),
651 reference_edge_quadrature);
652
653 this->generalized_face_support_points[face_no].resize(n_edge_points);
654
655 // Create face support points.
656 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
657 this->generalized_face_support_points[face_no][q_point] =
658 reference_edge_quadrature.point(q_point);
659
660 if (order > 0)
661 {
662 // If the polynomial degree is positive
663 // we have support points on the faces
664 // and in the interior of a cell.
665 const QGauss<dim> quadrature(order + 1);
666 const unsigned int n_interior_points = quadrature.size();
667
668 this->generalized_support_points.resize(n_boundary_points +
669 n_interior_points);
670 boundary_weights.reinit(n_edge_points, order);
671
672 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
673 {
674 for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
675 ++line)
676 this->generalized_support_points[line * n_edge_points + q_point] =
678 this->reference_cell(),
679 line,
681 n_edge_points) +
682 q_point);
683
684 for (unsigned int i = 0; i < order; ++i)
685 boundary_weights(q_point, i) =
686 reference_edge_quadrature.weight(q_point) *
687 lobatto_polynomials_grad[i + 1].value(
688 this->generalized_face_support_points[face_no][q_point][0]);
689 }
690
691 for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
692 this->generalized_support_points[q_point + n_boundary_points] =
693 quadrature.point(q_point);
694 }
695
696 else
697 {
698 // In this case we only need support points
699 // on the faces of a cell.
700 this->generalized_support_points.resize(n_boundary_points);
701
702 for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
703 ++line)
704 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
705 this->generalized_support_points[line * n_edge_points + q_point] =
707 this->reference_cell(),
708 line,
710 n_edge_points) +
711 q_point);
712 }
713}
714
715
716
717template <>
718void
720{
721 const int dim = 3;
722
723 // TODO: the implementation makes the assumption that all faces have the
724 // same number of dofs
725 AssertDimension(this->n_unique_faces(), 1);
726 const unsigned int face_no = 0;
727
728 // Create polynomial basis.
729 const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
731 std::vector<Polynomials::Polynomial<double>> lobatto_polynomials_grad(order +
732 1);
733
734 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
735 lobatto_polynomials_grad[i] = lobatto_polynomials[i + 1].derivative();
736
737 // Initialize quadratures to obtain
738 // quadrature points later on.
739 const QGauss<1> reference_edge_quadrature(order + 1);
740 const unsigned int n_edge_points = reference_edge_quadrature.size();
741 const Quadrature<dim - 1> &edge_quadrature =
743 ReferenceCells::get_hypercube<dim - 1>(), reference_edge_quadrature);
744
745 if (order > 0)
746 {
747 // If the polynomial order is positive
748 // we have support points on the edges,
749 // faces and in the interior of a cell.
750 const QGauss<dim - 1> reference_face_quadrature(order + 1);
751 const unsigned int n_face_points = reference_face_quadrature.size();
752 const unsigned int n_boundary_points =
753 GeometryInfo<dim>::lines_per_cell * n_edge_points +
754 GeometryInfo<dim>::faces_per_cell * n_face_points;
755 const QGauss<dim> quadrature(order + 1);
756 const unsigned int n_interior_points = quadrature.size();
757
758 boundary_weights.reinit(n_edge_points + n_face_points,
759 2 * (order + 1) * order);
760 this->generalized_face_support_points[face_no].resize(4 * n_edge_points +
761 n_face_points);
762 this->generalized_support_points.resize(n_boundary_points +
763 n_interior_points);
764
765 // Create support points on edges.
766 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
767 {
768 for (unsigned int line = 0;
769 line < GeometryInfo<dim - 1>::lines_per_cell;
770 ++line)
771 this
772 ->generalized_face_support_points[face_no][line * n_edge_points +
773 q_point] =
774 edge_quadrature.point(
776 ReferenceCells::get_hypercube<dim - 1>(),
777 line,
779 n_edge_points) +
780 q_point);
781
782 for (unsigned int i = 0; i < 2; ++i)
783 for (unsigned int j = 0; j < 2; ++j)
784 {
785 this->generalized_support_points[q_point +
786 (i + 4 * j) * n_edge_points] =
787 Point<dim>(i, reference_edge_quadrature.point(q_point)[0], j);
788 this->generalized_support_points[q_point + (i + 4 * j + 2) *
789 n_edge_points] =
790 Point<dim>(reference_edge_quadrature.point(q_point)[0], i, j);
791 this->generalized_support_points[q_point + (i + 2 * (j + 4)) *
792 n_edge_points] =
793 Point<dim>(i, j, reference_edge_quadrature.point(q_point)[0]);
794 }
795
796 for (unsigned int i = 0; i < order; ++i)
797 boundary_weights(q_point, i) =
798 reference_edge_quadrature.weight(q_point) *
799 lobatto_polynomials_grad[i + 1].value(
800 this->generalized_face_support_points[face_no][q_point][1]);
801 }
802
803 // Create support points on faces.
804 for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
805 {
806 this->generalized_face_support_points[face_no]
807 [q_point + 4 * n_edge_points] =
808 reference_face_quadrature.point(q_point);
809
810 for (unsigned int i = 0; i <= order; ++i)
811 for (unsigned int j = 0; j < order; ++j)
812 {
813 boundary_weights(q_point + n_edge_points, 2 * (i * order + j)) =
814 reference_face_quadrature.weight(q_point) *
815 lobatto_polynomials_grad[i].value(
816 this->generalized_face_support_points
817 [face_no][q_point + 4 * n_edge_points][0]) *
818 lobatto_polynomials[j + 2].value(
819 this->generalized_face_support_points
820 [face_no][q_point + 4 * n_edge_points][1]);
821 boundary_weights(q_point + n_edge_points,
822 2 * (i * order + j) + 1) =
823 reference_face_quadrature.weight(q_point) *
824 lobatto_polynomials_grad[i].value(
825 this->generalized_face_support_points
826 [face_no][q_point + 4 * n_edge_points][1]) *
827 lobatto_polynomials[j + 2].value(
828 this->generalized_face_support_points
829 [face_no][q_point + 4 * n_edge_points][0]);
830 }
831 }
832
833 const Quadrature<dim> &face_quadrature =
834 QProjector<dim>::project_to_all_faces(this->reference_cell(),
835 reference_face_quadrature);
836
837 for (const unsigned int face : GeometryInfo<dim>::face_indices())
838 for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
839 {
840 this->generalized_support_points[face * n_face_points + q_point +
842 n_edge_points] =
844 this->reference_cell(),
845 face,
847 n_face_points) +
848 q_point);
849 }
850
851 // Create support points in the interior.
852 for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
853 this->generalized_support_points[q_point + n_boundary_points] =
854 quadrature.point(q_point);
855 }
856
857 else
858 {
859 this->generalized_face_support_points[face_no].resize(4 * n_edge_points);
860 this->generalized_support_points.resize(
861 GeometryInfo<dim>::lines_per_cell * n_edge_points);
862
863 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
864 {
865 for (unsigned int line = 0;
866 line < GeometryInfo<dim - 1>::lines_per_cell;
867 ++line)
868 this
869 ->generalized_face_support_points[face_no][line * n_edge_points +
870 q_point] =
871 edge_quadrature.point(
873 ReferenceCells::get_hypercube<dim - 1>(),
874 line,
876 n_edge_points) +
877 q_point);
878
879 for (unsigned int i = 0; i < 2; ++i)
880 for (unsigned int j = 0; j < 2; ++j)
881 {
882 this->generalized_support_points[q_point +
883 (i + 4 * j) * n_edge_points] =
884 Point<dim>(i, reference_edge_quadrature.point(q_point)[0], j);
885 this->generalized_support_points[q_point + (i + 4 * j + 2) *
886 n_edge_points] =
887 Point<dim>(reference_edge_quadrature.point(q_point)[0], i, j);
888 this->generalized_support_points[q_point + (i + 2 * (j + 4)) *
889 n_edge_points] =
890 Point<dim>(i, j, reference_edge_quadrature.point(q_point)[0]);
891 }
892 }
893 }
894}
895
896
897
898// Set the restriction matrices.
899template <>
900void
902{
903 // there is only one refinement case in 1d,
904 // which is the isotropic one
905 for (unsigned int i = 0; i < GeometryInfo<1>::max_children_per_cell; ++i)
906 this->restriction[0][i].reinit(0, 0);
907}
908
909
910
911// Restriction operator
912template <int dim>
913void
915{
916 // This function does the same as the
917 // function interpolate further below.
918 // But since the functions, which we
919 // interpolate here, are discontinuous
920 // we have to use more quadrature
921 // points as in interpolate.
922 const QGauss<1> edge_quadrature(2 * this->degree);
923 const std::vector<Point<1>> &edge_quadrature_points =
924 edge_quadrature.get_points();
925 const unsigned int n_edge_quadrature_points = edge_quadrature.size();
926 const unsigned int index = RefinementCase<dim>::isotropic_refinement - 1;
927
928 switch (dim)
929 {
930 case 2:
931 {
932 // First interpolate the shape
933 // functions of the child cells
934 // to the lowest order shape
935 // functions of the parent cell.
936 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
937 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
938 ++q_point)
939 {
940 const double weight = 2.0 * edge_quadrature.weight(q_point);
941
942 if (edge_quadrature_points[q_point][0] < 0.5)
943 {
944 Point<dim> quadrature_point(
945 0.0, 2.0 * edge_quadrature_points[q_point][0]);
946
947 this->restriction[index][0](0, dof) +=
948 weight *
949 this->shape_value_component(dof, quadrature_point, 1);
950 quadrature_point[0] = 1.0;
951 this->restriction[index][1](this->degree, dof) +=
952 weight *
953 this->shape_value_component(dof, quadrature_point, 1);
954 quadrature_point[0] = quadrature_point[1];
955 quadrature_point[1] = 0.0;
956 this->restriction[index][0](2 * this->degree, dof) +=
957 weight *
958 this->shape_value_component(dof, quadrature_point, 0);
959 quadrature_point[1] = 1.0;
960 this->restriction[index][2](3 * this->degree, dof) +=
961 weight *
962 this->shape_value_component(dof, quadrature_point, 0);
963 }
964
965 else
966 {
967 Point<dim> quadrature_point(
968 0.0, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
969
970 this->restriction[index][2](0, dof) +=
971 weight *
972 this->shape_value_component(dof, quadrature_point, 1);
973 quadrature_point[0] = 1.0;
974 this->restriction[index][3](this->degree, dof) +=
975 weight *
976 this->shape_value_component(dof, quadrature_point, 1);
977 quadrature_point[0] = quadrature_point[1];
978 quadrature_point[1] = 0.0;
979 this->restriction[index][1](2 * this->degree, dof) +=
980 weight *
981 this->shape_value_component(dof, quadrature_point, 0);
982 quadrature_point[1] = 1.0;
983 this->restriction[index][3](3 * this->degree, dof) +=
984 weight *
985 this->shape_value_component(dof, quadrature_point, 0);
986 }
987 }
988
989 // Then project the shape functions
990 // of the child cells to the higher
991 // order shape functions of the
992 // parent cell.
993 if (this->degree > 1)
994 {
995 const unsigned int deg = this->degree - 1;
996 const std::vector<Polynomials::Polynomial<double>>
997 &legendre_polynomials =
999 FullMatrix<double> system_matrix_inv(deg, deg);
1000
1001 {
1002 FullMatrix<double> assembling_matrix(deg,
1003 n_edge_quadrature_points);
1004
1005 for (unsigned int q_point = 0;
1006 q_point < n_edge_quadrature_points;
1007 ++q_point)
1008 {
1009 const double weight =
1010 std::sqrt(edge_quadrature.weight(q_point));
1011
1012 for (unsigned int i = 0; i < deg; ++i)
1013 assembling_matrix(i, q_point) =
1014 weight * legendre_polynomials[i + 1].value(
1015 edge_quadrature_points[q_point][0]);
1016 }
1017
1018 FullMatrix<double> system_matrix(deg, deg);
1019
1020 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1021 system_matrix_inv.invert(system_matrix);
1022 }
1023
1024 FullMatrix<double> solution(this->degree - 1, 4);
1025 FullMatrix<double> system_rhs(this->degree - 1, 4);
1026 Vector<double> tmp(4);
1027
1028 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1029 for (unsigned int i = 0; i < 2; ++i)
1030 {
1031 system_rhs = 0.0;
1032
1033 for (unsigned int q_point = 0;
1034 q_point < n_edge_quadrature_points;
1035 ++q_point)
1036 {
1037 const double weight = edge_quadrature.weight(q_point);
1038 const Point<dim> quadrature_point_0(
1039 i, edge_quadrature_points[q_point][0]);
1040 const Point<dim> quadrature_point_1(
1041 edge_quadrature_points[q_point][0], i);
1042
1043 if (edge_quadrature_points[q_point][0] < 0.5)
1044 {
1045 Point<dim> quadrature_point_2(
1046 i, 2.0 * edge_quadrature_points[q_point][0]);
1047
1048 tmp(0) =
1049 weight *
1050 (2.0 * this->shape_value_component(
1051 dof, quadrature_point_2, 1) -
1052 this->restriction[index][i](i * this->degree,
1053 dof) *
1054 this->shape_value_component(i * this->degree,
1055 quadrature_point_0,
1056 1));
1057 tmp(1) =
1058 -1.0 * weight *
1059 this->restriction[index][i + 2](i * this->degree,
1060 dof) *
1061 this->shape_value_component(i * this->degree,
1062 quadrature_point_0,
1063 1);
1064 quadrature_point_2 = Point<dim>(
1065 2.0 * edge_quadrature_points[q_point][0], i);
1066 tmp(2) =
1067 weight *
1068 (2.0 * this->shape_value_component(
1069 dof, quadrature_point_2, 0) -
1070 this->restriction[index][2 * i]((i + 2) *
1071 this->degree,
1072 dof) *
1073 this->shape_value_component((i + 2) *
1074 this->degree,
1075 quadrature_point_1,
1076 0));
1077 tmp(3) =
1078 -1.0 * weight *
1079 this->restriction[index][2 * i + 1](
1080 (i + 2) * this->degree, dof) *
1081 this->shape_value_component(
1082 (i + 2) * this->degree, quadrature_point_1, 0);
1083 }
1084
1085 else
1086 {
1087 tmp(0) =
1088 -1.0 * weight *
1089 this->restriction[index][i](i * this->degree,
1090 dof) *
1091 this->shape_value_component(i * this->degree,
1092 quadrature_point_0,
1093 1);
1094
1095 Point<dim> quadrature_point_2(
1096 i,
1097 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1098
1099 tmp(1) =
1100 weight *
1101 (2.0 * this->shape_value_component(
1102 dof, quadrature_point_2, 1) -
1103 this->restriction[index][i + 2](i * this->degree,
1104 dof) *
1105 this->shape_value_component(i * this->degree,
1106 quadrature_point_0,
1107 1));
1108 tmp(2) =
1109 -1.0 * weight *
1110 this->restriction[index][2 * i]((i + 2) *
1111 this->degree,
1112 dof) *
1113 this->shape_value_component(
1114 (i + 2) * this->degree, quadrature_point_1, 0);
1115 quadrature_point_2 = Point<dim>(
1116 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1117 i);
1118 tmp(3) =
1119 weight *
1120 (2.0 * this->shape_value_component(
1121 dof, quadrature_point_2, 0) -
1122 this->restriction[index][2 * i + 1](
1123 (i + 2) * this->degree, dof) *
1124 this->shape_value_component((i + 2) *
1125 this->degree,
1126 quadrature_point_1,
1127 0));
1128 }
1129
1130 for (unsigned int j = 0; j < this->degree - 1; ++j)
1131 {
1132 const double L_j =
1133 legendre_polynomials[j + 1].value(
1134 edge_quadrature_points[q_point][0]);
1135
1136 for (unsigned int k = 0; k < tmp.size(); ++k)
1137 system_rhs(j, k) += tmp(k) * L_j;
1138 }
1139 }
1140
1141 system_matrix_inv.mmult(solution, system_rhs);
1142
1143 for (unsigned int j = 0; j < this->degree - 1; ++j)
1144 for (unsigned int k = 0; k < 2; ++k)
1145 {
1146 if (std::abs(solution(j, k)) > 1e-14)
1147 this->restriction[index][i + 2 * k](
1148 i * this->degree + j + 1, dof) = solution(j, k);
1149
1150 if (std::abs(solution(j, k + 2)) > 1e-14)
1151 this->restriction[index][2 * i + k](
1152 (i + 2) * this->degree + j + 1, dof) =
1153 solution(j, k + 2);
1154 }
1155 }
1156
1157 const QGauss<dim> quadrature(2 * this->degree);
1158 const std::vector<Point<dim>> &quadrature_points =
1159 quadrature.get_points();
1160 const std::vector<Polynomials::Polynomial<double>>
1161 &lobatto_polynomials =
1163 const unsigned int n_boundary_dofs =
1164 GeometryInfo<dim>::faces_per_cell * this->degree;
1165 const unsigned int n_quadrature_points = quadrature.size();
1166
1167 {
1168 FullMatrix<double> assembling_matrix((this->degree - 1) *
1169 this->degree,
1170 n_quadrature_points);
1171
1172 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1173 ++q_point)
1174 {
1175 const double weight = std::sqrt(quadrature.weight(q_point));
1176
1177 for (unsigned int i = 0; i < this->degree; ++i)
1178 {
1179 const double L_i =
1180 weight * legendre_polynomials[i].value(
1181 quadrature_points[q_point][0]);
1182
1183 for (unsigned int j = 0; j < this->degree - 1; ++j)
1184 assembling_matrix(i * (this->degree - 1) + j,
1185 q_point) =
1186 L_i * lobatto_polynomials[j + 2].value(
1187 quadrature_points[q_point][1]);
1188 }
1189 }
1190
1191 FullMatrix<double> system_matrix(assembling_matrix.m(),
1192 assembling_matrix.m());
1193
1194 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1195 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
1196 system_matrix_inv.invert(system_matrix);
1197 }
1198
1199 solution.reinit(system_matrix_inv.m(), 8);
1200 system_rhs.reinit(system_matrix_inv.m(), 8);
1201 tmp.reinit(8);
1202
1203 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1204 {
1205 system_rhs = 0.0;
1206
1207 for (unsigned int q_point = 0; q_point < n_quadrature_points;
1208 ++q_point)
1209 {
1210 tmp = 0.0;
1211
1212 if (quadrature_points[q_point][0] < 0.5)
1213 {
1214 if (quadrature_points[q_point][1] < 0.5)
1215 {
1216 const Point<dim> quadrature_point(
1217 2.0 * quadrature_points[q_point][0],
1218 2.0 * quadrature_points[q_point][1]);
1219
1220 tmp(0) += 2.0 * this->shape_value_component(
1221 dof, quadrature_point, 0);
1222 tmp(1) += 2.0 * this->shape_value_component(
1223 dof, quadrature_point, 1);
1224 }
1225
1226 else
1227 {
1228 const Point<dim> quadrature_point(
1229 2.0 * quadrature_points[q_point][0],
1230 2.0 * quadrature_points[q_point][1] - 1.0);
1231
1232 tmp(4) += 2.0 * this->shape_value_component(
1233 dof, quadrature_point, 0);
1234 tmp(5) += 2.0 * this->shape_value_component(
1235 dof, quadrature_point, 1);
1236 }
1237 }
1238
1239 else if (quadrature_points[q_point][1] < 0.5)
1240 {
1241 const Point<dim> quadrature_point(
1242 2.0 * quadrature_points[q_point][0] - 1.0,
1243 2.0 * quadrature_points[q_point][1]);
1244
1245 tmp(2) +=
1246 2.0 * this->shape_value_component(dof,
1247 quadrature_point,
1248 0);
1249 tmp(3) +=
1250 2.0 * this->shape_value_component(dof,
1251 quadrature_point,
1252 1);
1253 }
1254
1255 else
1256 {
1257 const Point<dim> quadrature_point(
1258 2.0 * quadrature_points[q_point][0] - 1.0,
1259 2.0 * quadrature_points[q_point][1] - 1.0);
1260
1261 tmp(6) +=
1262 2.0 * this->shape_value_component(dof,
1263 quadrature_point,
1264 0);
1265 tmp(7) +=
1266 2.0 * this->shape_value_component(dof,
1267 quadrature_point,
1268 1);
1269 }
1270
1271 for (unsigned int i = 0; i < 2; ++i)
1272 for (unsigned int j = 0; j < this->degree; ++j)
1273 {
1274 tmp(2 * i) -=
1275 this->restriction[index][i](j + 2 * this->degree,
1276 dof) *
1277 this->shape_value_component(
1278 j + 2 * this->degree,
1279 quadrature_points[q_point],
1280 0);
1281 tmp(2 * i + 1) -=
1282 this->restriction[index][i](i * this->degree + j,
1283 dof) *
1284 this->shape_value_component(
1285 i * this->degree + j,
1286 quadrature_points[q_point],
1287 1);
1288 tmp(2 * (i + 2)) -= this->restriction[index][i + 2](
1289 j + 3 * this->degree, dof) *
1290 this->shape_value_component(
1291 j + 3 * this->degree,
1292 quadrature_points[q_point],
1293 0);
1294 tmp(2 * i + 5) -= this->restriction[index][i + 2](
1295 i * this->degree + j, dof) *
1296 this->shape_value_component(
1297 i * this->degree + j,
1298 quadrature_points[q_point],
1299 1);
1300 }
1301
1302 tmp *= quadrature.weight(q_point);
1303
1304 for (unsigned int i = 0; i < this->degree; ++i)
1305 {
1306 const double L_i_0 = legendre_polynomials[i].value(
1307 quadrature_points[q_point][0]);
1308 const double L_i_1 = legendre_polynomials[i].value(
1309 quadrature_points[q_point][1]);
1310
1311 for (unsigned int j = 0; j < this->degree - 1; ++j)
1312 {
1313 const double l_j_0 =
1314 L_i_0 * lobatto_polynomials[j + 2].value(
1315 quadrature_points[q_point][1]);
1316 const double l_j_1 =
1317 L_i_1 * lobatto_polynomials[j + 2].value(
1318 quadrature_points[q_point][0]);
1319
1320 for (unsigned int k = 0; k < 4; ++k)
1321 {
1322 system_rhs(i * (this->degree - 1) + j,
1323 2 * k) += tmp(2 * k) * l_j_0;
1324 system_rhs(i * (this->degree - 1) + j,
1325 2 * k + 1) +=
1326 tmp(2 * k + 1) * l_j_1;
1327 }
1328 }
1329 }
1330 }
1331
1332 system_matrix_inv.mmult(solution, system_rhs);
1333
1334 for (unsigned int i = 0; i < this->degree; ++i)
1335 for (unsigned int j = 0; j < this->degree - 1; ++j)
1336 for (unsigned int k = 0; k < 4; ++k)
1337 {
1338 if (std::abs(solution(i * (this->degree - 1) + j,
1339 2 * k)) > 1e-14)
1340 this->restriction[index][k](i * (this->degree - 1) +
1341 j + n_boundary_dofs,
1342 dof) =
1343 solution(i * (this->degree - 1) + j, 2 * k);
1344
1345 if (std::abs(solution(i * (this->degree - 1) + j,
1346 2 * k + 1)) > 1e-14)
1347 this->restriction[index][k](
1348 i + (this->degree - 1 + j) * this->degree +
1349 n_boundary_dofs,
1350 dof) =
1351 solution(i * (this->degree - 1) + j, 2 * k + 1);
1352 }
1353 }
1354 }
1355
1356 break;
1357 }
1358
1359 case 3:
1360 {
1361 // First interpolate the shape
1362 // functions of the child cells
1363 // to the lowest order shape
1364 // functions of the parent cell.
1365 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1366 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
1367 ++q_point)
1368 {
1369 const double weight = 2.0 * edge_quadrature.weight(q_point);
1370
1371 if (edge_quadrature_points[q_point][0] < 0.5)
1372 for (unsigned int i = 0; i < 2; ++i)
1373 for (unsigned int j = 0; j < 2; ++j)
1374 {
1375 Point<dim> quadrature_point(
1376 i, 2.0 * edge_quadrature_points[q_point][0], j);
1377
1378 this->restriction[index][i + 4 * j]((i + 4 * j) *
1379 this->degree,
1380 dof) +=
1381 weight *
1382 this->shape_value_component(dof, quadrature_point, 1);
1383 quadrature_point =
1384 Point<dim>(2.0 * edge_quadrature_points[q_point][0],
1385 i,
1386 j);
1387 this->restriction[index][2 * (i + 2 * j)](
1388 (i + 4 * j + 2) * this->degree, dof) +=
1389 weight *
1390 this->shape_value_component(dof, quadrature_point, 0);
1391 quadrature_point =
1392 Point<dim>(i,
1393 j,
1394 2.0 * edge_quadrature_points[q_point][0]);
1395 this->restriction[index][i + 2 * j]((i + 2 * (j + 4)) *
1396 this->degree,
1397 dof) +=
1398 weight *
1399 this->shape_value_component(dof, quadrature_point, 2);
1400 }
1401
1402 else
1403 for (unsigned int i = 0; i < 2; ++i)
1404 for (unsigned int j = 0; j < 2; ++j)
1405 {
1406 Point<dim> quadrature_point(
1407 i, 2.0 * edge_quadrature_points[q_point][0] - 1.0, j);
1408
1409 this->restriction[index][i + 4 * j + 2]((i + 4 * j) *
1410 this->degree,
1411 dof) +=
1412 weight *
1413 this->shape_value_component(dof, quadrature_point, 1);
1414 quadrature_point = Point<dim>(
1415 2.0 * edge_quadrature_points[q_point][0] - 1.0, i, j);
1416 this->restriction[index][2 * (i + 2 * j) + 1](
1417 (i + 4 * j + 2) * this->degree, dof) +=
1418 weight *
1419 this->shape_value_component(dof, quadrature_point, 0);
1420 quadrature_point = Point<dim>(
1421 i, j, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1422 this->restriction[index][i + 2 * (j + 2)](
1423 (i + 2 * (j + 4)) * this->degree, dof) +=
1424 weight *
1425 this->shape_value_component(dof, quadrature_point, 2);
1426 }
1427 }
1428
1429 // Then project the shape functions
1430 // of the child cells to the higher
1431 // order shape functions of the
1432 // parent cell.
1433 if (this->degree > 1)
1434 {
1435 const unsigned int deg = this->degree - 1;
1436 const std::vector<Polynomials::Polynomial<double>>
1437 &legendre_polynomials =
1439 FullMatrix<double> system_matrix_inv(deg, deg);
1440
1441 {
1442 FullMatrix<double> assembling_matrix(deg,
1443 n_edge_quadrature_points);
1444
1445 for (unsigned int q_point = 0;
1446 q_point < n_edge_quadrature_points;
1447 ++q_point)
1448 {
1449 const double weight =
1450 std::sqrt(edge_quadrature.weight(q_point));
1451
1452 for (unsigned int i = 0; i < deg; ++i)
1453 assembling_matrix(i, q_point) =
1454 weight * legendre_polynomials[i + 1].value(
1455 edge_quadrature_points[q_point][0]);
1456 }
1457
1458 FullMatrix<double> system_matrix(deg, deg);
1459
1460 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1461 system_matrix_inv.invert(system_matrix);
1462 }
1463
1464 FullMatrix<double> solution(deg, 6);
1465 FullMatrix<double> system_rhs(deg, 6);
1466 Vector<double> tmp(6);
1467
1468 for (unsigned int i = 0; i < 2; ++i)
1469 for (unsigned int j = 0; j < 2; ++j)
1470 for (unsigned int dof = 0; dof < this->n_dofs_per_cell();
1471 ++dof)
1472 {
1473 system_rhs = 0.0;
1474
1475 for (unsigned int q_point = 0;
1476 q_point < n_edge_quadrature_points;
1477 ++q_point)
1478 {
1479 const double weight = edge_quadrature.weight(q_point);
1480 const Point<dim> quadrature_point_0(
1481 i, edge_quadrature_points[q_point][0], j);
1482 const Point<dim> quadrature_point_1(
1483 edge_quadrature_points[q_point][0], i, j);
1484 const Point<dim> quadrature_point_2(
1485 i, j, edge_quadrature_points[q_point][0]);
1486
1487 if (edge_quadrature_points[q_point][0] < 0.5)
1488 {
1489 Point<dim> quadrature_point_3(
1490 i, 2.0 * edge_quadrature_points[q_point][0], j);
1491
1492 tmp(0) =
1493 weight * (2.0 * this->shape_value_component(
1494 dof, quadrature_point_3, 1) -
1495 this->restriction[index][i + 4 * j](
1496 (i + 4 * j) * this->degree, dof) *
1497 this->shape_value_component(
1498 (i + 4 * j) * this->degree,
1499 quadrature_point_0,
1500 1));
1501 tmp(1) =
1502 -1.0 * weight *
1503 this->restriction[index][i + 4 * j + 2](
1504 (i + 4 * j) * this->degree, dof) *
1505 this->shape_value_component((i + 4 * j) *
1506 this->degree,
1507 quadrature_point_0,
1508 1);
1509 quadrature_point_3 = Point<dim>(
1510 2.0 * edge_quadrature_points[q_point][0], i, j);
1511 tmp(2) =
1512 weight *
1513 (2.0 * this->shape_value_component(
1514 dof, quadrature_point_3, 0) -
1515 this->restriction[index][2 * (i + 2 * j)](
1516 (i + 4 * j + 2) * this->degree, dof) *
1517 this->shape_value_component(
1518 (i + 4 * j + 2) * this->degree,
1519 quadrature_point_1,
1520 0));
1521 tmp(3) =
1522 -1.0 * weight *
1523 this->restriction[index][2 * (i + 2 * j) + 1](
1524 (i + 4 * j + 2) * this->degree, dof) *
1525 this->shape_value_component((i + 4 * j + 2) *
1526 this->degree,
1527 quadrature_point_1,
1528 0);
1529 quadrature_point_3 = Point<dim>(
1530 i, j, 2.0 * edge_quadrature_points[q_point][0]);
1531 tmp(4) =
1532 weight *
1533 (2.0 * this->shape_value_component(
1534 dof, quadrature_point_3, 2) -
1535 this->restriction[index][i + 2 * j](
1536 (i + 2 * (j + 4)) * this->degree, dof) *
1537 this->shape_value_component(
1538 (i + 2 * (j + 4)) * this->degree,
1539 quadrature_point_2,
1540 2));
1541 tmp(5) =
1542 -1.0 * weight *
1543 this->restriction[index][i + 2 * (j + 2)](
1544 (i + 2 * (j + 4)) * this->degree, dof) *
1545 this->shape_value_component((i + 2 * (j + 4)) *
1546 this->degree,
1547 quadrature_point_2,
1548 2);
1549 }
1550
1551 else
1552 {
1553 tmp(0) =
1554 -1.0 * weight *
1555 this->restriction[index][i + 4 * j](
1556 (i + 4 * j) * this->degree, dof) *
1557 this->shape_value_component((i + 4 * j) *
1558 this->degree,
1559 quadrature_point_0,
1560 1);
1561
1562 Point<dim> quadrature_point_3(
1563 i,
1564 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1565 j);
1566
1567 tmp(1) = weight *
1568 (2.0 * this->shape_value_component(
1569 dof, quadrature_point_3, 1) -
1570 this->restriction[index][i + 4 * j + 2](
1571 (i + 4 * j) * this->degree, dof) *
1572 this->shape_value_component(
1573 (i + 4 * j) * this->degree,
1574 quadrature_point_0,
1575 1));
1576 tmp(2) =
1577 -1.0 * weight *
1578 this->restriction[index][2 * (i + 2 * j)](
1579 (i + 4 * j + 2) * this->degree, dof) *
1580 this->shape_value_component((i + 4 * j + 2) *
1581 this->degree,
1582 quadrature_point_1,
1583 0);
1584 quadrature_point_3 = Point<dim>(
1585 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1586 i,
1587 j);
1588 tmp(3) =
1589 weight *
1590 (2.0 * this->shape_value_component(
1591 dof, quadrature_point_3, 0) -
1592 this->restriction[index][2 * (i + 2 * j) + 1](
1593 (i + 4 * j + 2) * this->degree, dof) *
1594 this->shape_value_component(
1595 (i + 4 * j + 2) * this->degree,
1596 quadrature_point_1,
1597 0));
1598 tmp(4) =
1599 -1.0 * weight *
1600 this->restriction[index][i + 2 * j](
1601 (i + 2 * (j + 4)) * this->degree, dof) *
1602 this->shape_value_component((i + 2 * (j + 4)) *
1603 this->degree,
1604 quadrature_point_2,
1605 2);
1606 quadrature_point_3 = Point<dim>(
1607 i,
1608 j,
1609 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1610 tmp(5) =
1611 weight *
1612 (2.0 * this->shape_value_component(
1613 dof, quadrature_point_3, 2) -
1614 this->restriction[index][i + 2 * (j + 2)](
1615 (i + 2 * (j + 4)) * this->degree, dof) *
1616 this->shape_value_component(
1617 (i + 2 * (j + 4)) * this->degree,
1618 quadrature_point_2,
1619 2));
1620 }
1621
1622 for (unsigned int k = 0; k < deg; ++k)
1623 {
1624 const double L_k =
1625 legendre_polynomials[k + 1].value(
1626 edge_quadrature_points[q_point][0]);
1627
1628 for (unsigned int l = 0; l < tmp.size(); ++l)
1629 system_rhs(k, l) += tmp(l) * L_k;
1630 }
1631 }
1632
1633 system_matrix_inv.mmult(solution, system_rhs);
1634
1635 for (unsigned int k = 0; k < 2; ++k)
1636 for (unsigned int l = 0; l < deg; ++l)
1637 {
1638 if (std::abs(solution(l, k)) > 1e-14)
1639 this->restriction[index][i + 2 * (2 * j + k)](
1640 (i + 4 * j) * this->degree + l + 1, dof) =
1641 solution(l, k);
1642
1643 if (std::abs(solution(l, k + 2)) > 1e-14)
1644 this->restriction[index][2 * (i + 2 * j) + k](
1645 (i + 4 * j + 2) * this->degree + l + 1, dof) =
1646 solution(l, k + 2);
1647
1648 if (std::abs(solution(l, k + 4)) > 1e-14)
1649 this->restriction[index][i + 2 * (j + 2 * k)](
1650 (i + 2 * (j + 4)) * this->degree + l + 1, dof) =
1651 solution(l, k + 4);
1652 }
1653 }
1654
1655 const QGauss<2> face_quadrature(2 * this->degree);
1656 const std::vector<Point<2>> &face_quadrature_points =
1657 face_quadrature.get_points();
1658 const std::vector<Polynomials::Polynomial<double>>
1659 &lobatto_polynomials =
1661 const unsigned int n_edge_dofs =
1662 GeometryInfo<dim>::lines_per_cell * this->degree;
1663 const unsigned int n_face_quadrature_points =
1664 face_quadrature.size();
1665
1666 {
1667 FullMatrix<double> assembling_matrix(deg * this->degree,
1668 n_face_quadrature_points);
1669
1670 for (unsigned int q_point = 0;
1671 q_point < n_face_quadrature_points;
1672 ++q_point)
1673 {
1674 const double weight =
1675 std::sqrt(face_quadrature.weight(q_point));
1676
1677 for (unsigned int i = 0; i <= deg; ++i)
1678 {
1679 const double L_i =
1680 weight * legendre_polynomials[i].value(
1681 face_quadrature_points[q_point][0]);
1682
1683 for (unsigned int j = 0; j < deg; ++j)
1684 assembling_matrix(i * deg + j, q_point) =
1685 L_i * lobatto_polynomials[j + 2].value(
1686 face_quadrature_points[q_point][1]);
1687 }
1688 }
1689
1690 FullMatrix<double> system_matrix(assembling_matrix.m(),
1691 assembling_matrix.m());
1692
1693 assembling_matrix.mTmult(system_matrix, assembling_matrix);
1694 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
1695 system_matrix_inv.invert(system_matrix);
1696 }
1697
1698 solution.reinit(system_matrix_inv.m(), 24);
1699 system_rhs.reinit(system_matrix_inv.m(), 24);
1700 tmp.reinit(24);
1701
1702 for (unsigned int i = 0; i < 2; ++i)
1703 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1704 {
1705 system_rhs = 0.0;
1706
1707 for (unsigned int q_point = 0;
1708 q_point < n_face_quadrature_points;
1709 ++q_point)
1710 {
1711 tmp = 0.0;
1712
1713 if (face_quadrature_points[q_point][0] < 0.5)
1714 {
1715 if (face_quadrature_points[q_point][1] < 0.5)
1716 {
1717 Point<dim> quadrature_point_0(
1718 i,
1719 2.0 * face_quadrature_points[q_point][0],
1720 2.0 * face_quadrature_points[q_point][1]);
1721
1722 tmp(0) += 2.0 * this->shape_value_component(
1723 dof, quadrature_point_0, 1);
1724 tmp(1) += 2.0 * this->shape_value_component(
1725 dof, quadrature_point_0, 2);
1726 quadrature_point_0 = Point<dim>(
1727 2.0 * face_quadrature_points[q_point][0],
1728 i,
1729 2.0 * face_quadrature_points[q_point][1]);
1730 tmp(8) += 2.0 * this->shape_value_component(
1731 dof, quadrature_point_0, 2);
1732 tmp(9) += 2.0 * this->shape_value_component(
1733 dof, quadrature_point_0, 0);
1734 quadrature_point_0 = Point<dim>(
1735 2.0 * face_quadrature_points[q_point][0],
1736 2.0 * face_quadrature_points[q_point][1],
1737 i);
1738 tmp(16) += 2.0 * this->shape_value_component(
1739 dof, quadrature_point_0, 0);
1740 tmp(17) += 2.0 * this->shape_value_component(
1741 dof, quadrature_point_0, 1);
1742 }
1743
1744 else
1745 {
1746 Point<dim> quadrature_point_0(
1747 i,
1748 2.0 * face_quadrature_points[q_point][0],
1749 2.0 * face_quadrature_points[q_point][1] -
1750 1.0);
1751
1752 tmp(2) += 2.0 * this->shape_value_component(
1753 dof, quadrature_point_0, 1);
1754 tmp(3) += 2.0 * this->shape_value_component(
1755 dof, quadrature_point_0, 2);
1756 quadrature_point_0 = Point<dim>(
1757 2.0 * face_quadrature_points[q_point][0],
1758 i,
1759 2.0 * face_quadrature_points[q_point][1] -
1760 1.0);
1761 tmp(10) += 2.0 * this->shape_value_component(
1762 dof, quadrature_point_0, 2);
1763 tmp(11) += 2.0 * this->shape_value_component(
1764 dof, quadrature_point_0, 0);
1765 quadrature_point_0 = Point<dim>(
1766 2.0 * face_quadrature_points[q_point][0],
1767 2.0 * face_quadrature_points[q_point][1] -
1768 1.0,
1769 i);
1770 tmp(18) += 2.0 * this->shape_value_component(
1771 dof, quadrature_point_0, 0);
1772 tmp(19) += 2.0 * this->shape_value_component(
1773 dof, quadrature_point_0, 1);
1774 }
1775 }
1776
1777 else if (face_quadrature_points[q_point][1] < 0.5)
1778 {
1779 Point<dim> quadrature_point_0(
1780 i,
1781 2.0 * face_quadrature_points[q_point][0] - 1.0,
1782 2.0 * face_quadrature_points[q_point][1]);
1783
1784 tmp(4) += 2.0 * this->shape_value_component(
1785 dof, quadrature_point_0, 1);
1786 tmp(5) += 2.0 * this->shape_value_component(
1787 dof, quadrature_point_0, 2);
1788 quadrature_point_0 = Point<dim>(
1789 2.0 * face_quadrature_points[q_point][0] - 1.0,
1790 i,
1791 2.0 * face_quadrature_points[q_point][1]);
1792 tmp(12) += 2.0 * this->shape_value_component(
1793 dof, quadrature_point_0, 2);
1794 tmp(13) += 2.0 * this->shape_value_component(
1795 dof, quadrature_point_0, 0);
1796 quadrature_point_0 = Point<dim>(
1797 2.0 * face_quadrature_points[q_point][0] - 1.0,
1798 2.0 * face_quadrature_points[q_point][1],
1799 i);
1800 tmp(20) += 2.0 * this->shape_value_component(
1801 dof, quadrature_point_0, 0);
1802 tmp(21) += 2.0 * this->shape_value_component(
1803 dof, quadrature_point_0, 1);
1804 }
1805
1806 else
1807 {
1808 Point<dim> quadrature_point_0(
1809 i,
1810 2.0 * face_quadrature_points[q_point][0] - 1.0,
1811 2.0 * face_quadrature_points[q_point][1] - 1.0);
1812
1813 tmp(6) += 2.0 * this->shape_value_component(
1814 dof, quadrature_point_0, 1);
1815 tmp(7) += 2.0 * this->shape_value_component(
1816 dof, quadrature_point_0, 2);
1817 quadrature_point_0 = Point<dim>(
1818 2.0 * face_quadrature_points[q_point][0] - 1.0,
1819 i,
1820 2.0 * face_quadrature_points[q_point][1] - 1.0);
1821 tmp(14) += 2.0 * this->shape_value_component(
1822 dof, quadrature_point_0, 2);
1823 tmp(15) += 2.0 * this->shape_value_component(
1824 dof, quadrature_point_0, 0);
1825 quadrature_point_0 = Point<dim>(
1826 2.0 * face_quadrature_points[q_point][0] - 1.0,
1827 2.0 * face_quadrature_points[q_point][1] - 1.0,
1828 i);
1829 tmp(22) += 2.0 * this->shape_value_component(
1830 dof, quadrature_point_0, 0);
1831 tmp(23) += 2.0 * this->shape_value_component(
1832 dof, quadrature_point_0, 1);
1833 }
1834
1835 const Point<dim> quadrature_point_0(
1836 i,
1837 face_quadrature_points[q_point][0],
1838 face_quadrature_points[q_point][1]);
1839 const Point<dim> quadrature_point_1(
1840 face_quadrature_points[q_point][0],
1841 i,
1842 face_quadrature_points[q_point][1]);
1843 const Point<dim> quadrature_point_2(
1844 face_quadrature_points[q_point][0],
1845 face_quadrature_points[q_point][1],
1846 i);
1847
1848 for (unsigned int j = 0; j < 2; ++j)
1849 for (unsigned int k = 0; k < 2; ++k)
1850 for (unsigned int l = 0; l <= deg; ++l)
1851 {
1852 tmp(2 * (j + 2 * k)) -=
1853 this->restriction[index][i + 2 * (2 * j + k)](
1854 (i + 4 * j) * this->degree + l, dof) *
1855 this->shape_value_component(
1856 (i + 4 * j) * this->degree + l,
1857 quadrature_point_0,
1858 1);
1859 tmp(2 * (j + 2 * k) + 1) -=
1860 this->restriction[index][i + 2 * (2 * j + k)](
1861 (i + 2 * (k + 4)) * this->degree + l, dof) *
1862 this->shape_value_component(
1863 (i + 2 * (k + 4)) * this->degree + l,
1864 quadrature_point_0,
1865 2);
1866 tmp(2 * (j + 2 * (k + 2))) -=
1867 this->restriction[index][2 * (i + 2 * j) + k](
1868 (2 * (i + 4) + k) * this->degree + l, dof) *
1869 this->shape_value_component(
1870 (2 * (i + 4) + k) * this->degree + l,
1871 quadrature_point_1,
1872 2);
1873 tmp(2 * (j + 2 * k) + 9) -=
1874 this->restriction[index][2 * (i + 2 * j) + k](
1875 (i + 4 * j + 2) * this->degree + l, dof) *
1876 this->shape_value_component(
1877 (i + 4 * j + 2) * this->degree + l,
1878 quadrature_point_1,
1879 0);
1880 tmp(2 * (j + 2 * (k + 4))) -=
1881 this->restriction[index][2 * (2 * i + j) + k](
1882 (4 * i + j + 2) * this->degree + l, dof) *
1883 this->shape_value_component(
1884 (4 * i + j + 2) * this->degree + l,
1885 quadrature_point_2,
1886 0);
1887 tmp(2 * (j + 2 * k) + 17) -=
1888 this->restriction[index][2 * (2 * i + j) + k](
1889 (4 * i + k) * this->degree + l, dof) *
1890 this->shape_value_component(
1891 (4 * i + k) * this->degree + l,
1892 quadrature_point_2,
1893 1);
1894 }
1895
1896 tmp *= face_quadrature.weight(q_point);
1897
1898 for (unsigned int j = 0; j <= deg; ++j)
1899 {
1900 const double L_j_0 = legendre_polynomials[j].value(
1901 face_quadrature_points[q_point][0]);
1902 const double L_j_1 = legendre_polynomials[j].value(
1903 face_quadrature_points[q_point][1]);
1904
1905 for (unsigned int k = 0; k < deg; ++k)
1906 {
1907 const double l_k_0 =
1908 L_j_0 * lobatto_polynomials[k + 2].value(
1909 face_quadrature_points[q_point][1]);
1910 const double l_k_1 =
1911 L_j_1 * lobatto_polynomials[k + 2].value(
1912 face_quadrature_points[q_point][0]);
1913
1914 for (unsigned int l = 0; l < 4; ++l)
1915 {
1916 system_rhs(j * deg + k, 2 * l) +=
1917 tmp(2 * l) * l_k_0;
1918 system_rhs(j * deg + k, 2 * l + 1) +=
1919 tmp(2 * l + 1) * l_k_1;
1920 system_rhs(j * deg + k, 2 * (l + 4)) +=
1921 tmp(2 * (l + 4)) * l_k_1;
1922 system_rhs(j * deg + k, 2 * l + 9) +=
1923 tmp(2 * l + 9) * l_k_0;
1924 system_rhs(j * deg + k, 2 * (l + 8)) +=
1925 tmp(2 * (l + 8)) * l_k_0;
1926 system_rhs(j * deg + k, 2 * l + 17) +=
1927 tmp(2 * l + 17) * l_k_1;
1928 }
1929 }
1930 }
1931 }
1932
1933 system_matrix_inv.mmult(solution, system_rhs);
1934
1935 for (unsigned int j = 0; j < 2; ++j)
1936 for (unsigned int k = 0; k < 2; ++k)
1937 for (unsigned int l = 0; l <= deg; ++l)
1938 for (unsigned int m = 0; m < deg; ++m)
1939 {
1940 if (std::abs(solution(l * deg + m,
1941 2 * (j + 2 * k))) > 1e-14)
1942 this->restriction[index][i + 2 * (2 * j + k)](
1943 (2 * i * this->degree + l) * deg + m +
1944 n_edge_dofs,
1945 dof) = solution(l * deg + m, 2 * (j + 2 * k));
1946
1947 if (std::abs(solution(l * deg + m,
1948 2 * (j + 2 * k) + 1)) >
1949 1e-14)
1950 this->restriction[index][i + 2 * (2 * j + k)](
1951 ((2 * i + 1) * deg + m) * this->degree + l +
1952 n_edge_dofs,
1953 dof) =
1954 solution(l * deg + m, 2 * (j + 2 * k) + 1);
1955
1956 if (std::abs(solution(l * deg + m,
1957 2 * (j + 2 * (k + 2)))) >
1958 1e-14)
1959 this->restriction[index][2 * (i + 2 * j) + k](
1960 (2 * (i + 2) * this->degree + l) * deg + m +
1961 n_edge_dofs,
1962 dof) =
1963 solution(l * deg + m, 2 * (j + 2 * (k + 2)));
1964
1965 if (std::abs(solution(l * deg + m,
1966 2 * (j + 2 * k) + 9)) >
1967 1e-14)
1968 this->restriction[index][2 * (i + 2 * j) + k](
1969 ((2 * i + 5) * deg + m) * this->degree + l +
1970 n_edge_dofs,
1971 dof) =
1972 solution(l * deg + m, 2 * (j + 2 * k) + 9);
1973
1974 if (std::abs(solution(l * deg + m,
1975 2 * (j + 2 * (k + 4)))) >
1976 1e-14)
1977 this->restriction[index][2 * (2 * i + j) + k](
1978 (2 * (i + 4) * this->degree + l) * deg + m +
1979 n_edge_dofs,
1980 dof) =
1981 solution(l * deg + m, 2 * (j + 2 * (k + 4)));
1982
1983 if (std::abs(solution(l * deg + m,
1984 2 * (j + 2 * k) + 17)) >
1985 1e-14)
1986 this->restriction[index][2 * (2 * i + j) + k](
1987 ((2 * i + 9) * deg + m) * this->degree + l +
1988 n_edge_dofs,
1989 dof) =
1990 solution(l * deg + m, 2 * (j + 2 * k) + 17);
1991 }
1992 }
1993
1994 const QGauss<dim> quadrature(2 * this->degree);
1995 const std::vector<Point<dim>> &quadrature_points =
1996 quadrature.get_points();
1997 const unsigned int n_boundary_dofs =
1998 2 * GeometryInfo<dim>::faces_per_cell * deg * this->degree +
1999 n_edge_dofs;
2000 const unsigned int n_quadrature_points = quadrature.size();
2001
2002 {
2003 FullMatrix<double> assembling_matrix(deg * deg * this->degree,
2004 n_quadrature_points);
2005
2006 for (unsigned int q_point = 0; q_point < n_quadrature_points;
2007 ++q_point)
2008 {
2009 const double weight = std::sqrt(quadrature.weight(q_point));
2010
2011 for (unsigned int i = 0; i <= deg; ++i)
2012 {
2013 const double L_i =
2014 weight * legendre_polynomials[i].value(
2015 quadrature_points[q_point][0]);
2016
2017 for (unsigned int j = 0; j < deg; ++j)
2018 {
2019 const double l_j =
2020 L_i * lobatto_polynomials[j + 2].value(
2021 quadrature_points[q_point][1]);
2022
2023 for (unsigned int k = 0; k < deg; ++k)
2024 assembling_matrix((i * deg + j) * deg + k,
2025 q_point) =
2026 l_j * lobatto_polynomials[k + 2].value(
2027 quadrature_points[q_point][2]);
2028 }
2029 }
2030 }
2031
2032 FullMatrix<double> system_matrix(assembling_matrix.m(),
2033 assembling_matrix.m());
2034
2035 assembling_matrix.mTmult(system_matrix, assembling_matrix);
2036 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
2037 system_matrix_inv.invert(system_matrix);
2038 }
2039
2040 solution.reinit(system_matrix_inv.m(), 24);
2041 system_rhs.reinit(system_matrix_inv.m(), 24);
2042 tmp.reinit(24);
2043
2044 for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2045 {
2046 system_rhs = 0.0;
2047
2048 for (unsigned int q_point = 0; q_point < n_quadrature_points;
2049 ++q_point)
2050 {
2051 tmp = 0.0;
2052
2053 if (quadrature_points[q_point][0] < 0.5)
2054 {
2055 if (quadrature_points[q_point][1] < 0.5)
2056 {
2057 if (quadrature_points[q_point][2] < 0.5)
2058 {
2059 const Point<dim> quadrature_point(
2060 2.0 * quadrature_points[q_point][0],
2061 2.0 * quadrature_points[q_point][1],
2062 2.0 * quadrature_points[q_point][2]);
2063
2064 tmp(0) += 2.0 * this->shape_value_component(
2065 dof, quadrature_point, 0);
2066 tmp(1) += 2.0 * this->shape_value_component(
2067 dof, quadrature_point, 1);
2068 tmp(2) += 2.0 * this->shape_value_component(
2069 dof, quadrature_point, 2);
2070 }
2071
2072 else
2073 {
2074 const Point<dim> quadrature_point(
2075 2.0 * quadrature_points[q_point][0],
2076 2.0 * quadrature_points[q_point][1],
2077 2.0 * quadrature_points[q_point][2] - 1.0);
2078
2079 tmp(3) += 2.0 * this->shape_value_component(
2080 dof, quadrature_point, 0);
2081 tmp(4) += 2.0 * this->shape_value_component(
2082 dof, quadrature_point, 1);
2083 tmp(5) += 2.0 * this->shape_value_component(
2084 dof, quadrature_point, 2);
2085 }
2086 }
2087
2088 else if (quadrature_points[q_point][2] < 0.5)
2089 {
2090 const Point<dim> quadrature_point(
2091 2.0 * quadrature_points[q_point][0],
2092 2.0 * quadrature_points[q_point][1] - 1.0,
2093 2.0 * quadrature_points[q_point][2]);
2094
2095 tmp(6) += 2.0 * this->shape_value_component(
2096 dof, quadrature_point, 0);
2097 tmp(7) += 2.0 * this->shape_value_component(
2098 dof, quadrature_point, 1);
2099 tmp(8) += 2.0 * this->shape_value_component(
2100 dof, quadrature_point, 2);
2101 }
2102
2103 else
2104 {
2105 const Point<dim> quadrature_point(
2106 2.0 * quadrature_points[q_point][0],
2107 2.0 * quadrature_points[q_point][1] - 1.0,
2108 2.0 * quadrature_points[q_point][2] - 1.0);
2109
2110 tmp(9) += 2.0 * this->shape_value_component(
2111 dof, quadrature_point, 0);
2112 tmp(10) += 2.0 * this->shape_value_component(
2113 dof, quadrature_point, 1);
2114 tmp(11) += 2.0 * this->shape_value_component(
2115 dof, quadrature_point, 2);
2116 }
2117 }
2118
2119 else if (quadrature_points[q_point][1] < 0.5)
2120 {
2121 if (quadrature_points[q_point][2] < 0.5)
2122 {
2123 const Point<dim> quadrature_point(
2124 2.0 * quadrature_points[q_point][0] - 1.0,
2125 2.0 * quadrature_points[q_point][1],
2126 2.0 * quadrature_points[q_point][2]);
2127
2128 tmp(12) += 2.0 * this->shape_value_component(
2129 dof, quadrature_point, 0);
2130 tmp(13) += 2.0 * this->shape_value_component(
2131 dof, quadrature_point, 1);
2132 tmp(14) += 2.0 * this->shape_value_component(
2133 dof, quadrature_point, 2);
2134 }
2135
2136 else
2137 {
2138 const Point<dim> quadrature_point(
2139 2.0 * quadrature_points[q_point][0] - 1.0,
2140 2.0 * quadrature_points[q_point][1],
2141 2.0 * quadrature_points[q_point][2] - 1.0);
2142
2143 tmp(15) += 2.0 * this->shape_value_component(
2144 dof, quadrature_point, 0);
2145 tmp(16) += 2.0 * this->shape_value_component(
2146 dof, quadrature_point, 1);
2147 tmp(17) += 2.0 * this->shape_value_component(
2148 dof, quadrature_point, 2);
2149 }
2150 }
2151
2152 else if (quadrature_points[q_point][2] < 0.5)
2153 {
2154 const Point<dim> quadrature_point(
2155 2.0 * quadrature_points[q_point][0] - 1.0,
2156 2.0 * quadrature_points[q_point][1] - 1.0,
2157 2.0 * quadrature_points[q_point][2]);
2158
2159 tmp(18) +=
2160 2.0 * this->shape_value_component(dof,
2161 quadrature_point,
2162 0);
2163 tmp(19) +=
2164 2.0 * this->shape_value_component(dof,
2165 quadrature_point,
2166 1);
2167 tmp(20) +=
2168 2.0 * this->shape_value_component(dof,
2169 quadrature_point,
2170 2);
2171 }
2172
2173 else
2174 {
2175 const Point<dim> quadrature_point(
2176 2.0 * quadrature_points[q_point][0] - 1.0,
2177 2.0 * quadrature_points[q_point][1] - 1.0,
2178 2.0 * quadrature_points[q_point][2] - 1.0);
2179
2180 tmp(21) +=
2181 2.0 * this->shape_value_component(dof,
2182 quadrature_point,
2183 0);
2184 tmp(22) +=
2185 2.0 * this->shape_value_component(dof,
2186 quadrature_point,
2187 1);
2188 tmp(23) +=
2189 2.0 * this->shape_value_component(dof,
2190 quadrature_point,
2191 2);
2192 }
2193
2194 for (unsigned int i = 0; i < 2; ++i)
2195 for (unsigned int j = 0; j < 2; ++j)
2196 for (unsigned int k = 0; k < 2; ++k)
2197 for (unsigned int l = 0; l <= deg; ++l)
2198 {
2199 tmp(3 * (i + 2 * (j + 2 * k))) -=
2200 this->restriction[index][2 * (2 * i + j) + k](
2201 (4 * i + j + 2) * this->degree + l, dof) *
2202 this->shape_value_component(
2203 (4 * i + j + 2) * this->degree + l,
2204 quadrature_points[q_point],
2205 0);
2206 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
2207 this->restriction[index][2 * (2 * i + j) + k](
2208 (4 * i + k) * this->degree + l, dof) *
2209 this->shape_value_component(
2210 (4 * i + k) * this->degree + l,
2211 quadrature_points[q_point],
2212 1);
2213 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
2214 this->restriction[index][2 * (2 * i + j) + k](
2215 (2 * (j + 4) + k) * this->degree + l, dof) *
2216 this->shape_value_component(
2217 (2 * (j + 4) + k) * this->degree + l,
2218 quadrature_points[q_point],
2219 2);
2220
2221 for (unsigned int m = 0; m < deg; ++m)
2222 {
2223 tmp(3 * (i + 2 * (j + 2 * k))) -=
2224 this->restriction[index][2 * (2 * i + j) +
2225 k](
2226 ((2 * j + 5) * deg + m) * this->degree +
2227 l + n_edge_dofs,
2228 dof) *
2229 this->shape_value_component(
2230 ((2 * j + 5) * deg + m) * this->degree +
2231 l + n_edge_dofs,
2232 quadrature_points[q_point],
2233 0);
2234 tmp(3 * (i + 2 * (j + 2 * k))) -=
2235 this->restriction[index][2 * (2 * i + j) +
2236 k](
2237 (2 * (i + 4) * this->degree + l) * deg +
2238 m + n_edge_dofs,
2239 dof) *
2240 this->shape_value_component(
2241 (2 * (i + 4) * this->degree + l) * deg +
2242 m + n_edge_dofs,
2243 quadrature_points[q_point],
2244 0);
2245 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
2246 this->restriction[index][2 * (2 * i + j) +
2247 k](
2248 (2 * k * this->degree + l) * deg + m +
2249 n_edge_dofs,
2250 dof) *
2251 this->shape_value_component(
2252 (2 * k * this->degree + l) * deg + m +
2253 n_edge_dofs,
2254 quadrature_points[q_point],
2255 1);
2256 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
2257 this->restriction[index][2 * (2 * i + j) +
2258 k](
2259 ((2 * i + 9) * deg + m) * this->degree +
2260 l + n_edge_dofs,
2261 dof) *
2262 this->shape_value_component(
2263 ((2 * i + 9) * deg + m) * this->degree +
2264 l + n_edge_dofs,
2265 quadrature_points[q_point],
2266 1);
2267 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
2268 this->restriction[index][2 * (2 * i + j) +
2269 k](
2270 ((2 * k + 1) * deg + m) * this->degree +
2271 l + n_edge_dofs,
2272 dof) *
2273 this->shape_value_component(
2274 ((2 * k + 1) * deg + m) * this->degree +
2275 l + n_edge_dofs,
2276 quadrature_points[q_point],
2277 2);
2278 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
2279 this->restriction[index][2 * (2 * i + j) +
2280 k](
2281 (2 * (j + 2) * this->degree + l) * deg +
2282 m + n_edge_dofs,
2283 dof) *
2284 this->shape_value_component(
2285 (2 * (j + 2) * this->degree + l) * deg +
2286 m + n_edge_dofs,
2287 quadrature_points[q_point],
2288 2);
2289 }
2290 }
2291
2292 tmp *= quadrature.weight(q_point);
2293
2294 for (unsigned int i = 0; i <= deg; ++i)
2295 {
2296 const double L_i_0 = legendre_polynomials[i].value(
2297 quadrature_points[q_point][0]);
2298 const double L_i_1 = legendre_polynomials[i].value(
2299 quadrature_points[q_point][1]);
2300 const double L_i_2 = legendre_polynomials[i].value(
2301 quadrature_points[q_point][2]);
2302
2303 for (unsigned int j = 0; j < deg; ++j)
2304 {
2305 const double l_j_0 =
2306 L_i_0 * lobatto_polynomials[j + 2].value(
2307 quadrature_points[q_point][1]);
2308 const double l_j_1 =
2309 L_i_1 * lobatto_polynomials[j + 2].value(
2310 quadrature_points[q_point][0]);
2311 const double l_j_2 =
2312 L_i_2 * lobatto_polynomials[j + 2].value(
2313 quadrature_points[q_point][0]);
2314
2315 for (unsigned int k = 0; k < deg; ++k)
2316 {
2317 const double l_k_0 =
2318 l_j_0 * lobatto_polynomials[k + 2].value(
2319 quadrature_points[q_point][2]);
2320 const double l_k_1 =
2321 l_j_1 * lobatto_polynomials[k + 2].value(
2322 quadrature_points[q_point][2]);
2323 const double l_k_2 =
2324 l_j_2 * lobatto_polynomials[k + 2].value(
2325 quadrature_points[q_point][1]);
2326
2327 for (unsigned int l = 0; l < 8; ++l)
2328 {
2329 system_rhs((i * deg + j) * deg + k,
2330 3 * l) += tmp(3 * l) * l_k_0;
2331 system_rhs((i * deg + j) * deg + k,
2332 3 * l + 1) +=
2333 tmp(3 * l + 1) * l_k_1;
2334 system_rhs((i * deg + j) * deg + k,
2335 3 * l + 2) +=
2336 tmp(3 * l + 2) * l_k_2;
2337 }
2338 }
2339 }
2340 }
2341 }
2342
2343 system_matrix_inv.mmult(solution, system_rhs);
2344
2345 for (unsigned int i = 0; i < 2; ++i)
2346 for (unsigned int j = 0; j < 2; ++j)
2347 for (unsigned int k = 0; k < 2; ++k)
2348 for (unsigned int l = 0; l <= deg; ++l)
2349 for (unsigned int m = 0; m < deg; ++m)
2350 for (unsigned int n = 0; n < deg; ++n)
2351 {
2352 if (std::abs(
2353 solution((l * deg + m) * deg + n,
2354 3 * (i + 2 * (j + 2 * k)))) >
2355 1e-14)
2356 this->restriction[index][2 * (2 * i + j) + k](
2357 (l * deg + m) * deg + n + n_boundary_dofs,
2358 dof) = solution((l * deg + m) * deg + n,
2359 3 * (i + 2 * (j + 2 * k)));
2360
2361 if (std::abs(
2362 solution((l * deg + m) * deg + n,
2363 3 * (i + 2 * (j + 2 * k)) + 1)) >
2364 1e-14)
2365 this->restriction[index][2 * (2 * i + j) + k](
2366 (l + (m + deg) * this->degree) * deg + n +
2367 n_boundary_dofs,
2368 dof) =
2369 solution((l * deg + m) * deg + n,
2370 3 * (i + 2 * (j + 2 * k)) + 1);
2371
2372 if (std::abs(
2373 solution((l * deg + m) * deg + n,
2374 3 * (i + 2 * (j + 2 * k)) + 2)) >
2375 1e-14)
2376 this->restriction[index][2 * (2 * i + j) + k](
2377 l +
2378 ((m + 2 * deg) * deg + n) * this->degree +
2379 n_boundary_dofs,
2380 dof) =
2381 solution((l * deg + m) * deg + n,
2382 3 * (i + 2 * (j + 2 * k)) + 2);
2383 }
2384 }
2385 }
2386
2387 break;
2388 }
2389
2390 default:
2392 }
2393}
2394
2395
2396
2397template <int dim>
2398std::vector<unsigned int>
2399FE_Nedelec<dim>::get_dpo_vector(const unsigned int degree, bool dg)
2400{
2401 std::vector<unsigned int> dpo;
2402
2403 if (dg)
2404 {
2405 dpo.resize(dim + 1);
2406 dpo[dim] = PolynomialsNedelec<dim>::n_polynomials(degree);
2407 }
2408 else
2409 {
2410 dpo.push_back(0);
2411 dpo.push_back(degree + 1);
2412 if (dim > 1)
2413 dpo.push_back(2 * degree * (degree + 1));
2414 if (dim > 2)
2415 dpo.push_back(3 * degree * degree * (degree + 1));
2416 }
2417
2418 return dpo;
2419}
2420
2421//---------------------------------------------------------------------------
2422// Data field initialization
2423//---------------------------------------------------------------------------
2424
2425// Check whether a given shape
2426// function has support on a
2427// given face.
2428
2429// We just switch through the
2430// faces of the cell and return
2431// true, if the shape function
2432// has support on the face
2433// and false otherwise.
2434template <int dim>
2435bool
2436FE_Nedelec<dim>::has_support_on_face(const unsigned int shape_index,
2437 const unsigned int face_index) const
2438{
2439 AssertIndexRange(shape_index, this->n_dofs_per_cell());
2441
2442 const unsigned int deg = this->degree - 1;
2443 switch (dim)
2444 {
2445 case 2:
2446 switch (face_index)
2447 {
2448 case 0:
2449 if (!((shape_index > deg) && (shape_index < 2 * this->degree)))
2450 return true;
2451
2452 else
2453 return false;
2454
2455 case 1:
2456 if ((shape_index > deg) &&
2457 (shape_index <
2458 GeometryInfo<2>::lines_per_cell * this->degree))
2459 return true;
2460
2461 else
2462 return false;
2463
2464 case 2:
2465 if (shape_index < 3 * this->degree)
2466 return true;
2467
2468 else
2469 return false;
2470
2471 case 3:
2472 if (!((shape_index >= 2 * this->degree) &&
2473 (shape_index < 3 * this->degree)))
2474 return true;
2475
2476 else
2477 return false;
2478
2479 default:
2480 {
2482 return false;
2483 }
2484 }
2485
2486 case 3:
2487 switch (face_index)
2488 {
2489 case 0:
2490 if (((shape_index > deg) && (shape_index < 2 * this->degree)) ||
2491 ((shape_index >= 5 * this->degree) &&
2492 (shape_index < 6 * this->degree)) ||
2493 ((shape_index >= 9 * this->degree) &&
2494 (shape_index < 10 * this->degree)) ||
2495 ((shape_index >= 11 * this->degree) &&
2496 (shape_index <
2497 GeometryInfo<3>::lines_per_cell * this->degree)) ||
2498 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2499 this->degree) &&
2500 (shape_index < (GeometryInfo<3>::lines_per_cell + 4 * deg) *
2501 this->degree)) ||
2502 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2503 this->degree) &&
2504 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2505 this->degree)) ||
2506 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2507 this->degree) &&
2508 (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2509 this->degree)) ||
2510 ((shape_index >=
2511 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2512 this->degree) &&
2513 (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2514 this->degree)))
2515 return false;
2516
2517 else
2518 return true;
2519
2520 case 1:
2521 if (((shape_index > deg) && (shape_index < 4 * this->degree)) ||
2522 ((shape_index >= 5 * this->degree) &&
2523 (shape_index < 8 * this->degree)) ||
2524 ((shape_index >= 9 * this->degree) &&
2525 (shape_index < 10 * this->degree)) ||
2526 ((shape_index >= 11 * this->degree) &&
2527 (shape_index <
2528 GeometryInfo<3>::lines_per_cell * this->degree)) ||
2529 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2530 this->degree) &&
2531 (shape_index < (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2532 this->degree)) ||
2533 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2534 this->degree) &&
2535 (shape_index < (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2536 this->degree)) ||
2537 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2538 this->degree) &&
2539 (shape_index < (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2540 this->degree)) ||
2541 ((shape_index >=
2542 (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2543 this->degree) &&
2544 (shape_index < (GeometryInfo<3>::lines_per_cell + 12 * deg) *
2545 this->degree)))
2546 return true;
2547
2548 else
2549 return false;
2550
2551 case 2:
2552 if ((shape_index < 3 * this->degree) ||
2553 ((shape_index >= 4 * this->degree) &&
2554 (shape_index < 7 * this->degree)) ||
2555 ((shape_index >= 8 * this->degree) &&
2556 (shape_index < 10 * this->degree)) ||
2557 ((shape_index >=
2558 (GeometryInfo<3>::lines_per_cell + deg) * this->degree) &&
2559 (shape_index < (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2560 this->degree)) ||
2561 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2562 this->degree) &&
2563 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2564 this->degree)) ||
2565 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 8 * deg) *
2566 this->degree) &&
2567 (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2568 this->degree)) ||
2569 ((shape_index >=
2570 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2571 this->degree) &&
2572 (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2573 this->degree)))
2574 return true;
2575
2576 else
2577 return false;
2578
2579 case 3:
2580 if ((shape_index < 2 * this->degree) ||
2581 ((shape_index >= 3 * this->degree) &&
2582 (shape_index < 6 * this->degree)) ||
2583 ((shape_index >= 7 * this->degree) &&
2584 (shape_index < 8 * this->degree)) ||
2585 ((shape_index >= 10 * this->degree) &&
2586 (shape_index <
2587 GeometryInfo<3>::lines_per_cell * this->degree)) ||
2588 ((shape_index >=
2589 (GeometryInfo<3>::lines_per_cell + deg) * this->degree) &&
2590 (shape_index < (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2591 this->degree)) ||
2592 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2593 this->degree) &&
2594 (shape_index < (GeometryInfo<3>::lines_per_cell + 4 * deg) *
2595 this->degree)) ||
2596 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2597 this->degree) &&
2598 (shape_index < (GeometryInfo<3>::lines_per_cell + 9 * deg) *
2599 this->degree)) ||
2600 ((shape_index >=
2601 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2602 this->degree) &&
2603 (shape_index < (GeometryInfo<3>::lines_per_cell + 11 * deg) *
2604 this->degree)))
2605 return true;
2606
2607 else
2608 return false;
2609
2610 case 4:
2611 if ((shape_index < 4 * this->degree) ||
2612 ((shape_index >= 8 * this->degree) &&
2613 (shape_index <
2614 (GeometryInfo<3>::lines_per_cell + deg) * this->degree)) ||
2615 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2616 this->degree) &&
2617 (shape_index < (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2618 this->degree)) ||
2619 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2620 this->degree) &&
2621 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2622 this->degree)) ||
2623 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2624 this->degree) &&
2625 (shape_index < (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2626 this->degree)))
2627 return true;
2628
2629 else
2630 return false;
2631
2632 case 5:
2633 if (((shape_index >= 4 * this->degree) &&
2634 (shape_index <
2635 (GeometryInfo<3>::lines_per_cell + deg) * this->degree)) ||
2636 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 2 * deg) *
2637 this->degree) &&
2638 (shape_index < (GeometryInfo<3>::lines_per_cell + 3 * deg) *
2639 this->degree)) ||
2640 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 5 * deg) *
2641 this->degree) &&
2642 (shape_index < (GeometryInfo<3>::lines_per_cell + 6 * deg) *
2643 this->degree)) ||
2644 ((shape_index >= (GeometryInfo<3>::lines_per_cell + 7 * deg) *
2645 this->degree) &&
2646 (shape_index < (GeometryInfo<3>::lines_per_cell + 8 * deg) *
2647 this->degree)) ||
2648 ((shape_index >=
2649 (GeometryInfo<3>::lines_per_cell + 10 * deg) *
2650 this->degree) &&
2651 (shape_index < (GeometryInfo<3>::lines_per_cell + 12 * deg) *
2652 this->degree)))
2653 return true;
2654
2655 else
2656 return false;
2657
2658 default:
2659 {
2661 return false;
2662 }
2663 }
2664
2665 default:
2666 {
2668 return false;
2669 }
2670 }
2671}
2672
2673template <int dim>
2676 const unsigned int codim) const
2677{
2678 Assert(codim <= dim, ExcImpossibleInDim(dim));
2679 (void)codim;
2680
2681 // vertex/line/face/cell domination
2682 // --------------------------------
2683 if (const FE_Nedelec<dim> *fe_nedelec_other =
2684 dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2685 {
2686 if (this->degree < fe_nedelec_other->degree)
2688 else if (this->degree == fe_nedelec_other->degree)
2690 else
2692 }
2693 else if (const FE_Nothing<dim> *fe_nothing =
2694 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
2695 {
2696 if (fe_nothing->is_dominating())
2698 else
2699 // the FE_Nothing has no degrees of freedom and it is typically used
2700 // in a context where we don't require any continuity along the
2701 // interface
2703 }
2704
2707}
2708
2709template <int dim>
2710bool
2712{
2713 return true;
2714}
2715
2716template <int dim>
2717std::vector<std::pair<unsigned int, unsigned int>>
2719{
2720 // Nedelec elements do not have any dofs
2721 // on vertices, hence return an empty vector.
2722 return std::vector<std::pair<unsigned int, unsigned int>>();
2723}
2724
2725template <int dim>
2726std::vector<std::pair<unsigned int, unsigned int>>
2728 const FiniteElement<dim> &fe_other) const
2729{
2730 // we can presently only compute these
2731 // identities if both FEs are
2732 // FE_Nedelec or if the other one is an
2733 // FE_Nothing
2734 if (const FE_Nedelec<dim> *fe_nedelec_other =
2735 dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2736 {
2737 // dofs are located on lines, so
2738 // two dofs are identical, if their
2739 // edge shape functions have the
2740 // same polynomial degree.
2741 std::vector<std::pair<unsigned int, unsigned int>> identities;
2742
2743 identities.reserve(std::min(fe_nedelec_other->degree, this->degree));
2744 for (unsigned int i = 0;
2745 i < std::min(fe_nedelec_other->degree, this->degree);
2746 ++i)
2747 identities.emplace_back(i, i);
2748
2749 return identities;
2750 }
2751
2752 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2753 {
2754 // the FE_Nothing has no
2755 // degrees of freedom, so there
2756 // are no equivalencies to be
2757 // recorded
2758 return std::vector<std::pair<unsigned int, unsigned int>>();
2759 }
2760
2761 else
2762 {
2764 return std::vector<std::pair<unsigned int, unsigned int>>();
2765 }
2766}
2767
2768template <int dim>
2769std::vector<std::pair<unsigned int, unsigned int>>
2771 const unsigned int) const
2772{
2773 // we can presently only compute
2774 // these identities if both FEs are
2775 // FE_Nedelec or if the other one is an
2776 // FE_Nothing
2777 if (const FE_Nedelec<dim> *fe_nedelec_other =
2778 dynamic_cast<const FE_Nedelec<dim> *>(&fe_other))
2779 {
2780 // dofs are located on the interior
2781 // of faces, so two dofs are identical,
2782 // if their face shape functions have
2783 // the same polynomial degree.
2784 const unsigned int p = fe_nedelec_other->degree;
2785 const unsigned int q = this->degree;
2786 const unsigned int p_min = std::min(p, q);
2787 std::vector<std::pair<unsigned int, unsigned int>> identities;
2788
2789 for (unsigned int i = 0; i < p_min; ++i)
2790 for (unsigned int j = 0; j < p_min - 1; ++j)
2791 {
2792 identities.emplace_back(i * (q - 1) + j, i * (p - 1) + j);
2793 identities.emplace_back(i + (j + q - 1) * q, i + (j + p - 1) * p);
2794 }
2795
2796 return identities;
2797 }
2798
2799 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
2800 {
2801 // the FE_Nothing has no
2802 // degrees of freedom, so there
2803 // are no equivalencies to be
2804 // recorded
2805 return std::vector<std::pair<unsigned int, unsigned int>>();
2806 }
2807
2808 else
2809 {
2811 return std::vector<std::pair<unsigned int, unsigned int>>();
2812 }
2813}
2814
2815// In this function we compute the face
2816// interpolation matrix. This is usually
2817// done by projection-based interpolation,
2818// but, since one can compute the entries
2819// easy per hand, we save some computation
2820// time at this point and just fill in the
2821// correct values.
2822template <int dim>
2823void
2825 const FiniteElement<dim> &source,
2826 FullMatrix<double> &interpolation_matrix,
2827 const unsigned int face_no) const
2828{
2829 (void)face_no;
2830 // this is only implemented, if the
2831 // source FE is also a
2832 // Nedelec element
2833 AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2834 (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2836 Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
2837 ExcDimensionMismatch(interpolation_matrix.m(),
2838 source.n_dofs_per_face(face_no)));
2839 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
2840 ExcDimensionMismatch(interpolation_matrix.n(),
2841 this->n_dofs_per_face(face_no)));
2842
2843 // ok, source is a Nedelec element, so
2844 // we will be able to do the work
2845 const FE_Nedelec<dim> &source_fe =
2846 dynamic_cast<const FE_Nedelec<dim> &>(source);
2847
2848 // Make sure, that the element,
2849 // for which the DoFs should be
2850 // constrained is the one with
2851 // the higher polynomial degree.
2852 // Actually the procedure will work
2853 // also if this assertion is not
2854 // satisfied. But the matrices
2855 // produced in that case might
2856 // lead to problems in the
2857 // hp-procedures, which use this
2858 // method.
2859 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
2861 interpolation_matrix = 0;
2862
2863 // On lines we can just identify
2864 // all degrees of freedom.
2865 for (unsigned int i = 0; i < this->degree; ++i)
2866 interpolation_matrix(i, i) = 1.0;
2867
2868 // In 3d we have some lines more
2869 // and a face. The procedure stays
2870 // the same as above, but we have
2871 // to take a bit more care of the
2872 // indices of the degrees of
2873 // freedom.
2874 if (dim == 3)
2875 {
2876 const unsigned int p = source_fe.degree;
2877 const unsigned int q = this->degree;
2878
2879 for (unsigned int i = 0; i < q; ++i)
2880 {
2881 for (unsigned int j = 1; j < GeometryInfo<dim>::lines_per_face; ++j)
2882 interpolation_matrix(j * p + i, j * q + i) = 1.0;
2883
2884 for (unsigned int j = 0; j < q - 1; ++j)
2885 {
2886 interpolation_matrix(GeometryInfo<dim>::lines_per_face * p +
2887 i * (p - 1) + j,
2889 i * (q - 1) + j) = 1.0;
2890 interpolation_matrix(GeometryInfo<dim>::lines_per_face * p + i +
2891 (j + p - 1) * p,
2893 (j + q - 1) * q) = 1.0;
2894 }
2895 }
2896 }
2897}
2898
2899
2900
2901template <>
2902void
2904 const unsigned int,
2906 const unsigned int) const
2907{
2909}
2910
2911
2912
2913// In this function we compute the
2914// subface interpolation matrix.
2915// This is done by a projection-
2916// based interpolation. Therefore
2917// we first interpolate the
2918// shape functions of the higher
2919// order element on the lowest
2920// order edge shape functions.
2921// Then the remaining part of
2922// the interpolated shape
2923// functions is projected on the
2924// higher order edge shape
2925// functions, the face shape
2926// functions and the interior
2927// shape functions (if they all
2928// exist).
2929template <int dim>
2930void
2932 const FiniteElement<dim> &source,
2933 const unsigned int subface,
2934 FullMatrix<double> &interpolation_matrix,
2935 const unsigned int face_no) const
2936{
2937 // this is only implemented, if the
2938 // source FE is also a
2939 // Nedelec element
2940 AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
2941 (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
2943 Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
2944 ExcDimensionMismatch(interpolation_matrix.m(),
2945 source.n_dofs_per_face(face_no)));
2946 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
2947 ExcDimensionMismatch(interpolation_matrix.n(),
2948 this->n_dofs_per_face(face_no)));
2949
2950 // ok, source is a Nedelec element, so
2951 // we will be able to do the work
2952 const FE_Nedelec<dim> &source_fe =
2953 dynamic_cast<const FE_Nedelec<dim> &>(source);
2954
2955 // Make sure, that the element,
2956 // for which the DoFs should be
2957 // constrained is the one with
2958 // the higher polynomial degree.
2959 // Actually the procedure will work
2960 // also if this assertion is not
2961 // satisfied. But the matrices
2962 // produced in that case might
2963 // lead to problems in the
2964 // hp-procedures, which use this
2965 // method.
2966 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
2968 interpolation_matrix = 0.0;
2969 // Perform projection-based interpolation
2970 // as usual.
2971 const QGauss<1> edge_quadrature(source_fe.degree);
2972 const std::vector<Point<1>> &edge_quadrature_points =
2973 edge_quadrature.get_points();
2974 const unsigned int n_edge_quadrature_points = edge_quadrature.size();
2975
2976 switch (dim)
2977 {
2978 case 2:
2979 {
2980 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
2981 ++dof)
2982 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
2983 ++q_point)
2984 {
2985 const Point<dim> quadrature_point(
2986 0.0, 0.5 * (edge_quadrature_points[q_point][0] + subface));
2987
2988 interpolation_matrix(0, dof) +=
2989 0.5 * edge_quadrature.weight(q_point) *
2990 this->shape_value_component(dof, quadrature_point, 1);
2991 }
2992
2993 if (source_fe.degree > 1)
2994 {
2995 const std::vector<Polynomials::Polynomial<double>>
2996 &legendre_polynomials =
2998 source_fe.degree - 1);
2999 FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
3000 source_fe.degree - 1);
3001
3002 {
3003 FullMatrix<double> assembling_matrix(source_fe.degree - 1,
3004 n_edge_quadrature_points);
3005
3006 for (unsigned int q_point = 0;
3007 q_point < n_edge_quadrature_points;
3008 ++q_point)
3009 {
3010 const double weight =
3011 std::sqrt(edge_quadrature.weight(q_point));
3012
3013 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
3014 assembling_matrix(i, q_point) =
3015 weight * legendre_polynomials[i + 1].value(
3016 edge_quadrature_points[q_point][0]);
3017 }
3018
3019 FullMatrix<double> system_matrix(source_fe.degree - 1,
3020 source_fe.degree - 1);
3021
3022 assembling_matrix.mTmult(system_matrix, assembling_matrix);
3023 system_matrix_inv.invert(system_matrix);
3024 }
3025
3026 Vector<double> solution(source_fe.degree - 1);
3027 Vector<double> system_rhs(source_fe.degree - 1);
3028
3029 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
3030 ++dof)
3031 {
3032 system_rhs = 0.0;
3033
3034 for (unsigned int q_point = 0;
3035 q_point < n_edge_quadrature_points;
3036 ++q_point)
3037 {
3038 const Point<dim> quadrature_point_0(
3039 0.0,
3040 0.5 * (edge_quadrature_points[q_point][0] + subface));
3041 const Point<dim> quadrature_point_1(
3042 0.0, edge_quadrature_points[q_point][0]);
3043 const double tmp =
3044 edge_quadrature.weight(q_point) *
3045 (0.5 * this->shape_value_component(dof,
3046 quadrature_point_0,
3047 1) -
3048 interpolation_matrix(0, dof) *
3049 source_fe.shape_value_component(0,
3050 quadrature_point_1,
3051 1));
3052
3053 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
3054 system_rhs(i) +=
3055 tmp * legendre_polynomials[i + 1].value(
3056 edge_quadrature_points[q_point][0]);
3057 }
3058
3059 system_matrix_inv.vmult(solution, system_rhs);
3060
3061 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
3062 if (std::abs(solution(i)) > 1e-14)
3063 interpolation_matrix(i + 1, dof) = solution(i);
3064 }
3065 }
3066
3067 break;
3068 }
3069
3070 case 3:
3071 {
3072 const double shifts[4][2] = {{0.0, 0.0},
3073 {1.0, 0.0},
3074 {0.0, 1.0},
3075 {1.0, 1.0}};
3076
3077 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
3078 ++dof)
3079 for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
3080 ++q_point)
3081 {
3082 const double weight = 0.5 * edge_quadrature.weight(q_point);
3083
3084 for (unsigned int i = 0; i < 2; ++i)
3085 {
3086 Point<dim> quadrature_point(
3087 0.5 * (i + shifts[subface][0]),
3088 0.5 * (edge_quadrature_points[q_point][0] +
3089 shifts[subface][1]),
3090 0.0);
3091
3092 interpolation_matrix(i * source_fe.degree, dof) +=
3093 weight *
3094 this->shape_value_component(
3095 this->face_to_cell_index(dof, 4), quadrature_point, 1);
3096 quadrature_point =
3097 Point<dim>(0.5 * (edge_quadrature_points[q_point][0] +
3098 shifts[subface][0]),
3099 0.5 * (i + shifts[subface][1]),
3100 0.0);
3101 interpolation_matrix((i + 2) * source_fe.degree, dof) +=
3102 weight *
3103 this->shape_value_component(
3104 this->face_to_cell_index(dof, 4), quadrature_point, 0);
3105 }
3106 }
3107
3108 if (source_fe.degree > 1)
3109 {
3110 const std::vector<Polynomials::Polynomial<double>>
3111 &legendre_polynomials =
3113 source_fe.degree - 1);
3114 FullMatrix<double> system_matrix_inv(source_fe.degree - 1,
3115 source_fe.degree - 1);
3116
3117 {
3118 FullMatrix<double> assembling_matrix(source_fe.degree - 1,
3119 n_edge_quadrature_points);
3120
3121 for (unsigned int q_point = 0;
3122 q_point < n_edge_quadrature_points;
3123 ++q_point)
3124 {
3125 const double weight =
3126 std::sqrt(edge_quadrature.weight(q_point));
3127
3128 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
3129 assembling_matrix(i, q_point) =
3130 weight * legendre_polynomials[i + 1].value(
3131 edge_quadrature_points[q_point][0]);
3132 }
3133
3134 FullMatrix<double> system_matrix(source_fe.degree - 1,
3135 source_fe.degree - 1);
3136
3137 assembling_matrix.mTmult(system_matrix, assembling_matrix);
3138 system_matrix_inv.invert(system_matrix);
3139 }
3140
3141 FullMatrix<double> solution(source_fe.degree - 1,
3143 FullMatrix<double> system_rhs(source_fe.degree - 1,
3146
3147 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
3148 ++dof)
3149 {
3150 system_rhs = 0.0;
3151
3152 for (unsigned int q_point = 0;
3153 q_point < n_edge_quadrature_points;
3154 ++q_point)
3155 {
3156 const double weight = edge_quadrature.weight(q_point);
3157
3158 for (unsigned int i = 0; i < 2; ++i)
3159 {
3160 Point<dim> quadrature_point_0(
3161 0.5 * (i + shifts[subface][0]),
3162 0.5 * (edge_quadrature_points[q_point][0] +
3163 shifts[subface][1]),
3164 0.0);
3165 Point<dim> quadrature_point_1(
3166 i, edge_quadrature_points[q_point][0], 0.0);
3167
3168 tmp(i) =
3169 weight *
3170 (0.5 * this->shape_value_component(
3171 this->face_to_cell_index(dof, 4),
3172 quadrature_point_0,
3173 1) -
3174 interpolation_matrix(i * source_fe.degree, dof) *
3175 source_fe.shape_value_component(
3176 i * source_fe.degree, quadrature_point_1, 1));
3177 quadrature_point_0 =
3178 Point<dim>(0.5 *
3179 (edge_quadrature_points[q_point][0] +
3180 shifts[subface][0]),
3181 0.5 * (i + shifts[subface][1]),
3182 0.0);
3183 quadrature_point_1 =
3184 Point<dim>(edge_quadrature_points[q_point][0],
3185 i,
3186 0.0);
3187 tmp(i + 2) =
3188 weight *
3189 (0.5 * this->shape_value_component(
3190 this->face_to_cell_index(dof, 4),
3191 quadrature_point_0,
3192 0) -
3193 interpolation_matrix((i + 2) * source_fe.degree,
3194 dof) *
3195 source_fe.shape_value_component(
3196 (i + 2) * source_fe.degree,
3197 quadrature_point_1,
3198 0));
3199 }
3200
3201 for (unsigned int i = 0; i < source_fe.degree - 1; ++i)
3202 {
3203 const double L_i = legendre_polynomials[i + 1].value(
3204 edge_quadrature_points[q_point][0]);
3205
3206 for (unsigned int j = 0;
3207 j < GeometryInfo<dim>::lines_per_face;
3208 ++j)
3209 system_rhs(i, j) += tmp(j) * L_i;
3210 }
3211 }
3212
3213 system_matrix_inv.mmult(solution, system_rhs);
3214
3215 for (unsigned int i = 0;
3216 i < GeometryInfo<dim>::lines_per_face;
3217 ++i)
3218 for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
3219 if (std::abs(solution(j, i)) > 1e-14)
3220 interpolation_matrix(i * source_fe.degree + j + 1,
3221 dof) = solution(j, i);
3222 }
3223
3224 const QGauss<2> quadrature(source_fe.degree);
3225 const std::vector<Point<2>> &quadrature_points =
3226 quadrature.get_points();
3227 const std::vector<Polynomials::Polynomial<double>>
3228 &lobatto_polynomials =
3230 source_fe.degree);
3231 const unsigned int n_boundary_dofs =
3233 const unsigned int n_quadrature_points = quadrature.size();
3234
3235 {
3236 FullMatrix<double> assembling_matrix(source_fe.degree *
3237 (source_fe.degree - 1),
3238 n_quadrature_points);
3239
3240 for (unsigned int q_point = 0; q_point < n_quadrature_points;
3241 ++q_point)
3242 {
3243 const double weight = std::sqrt(quadrature.weight(q_point));
3244
3245 for (unsigned int i = 0; i < source_fe.degree; ++i)
3246 {
3247 const double L_i =
3248 weight * legendre_polynomials[i].value(
3249 quadrature_points[q_point][0]);
3250
3251 for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
3252 assembling_matrix(i * (source_fe.degree - 1) + j,
3253 q_point) =
3254 L_i * lobatto_polynomials[j + 2].value(
3255 quadrature_points[q_point][1]);
3256 }
3257 }
3258
3259 FullMatrix<double> system_matrix(assembling_matrix.m(),
3260 assembling_matrix.m());
3261
3262 assembling_matrix.mTmult(system_matrix, assembling_matrix);
3263 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3264 system_matrix_inv.invert(system_matrix);
3265 }
3266
3267 solution.reinit(system_matrix_inv.m(), 2);
3268 system_rhs.reinit(system_matrix_inv.m(), 2);
3269 tmp.reinit(2);
3270
3271 for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
3272 ++dof)
3273 {
3274 system_rhs = 0.0;
3275
3276 for (unsigned int q_point = 0; q_point < n_quadrature_points;
3277 ++q_point)
3278 {
3279 Point<dim> quadrature_point(
3280 0.5 *
3281 (quadrature_points[q_point][0] + shifts[subface][0]),
3282 0.5 *
3283 (quadrature_points[q_point][1] + shifts[subface][1]),
3284 0.0);
3285 tmp(0) = 0.5 * this->shape_value_component(
3286 this->face_to_cell_index(dof, 4),
3287 quadrature_point,
3288 0);
3289 tmp(1) = 0.5 * this->shape_value_component(
3290 this->face_to_cell_index(dof, 4),
3291 quadrature_point,
3292 1);
3293 quadrature_point =
3294 Point<dim>(quadrature_points[q_point][0],
3295 quadrature_points[q_point][1],
3296 0.0);
3297
3298 for (unsigned int i = 0; i < 2; ++i)
3299 for (unsigned int j = 0; j < source_fe.degree; ++j)
3300 {
3301 tmp(0) -= interpolation_matrix(
3302 (i + 2) * source_fe.degree + j, dof) *
3303 source_fe.shape_value_component(
3304 (i + 2) * source_fe.degree + j,
3305 quadrature_point,
3306 0);
3307 tmp(1) -=
3308 interpolation_matrix(i * source_fe.degree + j,
3309 dof) *
3310 source_fe.shape_value_component(
3311 i * source_fe.degree + j, quadrature_point, 1);
3312 }
3313
3314 tmp *= quadrature.weight(q_point);
3315
3316 for (unsigned int i = 0; i < source_fe.degree; ++i)
3317 {
3318 const double L_i_0 = legendre_polynomials[i].value(
3319 quadrature_points[q_point][0]);
3320 const double L_i_1 = legendre_polynomials[i].value(
3321 quadrature_points[q_point][1]);
3322
3323 for (unsigned int j = 0; j < source_fe.degree - 1;
3324 ++j)
3325 {
3326 system_rhs(i * (source_fe.degree - 1) + j, 0) +=
3327 tmp(0) * L_i_0 *
3328 lobatto_polynomials[j + 2].value(
3329 quadrature_points[q_point][1]);
3330 system_rhs(i * (source_fe.degree - 1) + j, 1) +=
3331 tmp(1) * L_i_1 *
3332 lobatto_polynomials[j + 2].value(
3333 quadrature_points[q_point][0]);
3334 }
3335 }
3336 }
3337
3338 system_matrix_inv.mmult(solution, system_rhs);
3339
3340 for (unsigned int i = 0; i < source_fe.degree; ++i)
3341 for (unsigned int j = 0; j < source_fe.degree - 1; ++j)
3342 {
3343 if (std::abs(solution(i * (source_fe.degree - 1) + j,
3344 0)) > 1e-14)
3345 interpolation_matrix(i * (source_fe.degree - 1) + j +
3346 n_boundary_dofs,
3347 dof) =
3348 solution(i * (source_fe.degree - 1) + j, 0);
3349
3350 if (std::abs(solution(i * (source_fe.degree - 1) + j,
3351 1)) > 1e-14)
3352 interpolation_matrix(
3353 i + (j + source_fe.degree - 1) * source_fe.degree +
3354 n_boundary_dofs,
3355 dof) = solution(i * (source_fe.degree - 1) + j, 1);
3356 }
3357 }
3358 }
3359
3360 break;
3361 }
3362
3363 default:
3365 }
3366}
3367
3368template <int dim>
3369const FullMatrix<double> &
3371 const unsigned int child,
3372 const RefinementCase<dim> &refinement_case) const
3373{
3374 AssertIndexRange(refinement_case,
3376 Assert(refinement_case != RefinementCase<dim>::no_refinement,
3377 ExcMessage(
3378 "Prolongation matrices are only available for refined cells!"));
3380 child, this->reference_cell().template n_children<dim>(refinement_case));
3381
3382 // initialization upon first request
3383 if (this->prolongation[refinement_case - 1][child].n() == 0)
3384 {
3385 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
3386
3387 // if matrix got updated while waiting for the lock
3388 if (this->prolongation[refinement_case - 1][child].n() ==
3389 this->n_dofs_per_cell())
3390 return this->prolongation[refinement_case - 1][child];
3391
3392 // now do the work. need to get a non-const version of data in order to
3393 // be able to modify them inside a const function
3394 FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
3395
3396 // Reinit the vectors of
3397 // restriction and prolongation
3398 // matrices to the right sizes.
3399 // Restriction only for isotropic
3400 // refinement
3401#ifdef DEBUG_NEDELEC
3402 deallog << "Embedding" << std::endl;
3403#endif
3405 // Fill prolongation matrices with embedding operators
3407 this_nonconst,
3408 this_nonconst.prolongation,
3409 true,
3410 internal::FE_Nedelec::get_embedding_computation_tolerance(
3411 this->degree));
3412#ifdef DEBUG_NEDELEC
3413 deallog << "Restriction" << std::endl;
3414#endif
3415 this_nonconst.initialize_restriction();
3416 }
3417
3418 // we use refinement_case-1 here. the -1 takes care of the origin of the
3419 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3420 // available and so the vector indices are shifted
3421 return this->prolongation[refinement_case - 1][child];
3422}
3423
3424template <int dim>
3425const FullMatrix<double> &
3427 const unsigned int child,
3428 const RefinementCase<dim> &refinement_case) const
3429{
3430 AssertIndexRange(refinement_case,
3432 Assert(refinement_case != RefinementCase<dim>::no_refinement,
3433 ExcMessage(
3434 "Restriction matrices are only available for refined cells!"));
3435 AssertIndexRange(child,
3436 this->reference_cell().template n_children<dim>(
3437 RefinementCase<dim>(refinement_case)));
3438
3439 // initialization upon first request
3440 if (this->restriction[refinement_case - 1][child].n() == 0)
3441 {
3442 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
3443
3444 // if matrix got updated while waiting for the lock...
3445 if (this->restriction[refinement_case - 1][child].n() ==
3446 this->n_dofs_per_cell())
3447 return this->restriction[refinement_case - 1][child];
3448
3449 // now do the work. need to get a non-const version of data in order to
3450 // be able to modify them inside a const function
3451 FE_Nedelec<dim> &this_nonconst = const_cast<FE_Nedelec<dim> &>(*this);
3452
3453 // Reinit the vectors of
3454 // restriction and prolongation
3455 // matrices to the right sizes.
3456 // Restriction only for isotropic
3457 // refinement
3458#ifdef DEBUG_NEDELEC
3459 deallog << "Embedding" << std::endl;
3460#endif
3462 // Fill prolongation matrices with embedding operators
3464 this_nonconst,
3465 this_nonconst.prolongation,
3466 true,
3467 internal::FE_Nedelec::get_embedding_computation_tolerance(
3468 this->degree));
3469#ifdef DEBUG_NEDELEC
3470 deallog << "Restriction" << std::endl;
3471#endif
3472 this_nonconst.initialize_restriction();
3473 }
3474
3475 // we use refinement_case-1 here. the -1 takes care of the origin of the
3476 // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
3477 // available and so the vector indices are shifted
3478 return this->restriction[refinement_case - 1][child];
3479}
3480
3481
3482// Interpolate a function, which is given by
3483// its values at the generalized support
3484// points in the finite element space on the
3485// reference cell.
3486// This is done as usual by projection-based
3487// interpolation.
3488template <int dim>
3489void
3491 const std::vector<Vector<double>> &support_point_values,
3492 std::vector<double> &nodal_values) const
3493{
3494 // TODO: the implementation makes the assumption that all faces have the
3495 // same number of dofs
3496 AssertDimension(this->n_unique_faces(), 1);
3497 const unsigned int face_no = 0;
3498
3499 const unsigned int deg = this->degree - 1;
3500 Assert(support_point_values.size() == this->generalized_support_points.size(),
3501 ExcDimensionMismatch(support_point_values.size(),
3502 this->generalized_support_points.size()));
3503 Assert(support_point_values[0].size() == this->n_components(),
3504 ExcDimensionMismatch(support_point_values[0].size(),
3505 this->n_components()));
3506 Assert(nodal_values.size() == this->n_dofs_per_cell(),
3507 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
3508 std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3509
3510 switch (dim)
3511 {
3512 case 2:
3513 {
3514 // Let us begin with the
3515 // interpolation part.
3516 const QGauss<1> reference_edge_quadrature(this->degree);
3517 const unsigned int n_edge_points = reference_edge_quadrature.size();
3518
3519 for (unsigned int i = 0; i < 2; ++i)
3520 for (unsigned int j = 0; j < 2; ++j)
3521 {
3522 for (unsigned int q_point = 0; q_point < n_edge_points;
3523 ++q_point)
3524 nodal_values[(i + 2 * j) * this->degree] +=
3525 reference_edge_quadrature.weight(q_point) *
3526 support_point_values[q_point + (i + 2 * j) * n_edge_points]
3527 [1 - j];
3528
3529 // Add the computed support_point_values to the resulting vector
3530 // only, if they are not
3531 // too small.
3532 if (std::abs(nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3533 nodal_values[(i + 2 * j) * this->degree] = 0.0;
3534 }
3535
3536 // If the Nedelec element degree is greater
3537 // than 0 (i.e., the polynomial degree is greater than 1),
3538 // then we have still some higher order edge
3539 // shape functions to consider.
3540 // Note that this->degree returns the polynomial
3541 // degree.
3542 // Here the projection part starts.
3543 // The dof support_point_values are obtained by solving
3544 // a linear system of
3545 // equations.
3546 if (this->degree > 1)
3547 {
3548 // We start with projection
3549 // on the higher order edge
3550 // shape function.
3551 const std::vector<Polynomials::Polynomial<double>>
3552 &lobatto_polynomials =
3554 FullMatrix<double> system_matrix(this->degree - 1,
3555 this->degree - 1);
3556 std::vector<Polynomials::Polynomial<double>>
3557 lobatto_polynomials_grad(this->degree);
3558
3559 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3560 lobatto_polynomials_grad[i] =
3561 lobatto_polynomials[i + 1].derivative();
3562
3563 // Set up the system matrix.
3564 // This can be used for all
3565 // edges.
3566 for (unsigned int i = 0; i < system_matrix.m(); ++i)
3567 for (unsigned int j = 0; j < system_matrix.n(); ++j)
3568 for (unsigned int q_point = 0; q_point < n_edge_points;
3569 ++q_point)
3570 system_matrix(i, j) +=
3571 boundary_weights(q_point, j) *
3572 lobatto_polynomials_grad[i + 1].value(
3573 this->generalized_face_support_points[face_no][q_point]
3574 [0]);
3575
3576 FullMatrix<double> system_matrix_inv(this->degree - 1,
3577 this->degree - 1);
3578
3579 system_matrix_inv.invert(system_matrix);
3580
3581 const unsigned int
3582 line_coordinate[GeometryInfo<2>::lines_per_cell] = {1, 1, 0, 0};
3583 Vector<double> system_rhs(system_matrix.m());
3584 Vector<double> solution(system_rhs.size());
3585
3586 for (unsigned int line = 0;
3587 line < GeometryInfo<dim>::lines_per_cell;
3588 ++line)
3589 {
3590 // Set up the right hand side.
3591 system_rhs = 0;
3592
3593 for (unsigned int q_point = 0; q_point < n_edge_points;
3594 ++q_point)
3595 {
3596 const double tmp =
3597 support_point_values[line * n_edge_points + q_point]
3598 [line_coordinate[line]] -
3599 nodal_values[line * this->degree] *
3600 this->shape_value_component(
3601 line * this->degree,
3602 this->generalized_support_points[line *
3603 n_edge_points +
3604 q_point],
3605 line_coordinate[line]);
3606
3607 for (unsigned int i = 0; i < system_rhs.size(); ++i)
3608 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3609 }
3610
3611 system_matrix_inv.vmult(solution, system_rhs);
3612
3613 // Add the computed support_point_values
3614 // to the resulting vector
3615 // only, if they are not
3616 // too small.
3617 for (unsigned int i = 0; i < solution.size(); ++i)
3618 if (std::abs(solution(i)) > 1e-14)
3619 nodal_values[line * this->degree + i + 1] = solution(i);
3620 }
3621
3622 // Then we go on to the
3623 // interior shape
3624 // functions. Again we
3625 // set up the system
3626 // matrix and use it
3627 // for both, the
3628 // horizontal and the
3629 // vertical, interior
3630 // shape functions.
3631 const QGauss<dim> reference_quadrature(this->degree);
3632 const unsigned int n_interior_points =
3633 reference_quadrature.size();
3634 const std::vector<Polynomials::Polynomial<double>>
3635 &legendre_polynomials =
3637 1);
3638
3639 system_matrix.reinit((this->degree - 1) * this->degree,
3640 (this->degree - 1) * this->degree);
3641 system_matrix = 0;
3642
3643 for (unsigned int i = 0; i < this->degree; ++i)
3644 for (unsigned int j = 0; j < this->degree - 1; ++j)
3645 for (unsigned int k = 0; k < this->degree; ++k)
3646 for (unsigned int l = 0; l < this->degree - 1; ++l)
3647 for (unsigned int q_point = 0;
3648 q_point < n_interior_points;
3649 ++q_point)
3650 system_matrix(i * (this->degree - 1) + j,
3651 k * (this->degree - 1) + l) +=
3652 reference_quadrature.weight(q_point) *
3653 legendre_polynomials[i].value(
3654 this->generalized_support_points
3656 n_edge_points][0]) *
3657 lobatto_polynomials[j + 2].value(
3658 this->generalized_support_points
3660 n_edge_points][1]) *
3661 lobatto_polynomials_grad[k].value(
3662 this->generalized_support_points
3664 n_edge_points][0]) *
3665 lobatto_polynomials[l + 2].value(
3666 this->generalized_support_points
3668 n_edge_points][1]);
3669
3670 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3671 system_matrix_inv.invert(system_matrix);
3672 // Set up the right hand side
3673 // for the horizontal shape
3674 // functions.
3675 system_rhs.reinit(system_matrix_inv.m());
3676 system_rhs = 0;
3677
3678 for (unsigned int q_point = 0; q_point < n_interior_points;
3679 ++q_point)
3680 {
3681 double tmp =
3682 support_point_values[q_point +
3684 n_edge_points][0];
3685
3686 for (unsigned int i = 0; i < 2; ++i)
3687 for (unsigned int j = 0; j <= deg; ++j)
3688 tmp -= nodal_values[(i + 2) * this->degree + j] *
3689 this->shape_value_component(
3690 (i + 2) * this->degree + j,
3691 this->generalized_support_points
3693 n_edge_points],
3694 0);
3695
3696 for (unsigned int i = 0; i <= deg; ++i)
3697 for (unsigned int j = 0; j < deg; ++j)
3698 system_rhs(i * deg + j) +=
3699 reference_quadrature.weight(q_point) * tmp *
3700 lobatto_polynomials_grad[i].value(
3701 this->generalized_support_points
3703 n_edge_points][0]) *
3704 lobatto_polynomials[j + 2].value(
3705 this->generalized_support_points
3707 n_edge_points][1]);
3708 }
3709
3710 solution.reinit(system_matrix.m());
3711 system_matrix_inv.vmult(solution, system_rhs);
3712
3713 // Add the computed support_point_values
3714 // to the resulting vector
3715 // only, if they are not
3716 // too small.
3717 for (unsigned int i = 0; i <= deg; ++i)
3718 for (unsigned int j = 0; j < deg; ++j)
3719 if (std::abs(solution(i * deg + j)) > 1e-14)
3720 nodal_values[(i + GeometryInfo<dim>::lines_per_cell) * deg +
3722 solution(i * deg + j);
3723
3724 system_rhs = 0;
3725 // Set up the right hand side
3726 // for the vertical shape
3727 // functions.
3728
3729 for (unsigned int q_point = 0; q_point < n_interior_points;
3730 ++q_point)
3731 {
3732 double tmp =
3733 support_point_values[q_point +
3735 n_edge_points][1];
3736
3737 for (unsigned int i = 0; i < 2; ++i)
3738 for (unsigned int j = 0; j <= deg; ++j)
3739 tmp -= nodal_values[i * this->degree + j] *
3740 this->shape_value_component(
3741 i * this->degree + j,
3742 this->generalized_support_points
3744 n_edge_points],
3745 1);
3746
3747 for (unsigned int i = 0; i <= deg; ++i)
3748 for (unsigned int j = 0; j < deg; ++j)
3749 system_rhs(i * deg + j) +=
3750 reference_quadrature.weight(q_point) * tmp *
3751 lobatto_polynomials_grad[i].value(
3752 this->generalized_support_points
3754 n_edge_points][1]) *
3755 lobatto_polynomials[j + 2].value(
3756 this->generalized_support_points
3758 n_edge_points][0]);
3759 }
3760
3761 system_matrix_inv.vmult(solution, system_rhs);
3762
3763 // Add the computed support_point_values
3764 // to the resulting vector
3765 // only, if they are not
3766 // too small.
3767 for (unsigned int i = 0; i <= deg; ++i)
3768 for (unsigned int j = 0; j < deg; ++j)
3769 if (std::abs(solution(i * deg + j)) > 1e-14)
3770 nodal_values[i +
3772 this->degree] = solution(i * deg + j);
3773 }
3774
3775 break;
3776 }
3777
3778 case 3:
3779 {
3780 // Let us begin with the
3781 // interpolation part.
3782 const QGauss<1> reference_edge_quadrature(this->degree);
3783 const unsigned int n_edge_points = reference_edge_quadrature.size();
3784
3785 for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3786 {
3787 for (unsigned int i = 0; i < 4; ++i)
3788 nodal_values[(i + 8) * this->degree] +=
3789 reference_edge_quadrature.weight(q_point) *
3790 support_point_values[q_point + (i + 8) * n_edge_points][2];
3791
3792 for (unsigned int i = 0; i < 2; ++i)
3793 for (unsigned int j = 0; j < 2; ++j)
3794 for (unsigned int k = 0; k < 2; ++k)
3795 nodal_values[(i + 2 * (2 * j + k)) * this->degree] +=
3796 reference_edge_quadrature.weight(q_point) *
3797 support_point_values[q_point + (i + 2 * (2 * j + k)) *
3798 n_edge_points][1 - k];
3799 }
3800
3801 // Add the computed support_point_values
3802 // to the resulting vector
3803 // only, if they are not
3804 // too small.
3805 for (unsigned int i = 0; i < 4; ++i)
3806 if (std::abs(nodal_values[(i + 8) * this->degree]) < 1e-14)
3807 nodal_values[(i + 8) * this->degree] = 0.0;
3808
3809 for (unsigned int i = 0; i < 2; ++i)
3810 for (unsigned int j = 0; j < 2; ++j)
3811 for (unsigned int k = 0; k < 2; ++k)
3812 if (std::abs(
3813 nodal_values[(i + 2 * (2 * j + k)) * this->degree]) <
3814 1e-14)
3815 nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3816
3817 // If the degree is greater
3818 // than 0, then we have still
3819 // some higher order shape
3820 // functions to consider.
3821 // Here the projection part
3822 // starts. The dof support_point_values
3823 // are obtained by solving
3824 // a linear system of
3825 // equations.
3826 if (this->degree > 1)
3827 {
3828 // We start with projection
3829 // on the higher order edge
3830 // shape function.
3831 const std::vector<Polynomials::Polynomial<double>>
3832 &lobatto_polynomials =
3834 FullMatrix<double> system_matrix(this->degree - 1,
3835 this->degree - 1);
3836 std::vector<Polynomials::Polynomial<double>>
3837 lobatto_polynomials_grad(this->degree);
3838
3839 for (unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3840 lobatto_polynomials_grad[i] =
3841 lobatto_polynomials[i + 1].derivative();
3842
3843 // Set up the system matrix.
3844 // This can be used for all
3845 // edges.
3846 for (unsigned int i = 0; i < system_matrix.m(); ++i)
3847 for (unsigned int j = 0; j < system_matrix.n(); ++j)
3848 for (unsigned int q_point = 0; q_point < n_edge_points;
3849 ++q_point)
3850 system_matrix(i, j) +=
3851 boundary_weights(q_point, j) *
3852 lobatto_polynomials_grad[i + 1].value(
3853 this->generalized_face_support_points[face_no][q_point]
3854 [1]);
3855
3856 FullMatrix<double> system_matrix_inv(this->degree - 1,
3857 this->degree - 1);
3858
3859 system_matrix_inv.invert(system_matrix);
3860
3861 const unsigned int
3862 line_coordinate[GeometryInfo<3>::lines_per_cell] = {
3863 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3864 Vector<double> system_rhs(system_matrix.m());
3865 Vector<double> solution(system_rhs.size());
3866
3867 for (unsigned int line = 0;
3868 line < GeometryInfo<dim>::lines_per_cell;
3869 ++line)
3870 {
3871 // Set up the right hand side.
3872 system_rhs = 0;
3873
3874 for (unsigned int q_point = 0; q_point < this->degree;
3875 ++q_point)
3876 {
3877 const double tmp =
3878 support_point_values[line * this->degree + q_point]
3879 [line_coordinate[line]] -
3880 nodal_values[line * this->degree] *
3881 this->shape_value_component(
3882 line * this->degree,
3883 this
3884 ->generalized_support_points[line * this->degree +
3885 q_point],
3886 line_coordinate[line]);
3887
3888 for (unsigned int i = 0; i < system_rhs.size(); ++i)
3889 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3890 }
3891
3892 system_matrix_inv.vmult(solution, system_rhs);
3893
3894 // Add the computed values
3895 // to the resulting vector
3896 // only, if they are not
3897 // too small.
3898 for (unsigned int i = 0; i < solution.size(); ++i)
3899 if (std::abs(solution(i)) > 1e-14)
3900 nodal_values[line * this->degree + i + 1] = solution(i);
3901 }
3902
3903 // Then we go on to the
3904 // face shape functions.
3905 // Again we set up the
3906 // system matrix and
3907 // use it for both, the
3908 // horizontal and the
3909 // vertical, shape
3910 // functions.
3911 const std::vector<Polynomials::Polynomial<double>>
3912 &legendre_polynomials =
3914 1);
3915 const unsigned int n_face_points = n_edge_points * n_edge_points;
3916
3917 system_matrix.reinit((this->degree - 1) * this->degree,
3918 (this->degree - 1) * this->degree);
3919 system_matrix = 0;
3920
3921 for (unsigned int i = 0; i < this->degree; ++i)
3922 for (unsigned int j = 0; j < this->degree - 1; ++j)
3923 for (unsigned int k = 0; k < this->degree; ++k)
3924 for (unsigned int l = 0; l < this->degree - 1; ++l)
3925 for (unsigned int q_point = 0; q_point < n_face_points;
3926 ++q_point)
3927 system_matrix(i * (this->degree - 1) + j,
3928 k * (this->degree - 1) + l) +=
3929 boundary_weights(q_point + n_edge_points,
3930 2 * (k * (this->degree - 1) + l)) *
3931 legendre_polynomials[i].value(
3932 this->generalized_face_support_points
3933 [face_no][q_point + 4 * n_edge_points][0]) *
3934 lobatto_polynomials[j + 2].value(
3935 this->generalized_face_support_points
3936 [face_no][q_point + 4 * n_edge_points][1]);
3937
3938 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
3939 system_matrix_inv.invert(system_matrix);
3940 solution.reinit(system_matrix.m());
3941 system_rhs.reinit(system_matrix.m());
3942
3943 const unsigned int
3944 face_coordinates[GeometryInfo<3>::faces_per_cell][2] = {
3945 {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3948 {{0, 4, 8, 10},
3949 {1, 5, 9, 11},
3950 {8, 9, 2, 6},
3951 {10, 11, 3, 7},
3952 {2, 3, 0, 1},
3953 {6, 7, 4, 5}};
3954
3955 for (const unsigned int face : GeometryInfo<dim>::face_indices())
3956 {
3957 // Set up the right hand side
3958 // for the horizontal shape
3959 // functions.
3960 system_rhs = 0;
3961
3962 for (unsigned int q_point = 0; q_point < n_face_points;
3963 ++q_point)
3964 {
3965 double tmp =
3966 support_point_values[q_point +
3968 n_edge_points]
3969 [face_coordinates[face][0]];
3970
3971 for (unsigned int i = 0; i < 2; ++i)
3972 for (unsigned int j = 0; j <= deg; ++j)
3973 tmp -=
3974 nodal_values[edge_indices[face][i] * this->degree +
3975 j] *
3976 this->shape_value_component(
3977 edge_indices[face][i] * this->degree + j,
3978 this->generalized_support_points
3980 n_edge_points],
3981 face_coordinates[face][0]);
3982
3983 for (unsigned int i = 0; i <= deg; ++i)
3984 for (unsigned int j = 0; j < deg; ++j)
3985 system_rhs(i * deg + j) +=
3986 boundary_weights(q_point + n_edge_points,
3987 2 * (i * deg + j)) *
3988 tmp;
3989 }
3990
3991 system_matrix_inv.vmult(solution, system_rhs);
3992
3993 // Add the computed support_point_values
3994 // to the resulting vector
3995 // only, if they are not
3996 // too small.
3997 for (unsigned int i = 0; i <= deg; ++i)
3998 for (unsigned int j = 0; j < deg; ++j)
3999 if (std::abs(solution(i * deg + j)) > 1e-14)
4000 nodal_values[(2 * face * this->degree + i +
4002 deg +
4004 solution(i * deg + j);
4005
4006 // Set up the right hand side
4007 // for the vertical shape
4008 // functions.
4009 system_rhs = 0;
4010
4011 for (unsigned int q_point = 0; q_point < n_face_points;
4012 ++q_point)
4013 {
4014 double tmp =
4015 support_point_values[q_point +
4017 n_edge_points]
4018 [face_coordinates[face][1]];
4019
4020 for (unsigned int i = 2;
4021 i < GeometryInfo<dim>::lines_per_face;
4022 ++i)
4023 for (unsigned int j = 0; j <= deg; ++j)
4024 tmp -=
4025 nodal_values[edge_indices[face][i] * this->degree +
4026 j] *
4027 this->shape_value_component(
4028 edge_indices[face][i] * this->degree + j,
4029 this->generalized_support_points
4031 n_edge_points],
4032 face_coordinates[face][1]);
4033
4034 for (unsigned int i = 0; i <= deg; ++i)
4035 for (unsigned int j = 0; j < deg; ++j)
4036 system_rhs(i * deg + j) +=
4037 boundary_weights(q_point + n_edge_points,
4038 2 * (i * deg + j) + 1) *
4039 tmp;
4040 }
4041
4042 system_matrix_inv.vmult(solution, system_rhs);
4043
4044 // Add the computed support_point_values
4045 // to the resulting vector
4046 // only, if they are not
4047 // too small.
4048 for (unsigned int i = 0; i <= deg; ++i)
4049 for (unsigned int j = 0; j < deg; ++j)
4050 if (std::abs(solution(i * deg + j)) > 1e-14)
4051 nodal_values[((2 * face + 1) * deg + j +
4053 this->degree +
4054 i] = solution(i * deg + j);
4055 }
4056
4057 // Finally we project
4058 // the remaining parts
4059 // of the function on
4060 // the interior shape
4061 // functions.
4062 const QGauss<dim> reference_quadrature(this->degree);
4063 const unsigned int n_interior_points =
4064 reference_quadrature.size();
4065
4066 // We create the
4067 // system matrix.
4068 system_matrix.reinit(this->degree * deg * deg,
4069 this->degree * deg * deg);
4070 system_matrix = 0;
4071
4072 for (unsigned int i = 0; i <= deg; ++i)
4073 for (unsigned int j = 0; j < deg; ++j)
4074 for (unsigned int k = 0; k < deg; ++k)
4075 for (unsigned int l = 0; l <= deg; ++l)
4076 for (unsigned int m = 0; m < deg; ++m)
4077 for (unsigned int n = 0; n < deg; ++n)
4078 for (unsigned int q_point = 0;
4079 q_point < n_interior_points;
4080 ++q_point)
4081 system_matrix((i * deg + j) * deg + k,
4082 (l * deg + m) * deg + n) +=
4083 reference_quadrature.weight(q_point) *
4084 legendre_polynomials[i].value(
4085 this->generalized_support_points
4086 [q_point +
4088 n_edge_points +
4090 n_face_points][0]) *
4091 lobatto_polynomials[j + 2].value(
4092 this->generalized_support_points
4093 [q_point +
4095 n_edge_points +
4097 n_face_points][1]) *
4098 lobatto_polynomials[k + 2].value(
4099 this->generalized_support_points
4100 [q_point +
4102 n_edge_points +
4104 n_face_points][2]) *
4105 lobatto_polynomials_grad[l].value(
4106 this->generalized_support_points
4107 [q_point +
4109 n_edge_points +
4111 n_face_points][0]) *
4112 lobatto_polynomials[m + 2].value(
4113 this->generalized_support_points
4114 [q_point +
4116 n_edge_points +
4118 n_face_points][1]) *
4119 lobatto_polynomials[n + 2].value(
4120 this->generalized_support_points
4121 [q_point +
4123 n_edge_points +
4125 n_face_points][2]);
4126
4127 system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
4128 system_matrix_inv.invert(system_matrix);
4129 // Set up the right hand side.
4130 system_rhs.reinit(system_matrix.m());
4131 system_rhs = 0;
4132
4133 for (unsigned int q_point = 0; q_point < n_interior_points;
4134 ++q_point)
4135 {
4136 double tmp =
4137 support_point_values[q_point +
4139 n_edge_points +
4141 n_face_points][0];
4142
4143 for (unsigned int i = 0; i <= deg; ++i)
4144 {
4145 for (unsigned int j = 0; j < 2; ++j)
4146 for (unsigned int k = 0; k < 2; ++k)
4147 tmp -=
4148 nodal_values[i + (j + 4 * k + 2) * this->degree] *
4149 this->shape_value_component(
4150 i + (j + 4 * k + 2) * this->degree,
4151 this->generalized_support_points
4152 [q_point +
4154 n_edge_points +
4156 n_face_points],
4157 0);
4158
4159 for (unsigned int j = 0; j < deg; ++j)
4160 for (unsigned int k = 0; k < 4; ++k)
4161 tmp -=
4162 nodal_values[(i + 2 * (k + 2) * this->degree +
4164 deg +
4165 j +
4167 this->shape_value_component(
4168 (i + 2 * (k + 2) * this->degree +
4170 deg +
4172 this->generalized_support_points
4173 [q_point +
4175 n_edge_points +
4177 n_face_points],
4178 0);
4179 }
4180
4181 for (unsigned int i = 0; i <= deg; ++i)
4182 for (unsigned int j = 0; j < deg; ++j)
4183 for (unsigned int k = 0; k < deg; ++k)
4184 system_rhs((i * deg + j) * deg + k) +=
4185 reference_quadrature.weight(q_point) * tmp *
4186 lobatto_polynomials_grad[i].value(
4187 this->generalized_support_points
4188 [q_point +
4190 n_edge_points +
4192 n_face_points][0]) *
4193 lobatto_polynomials[j + 2].value(
4194 this->generalized_support_points
4195 [q_point +
4197 n_edge_points +
4199 n_face_points][1]) *
4200 lobatto_polynomials[k + 2].value(
4201 this->generalized_support_points
4202 [q_point +
4204 n_edge_points +
4206 n_face_points][2]);
4207 }
4208
4209 solution.reinit(system_rhs.size());
4210 system_matrix_inv.vmult(solution, system_rhs);
4211
4212 // Add the computed values
4213 // to the resulting vector
4214 // only, if they are not
4215 // too small.
4216 for (unsigned int i = 0; i <= deg; ++i)
4217 for (unsigned int j = 0; j < deg; ++j)
4218 for (unsigned int k = 0; k < deg; ++k)
4219 if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4220 nodal_values
4221 [((i + 2 * GeometryInfo<dim>::faces_per_cell) * deg +
4224 deg +
4226 solution((i * deg + j) * deg + k);
4227
4228 // Set up the right hand side.
4229 system_rhs = 0;
4230
4231 for (unsigned int q_point = 0; q_point < n_interior_points;
4232 ++q_point)
4233 {
4234 double tmp =
4235 support_point_values[q_point +
4237 n_edge_points +
4239 n_face_points][1];
4240
4241 for (unsigned int i = 0; i <= deg; ++i)
4242 for (unsigned int j = 0; j < 2; ++j)
4243 {
4244 for (unsigned int k = 0; k < 2; ++k)
4245 tmp -= nodal_values[i + (4 * j + k) * this->degree] *
4246 this->shape_value_component(
4247 i + (4 * j + k) * this->degree,
4248 this->generalized_support_points
4249 [q_point +
4251 n_edge_points +
4253 n_face_points],
4254 1);
4255
4256 for (unsigned int k = 0; k < deg; ++k)
4257 tmp -=
4258 nodal_values[(i + 2 * j * this->degree +
4260 deg +
4261 k +
4263 this->shape_value_component(
4264 (i + 2 * j * this->degree +
4266 deg +
4268 this->generalized_support_points
4269 [q_point +
4271 n_edge_points +
4273 n_face_points],
4274 1) +
4275 nodal_values[i +
4276 ((2 * j + 9) * deg + k +
4278 this->degree] *
4279 this->shape_value_component(
4280 i + ((2 * j + 9) * deg + k +
4282 this->degree,
4283 this->generalized_support_points
4284 [q_point +
4286 n_edge_points +
4288 n_face_points],
4289 1);
4290 }
4291
4292 for (unsigned int i = 0; i <= deg; ++i)
4293 for (unsigned int j = 0; j < deg; ++j)
4294 for (unsigned int k = 0; k < deg; ++k)
4295 system_rhs((i * deg + j) * deg + k) +=
4296 reference_quadrature.weight(q_point) * tmp *
4297 lobatto_polynomials_grad[i].value(
4298 this->generalized_support_points
4299 [q_point +
4301 n_edge_points +
4303 n_face_points][1]) *
4304 lobatto_polynomials[j + 2].value(
4305 this->generalized_support_points
4306 [q_point +
4308 n_edge_points +
4310 n_face_points][0]) *
4311 lobatto_polynomials[k + 2].value(
4312 this->generalized_support_points
4313 [q_point +
4315 n_edge_points +
4317 n_face_points][2]);
4318 }
4319
4320 system_matrix_inv.vmult(solution, system_rhs);
4321
4322 // Add the computed support_point_values
4323 // to the resulting vector
4324 // only, if they are not
4325 // too small.
4326 for (unsigned int i = 0; i <= deg; ++i)
4327 for (unsigned int j = 0; j < deg; ++j)
4328 for (unsigned int k = 0; k < deg; ++k)
4329 if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4330 nodal_values[((i + this->degree +
4332 deg +
4335 deg +
4337 solution((i * deg + j) * deg + k);
4338
4339 // Set up the right hand side.
4340 system_rhs = 0;
4341
4342 for (unsigned int q_point = 0; q_point < n_interior_points;
4343 ++q_point)
4344 {
4345 double tmp =
4346 support_point_values[q_point +
4348 n_edge_points +
4350 n_face_points][2];
4351
4352 for (unsigned int i = 0; i <= deg; ++i)
4353 for (unsigned int j = 0; j < 4; ++j)
4354 {
4355 tmp -= nodal_values[i + (j + 8) * this->degree] *
4356 this->shape_value_component(
4357 i + (j + 8) * this->degree,
4358 this->generalized_support_points
4359 [q_point +
4361 n_edge_points +
4363 n_face_points],
4364 2);
4365
4366 for (unsigned int k = 0; k < deg; ++k)
4367 tmp -=
4368 nodal_values[i +
4369 ((2 * j + 1) * deg + k +
4371 this->degree] *
4372 this->shape_value_component(
4373 i + ((2 * j + 1) * deg + k +
4375 this->degree,
4376 this->generalized_support_points
4377 [q_point +
4379 n_edge_points +
4381 n_face_points],
4382 2);
4383 }
4384
4385 for (unsigned int i = 0; i <= deg; ++i)
4386 for (unsigned int j = 0; j < deg; ++j)
4387 for (unsigned int k = 0; k < deg; ++k)
4388 system_rhs((i * deg + j) * deg + k) +=
4389 reference_quadrature.weight(q_point) * tmp *
4390 lobatto_polynomials_grad[i].value(
4391 this->generalized_support_points
4392 [q_point +
4394 n_edge_points +
4396 n_face_points][2]) *
4397 lobatto_polynomials[j + 2].value(
4398 this->generalized_support_points
4399 [q_point +
4401 n_edge_points +
4403 n_face_points][0]) *
4404 lobatto_polynomials[k + 2].value(
4405 this->generalized_support_points
4406 [q_point +
4408 n_edge_points +
4410 n_face_points][1]);
4411 }
4412
4413 system_matrix_inv.vmult(solution, system_rhs);
4414
4415 // Add the computed support_point_values
4416 // to the resulting vector
4417 // only, if they are not
4418 // too small.
4419 for (unsigned int i = 0; i <= deg; ++i)
4420 for (unsigned int j = 0; j < deg; ++j)
4421 for (unsigned int k = 0; k < deg; ++k)
4422 if (std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4423 nodal_values
4424 [i +
4425 ((j + 2 * (deg + GeometryInfo<dim>::faces_per_cell)) *
4426 deg +
4428 this->degree] = solution((i * deg + j) * deg + k);
4429 }
4430
4431 break;
4432 }
4433
4434 default:
4436 }
4437}
4438
4439
4440
4441template <int dim>
4442std::pair<Table<2, bool>, std::vector<unsigned int>>
4444{
4445 Table<2, bool> constant_modes(dim, this->n_dofs_per_cell());
4446 for (unsigned int d = 0; d < dim; ++d)
4447 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
4448 constant_modes(d, i) = true;
4449 std::vector<unsigned int> components;
4450 components.reserve(dim);
4451 for (unsigned int d = 0; d < dim; ++d)
4452 components.push_back(d);
4453 return std::pair<Table<2, bool>, std::vector<unsigned int>>(constant_modes,
4454 components);
4455}
4456
4457
4458template <int dim>
4459std::size_t
4461{
4463 return 0;
4464}
4465
4466template <int dim>
4467std::vector<unsigned int>
4468FE_Nedelec<dim>::get_embedding_dofs(const unsigned int sub_degree) const
4469{
4470 Assert((sub_degree > 0) && (sub_degree <= this->degree),
4471 ExcIndexRange(sub_degree, 1, this->degree));
4472
4473 switch (dim)
4474 {
4475 case 2:
4476 {
4477 // The Nedelec cell has only Face (Line) and Cell DoFs...
4478 const unsigned int n_face_dofs_sub =
4480 const unsigned int n_cell_dofs_sub =
4481 2 * (sub_degree - 1) * sub_degree;
4482
4483 std::vector<unsigned int> embedding_dofs(n_face_dofs_sub +
4484 n_cell_dofs_sub);
4485
4486 unsigned int i = 0;
4487
4488 // Identify the Face/Line DoFs
4489 while (i < n_face_dofs_sub)
4490 {
4491 const unsigned int face_index = i / sub_degree;
4492 embedding_dofs[i] = i % sub_degree + face_index * this->degree;
4493 ++i;
4494 }
4495
4496 // Identify the Cell DoFs
4497 if (sub_degree >= 2)
4498 {
4499 const unsigned int n_face_dofs =
4500 GeometryInfo<dim>::lines_per_cell * this->degree;
4501
4502 // For the first component
4503 for (unsigned ku = 0; ku < sub_degree; ++ku)
4504 for (unsigned kv = 2; kv <= sub_degree; ++kv)
4505 embedding_dofs[i++] =
4506 n_face_dofs + ku * (this->degree - 1) + (kv - 2);
4507
4508 // For the second component
4509 for (unsigned ku = 2; ku <= sub_degree; ++ku)
4510 for (unsigned kv = 0; kv < sub_degree; ++kv)
4511 embedding_dofs[i++] = n_face_dofs +
4512 this->degree * (this->degree - 1) +
4513 (ku - 2) * (this->degree) + kv;
4514 }
4515 Assert(i == (n_face_dofs_sub + n_cell_dofs_sub), ExcInternalError());
4516 return embedding_dofs;
4517 }
4518 default:
4520 return std::vector<unsigned int>();
4521 }
4522}
4523//----------------------------------------------------------------------//
4524
4525
4526// explicit instantiations
4527#include "fe/fe_nedelec.inst"
4528
4529
void initialize_quad_dof_index_permutation_and_sign_change()
std::vector< unsigned int > get_embedding_dofs(const unsigned int sub_degree) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual bool hp_constraints_are_implemented() const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
void initialize_restriction()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree, bool dg=false)
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) 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::size_t memory_consumption() const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual std::string get_name() const override
friend class FE_Nedelec
Definition fe_nedelec.h:655
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
void initialize_support_points(const unsigned int order)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
FullMatrix< double > inverse_node_matrix
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const override
std::vector< MappingKind > mapping_kind
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_unique_faces() const
virtual std::string get_name() const =0
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
Definition fe.h:2549
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2537
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
size_type n() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void invert(const FullMatrix< number2 > &M)
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
size_type m() const
Definition point.h:113
static unsigned int n_polynomials(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:69
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
const Point< dim > & point(const unsigned int i) const
virtual size_type size() const override
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int edge_indices[GeometryInfo< dim >::lines_per_cell]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ mapping_nedelec
Definition mapping.h:131
LogStream deallog
Definition logstream.cc:38
std::size_t size
Definition mpi.cc:745
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)
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()