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