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