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