deal.II version GIT relicensing-6809-ge913b9bb34 2026-09-25 17:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor_product_matrix_creator.h
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2022 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#ifndef dealii_tensor_product_matrix_creator_h
14#define dealii_tensor_product_matrix_creator_h
15
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/quadrature_lib.h>
20
21#include <deal.II/dofs/dof_handler.h>
22
23#include <deal.II/fe/fe.h>
24#include <deal.II/fe/fe_dgq.h>
25#include <deal.II/fe/fe_q.h>
26#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/fe_values.h>
28#include <deal.II/fe/mapping_q1.h>
29
30#include <deal.II/grid/grid_generator.h>
31#include <deal.II/grid/tria.h>
32
33#include <set>
34
35DEAL_II_NAMESPACE_OPEN
36
37
43namespace TensorProductMatrixCreator
44{
54
64 template <int dim, typename Number>
65 std::pair<std::array<FullMatrix<Number>, dim>,
66 std::array<FullMatrix<Number>, dim>>
67 create_laplace_tensor_product_matrix(
68 const FiniteElement<1> &fe,
69 const Quadrature<1> &quadrature,
70 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
71 const ::ndarray<double, dim, 3> &cell_extent,
72 const unsigned int n_overlap = 1);
73
78 template <int dim, typename Number>
79 std::pair<std::array<FullMatrix<Number>, dim>,
80 std::array<FullMatrix<Number>, dim>>
81 create_laplace_tensor_product_matrix(
82 const typename Triangulation<dim>::cell_iterator &cell,
83 const std::set<types::boundary_id> &dirichlet_boundaries,
84 const std::set<types::boundary_id> &neumann_boundaries,
85 const FiniteElement<1> &fe,
86 const Quadrature<1> &quadrature,
87 const ::ndarray<double, dim, 3> &cell_extent,
88 const unsigned int n_overlap = 1);
89
90
91
114 template <typename Number = double>
115 FullMatrix<Number>
116 create_1d_cell_mass_matrix(
117 const FiniteElement<1> &fe,
118 const Number &h,
119 const std::pair<bool, bool> include_endpoints = {true, true},
120 std::vector<unsigned int> numbering = std::vector<unsigned int>());
121
122
123
142 template <typename Number = double>
143 FullMatrix<Number>
144 create_1d_cell_laplace_matrix(
145 const FiniteElement<1> &fe,
146 const Number &h,
147 const std::pair<bool, bool> include_endpoints = {true, true},
148 std::vector<unsigned int> numbering = std::vector<unsigned int>());
149
150
182 template <typename Number = double>
183 FullMatrix<Number>
184 create_1D_discretization_matrix(
185 FullMatrix<Number> &cell_matrix,
186 const unsigned int &n_cells,
187 const unsigned int &overlap,
188 const std::pair<bool, bool> include_endpoints = {true, true});
189
190
207 template <typename Number = double>
208 FullMatrix<Number>
209 create_1d_ghost_penalty_matrix(
210 const FiniteElement<1> &fe,
211 const Number h,
212 std::vector<Number> coefficients = std::vector<Number>());
213
214
215
233 template <typename Number = double>
234 FullMatrix<Number>
235 create_1d_ghost_penalty_matrix(
236 const std::vector<Polynomials::Polynomial<double>>
237 &polynomial_basis_derivative,
238 const unsigned int overlap = 1);
239
240} // namespace TensorProductMatrixCreator
241
242
243
244/*----------------------- Inline functions ----------------------------------*/
245
246
247namespace TensorProductMatrixCreator
248{
249 namespace internal
250 {
251 template <typename Number>
252 void
253 clear_row_and_column(const unsigned int n_dofs_1D_with_overlap,
254 const unsigned int n,
255 FullMatrix<Number> &matrix)
256 {
257 for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i)
258 {
259 matrix[i][n] = 0.0;
260 matrix[n][i] = 0.0;
261 }
262 }
263
264 template <typename Number>
265 std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>
266 create_reference_mass_and_stiffness_matrices(
267 const FiniteElement<1> &fe,
268 const Quadrature<1> &quadrature)
269 {
270 Triangulation<1> tria;
271 GridGenerator::hyper_cube(tria);
272
273 DoFHandler<1> dof_handler(tria);
274 dof_handler.distribute_dofs(fe);
275
276 MappingQ1<1> mapping;
277
278 const unsigned int n_dofs_1D = fe.n_dofs_per_cell();
279
280 FullMatrix<Number> mass_matrix_reference(n_dofs_1D, n_dofs_1D);
281 FullMatrix<Number> derivative_matrix_reference(n_dofs_1D, n_dofs_1D);
282
283 FEValues<1> fe_values(mapping,
284 fe,
285 quadrature,
286 update_values | update_gradients |
287 update_JxW_values);
288
289 fe_values.reinit(tria.begin());
290
291 const auto lexicographic_to_hierarchic_numbering =
292 Utilities::invert_permutation(
293 FETools::hierarchic_to_lexicographic_numbering<1>(
294 fe.tensor_degree()));
295
296 for (const unsigned int q_index : fe_values.quadrature_point_indices())
297 for (const unsigned int i : fe_values.dof_indices())
298 for (const unsigned int j : fe_values.dof_indices())
299 {
300 mass_matrix_reference(i, j) +=
301 (fe_values.shape_value(lexicographic_to_hierarchic_numbering[i],
302 q_index) *
303 fe_values.shape_value(lexicographic_to_hierarchic_numbering[j],
304 q_index) *
305 fe_values.JxW(q_index));
306
307 derivative_matrix_reference(i, j) +=
308 (fe_values.shape_grad(lexicographic_to_hierarchic_numbering[i],
309 q_index) *
310 fe_values.shape_grad(lexicographic_to_hierarchic_numbering[j],
311 q_index) *
312 fe_values.JxW(q_index));
313 }
314
315 return std::tuple<FullMatrix<Number>, FullMatrix<Number>, bool>{
316 mass_matrix_reference, derivative_matrix_reference, false};
317 }
318 } // namespace internal
319
320
321
322 template <int dim, typename Number>
323 std::pair<std::array<FullMatrix<Number>, dim>,
324 std::array<FullMatrix<Number>, dim>>
325 create_laplace_tensor_product_matrix(
326 const FiniteElement<1> &fe,
327 const Quadrature<1> &quadrature,
328 const ::ndarray<LaplaceBoundaryType, dim, 2> &boundary_ids,
329 const ::ndarray<double, dim, 3> &cell_extent,
330 const unsigned int n_overlap)
331 {
332 // 1) create element mass and siffness matrix (without overlap)
333 const auto create_reference_mass_and_stiffness_matrices =
334 internal::create_reference_mass_and_stiffness_matrices<Number>(
335 fe, quadrature);
336
337 const auto &M_ref =
338 std::get<0>(create_reference_mass_and_stiffness_matrices);
339 const auto &K_ref =
340 std::get<1>(create_reference_mass_and_stiffness_matrices);
341 const auto &is_dg =
342 std::get<2>(create_reference_mass_and_stiffness_matrices);
343
344 AssertIndexRange(n_overlap, M_ref.n());
345 AssertIndexRange(0, n_overlap);
346 AssertThrow(is_dg == false, ExcNotImplemented());
347
348 // 2) loop over all dimensions and create 1d mass and stiffness
349 // matrices so that boundary conditions and overlap are considered
350
351 const unsigned int n_dofs_1D = M_ref.n();
352 const unsigned int n_dofs_1D_with_overlap = M_ref.n() - 2 + 2 * n_overlap;
353
354 std::array<FullMatrix<Number>, dim> Ms;
355 std::array<FullMatrix<Number>, dim> Ks;
356
357 for (unsigned int d = 0; d < dim; ++d)
358 {
359 Ms[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
360 Ks[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap);
361
362 // inner cell
363 for (unsigned int i = 0; i < n_dofs_1D; ++i)
364 for (unsigned int j = 0; j < n_dofs_1D; ++j)
365 {
366 const unsigned int i0 = i + n_overlap - 1;
367 const unsigned int j0 = j + n_overlap - 1;
368 Ms[d][i0][j0] = M_ref[i][j] * cell_extent[d][1];
369 Ks[d][i0][j0] = K_ref[i][j] / cell_extent[d][1];
370 }
371
372 // left neighbor or left boundary
373 if (boundary_ids[d][0] == LaplaceBoundaryType::internal_boundary)
374 {
375 // left neighbor
376 Assert(cell_extent[d][0] > 0.0, ExcInternalError());
377
378 for (unsigned int i = 0; i < n_overlap; ++i)
379 for (unsigned int j = 0; j < n_overlap; ++j)
380 {
381 const unsigned int i0 = n_dofs_1D - n_overlap + i;
382 const unsigned int j0 = n_dofs_1D - n_overlap + j;
383 Ms[d][i][j] += M_ref[i0][j0] * cell_extent[d][0];
384 Ks[d][i][j] += K_ref[i0][j0] / cell_extent[d][0];
385 }
386 }
387 else
388 {
389 if (boundary_ids[d][0] == LaplaceBoundaryType::dirichlet)
390 {
391 // left DBC
392 const unsigned i0 = n_overlap - 1;
393 internal::clear_row_and_column(n_dofs_1D_with_overlap,
394 i0,
395 Ms[d]);
396 internal::clear_row_and_column(n_dofs_1D_with_overlap,
397 i0,
398 Ks[d]);
399 }
400 else if (boundary_ids[d][0] == LaplaceBoundaryType::neumann)
401 {
402 // left NBC -> nothing to do
403 }
404 else
405 {
406 AssertThrow(false, ExcNotImplemented());
407 }
408 }
409
410 // right neighbor or right boundary
411 if (boundary_ids[d][1] == LaplaceBoundaryType::internal_boundary)
412 {
413 Assert(cell_extent[d][2] > 0.0, ExcInternalError());
414
415 for (unsigned int i = 0; i < n_overlap; ++i)
416 for (unsigned int j = 0; j < n_overlap; ++j)
417 {
418 const unsigned int i0 = n_overlap + n_dofs_1D + i - 2;
419 const unsigned int j0 = n_overlap + n_dofs_1D + j - 2;
420 Ms[d][i0][j0] += M_ref[i][j] * cell_extent[d][2];
421 Ks[d][i0][j0] += K_ref[i][j] / cell_extent[d][2];
422 }
423 }
424 else
425 {
426 if (boundary_ids[d][1] == LaplaceBoundaryType::dirichlet)
427 {
428 // right DBC
429 const unsigned i0 = n_overlap + n_dofs_1D - 2;
430 internal::clear_row_and_column(n_dofs_1D_with_overlap,
431 i0,
432 Ms[d]);
433 internal::clear_row_and_column(n_dofs_1D_with_overlap,
434 i0,
435 Ks[d]);
436 }
437 else if (boundary_ids[d][1] == LaplaceBoundaryType::neumann)
438 {
439 // right NBC -> nothing to do
440 }
441 else
442 {
443 AssertThrow(false, ExcNotImplemented());
444 }
445 }
446 }
447
448 return {Ms, Ks};
449 }
450
451
452 template <int dim, typename Number>
453 std::pair<std::array<FullMatrix<Number>, dim>,
454 std::array<FullMatrix<Number>, dim>>
455 create_laplace_tensor_product_matrix(
456 const typename Triangulation<dim>::cell_iterator &cell,
457 const std::set<types::boundary_id> &dirichlet_boundaries,
458 const std::set<types::boundary_id> &neumann_boundaries,
459 const FiniteElement<1> &fe,
460 const Quadrature<1> &quadrature,
461 const ::ndarray<double, dim, 3> &cell_extent,
462 const unsigned int n_overlap)
463 {
464 ::ndarray<LaplaceBoundaryType, dim, 2> boundary_ids;
465
466 for (unsigned int d = 0; d < dim; ++d)
467 {
468 // left neighbor or left boundary
469 if ((cell->at_boundary(2 * d) == false) ||
470 cell->has_periodic_neighbor(2 * d))
471 {
472 // left neighbor
473 Assert(cell_extent[d][0] > 0.0, ExcInternalError());
474
475 boundary_ids[d][0] = LaplaceBoundaryType::internal_boundary;
476 }
477 else
478 {
479 const auto bid = cell->face(2 * d)->boundary_id();
480 if (dirichlet_boundaries.find(bid) !=
481 dirichlet_boundaries.end() /*DBC*/)
482 {
483 // left DBC
484 boundary_ids[d][0] = LaplaceBoundaryType::dirichlet;
485 }
486 else if (neumann_boundaries.find(bid) !=
487 neumann_boundaries.end() /*NBC*/)
488 {
489 // left NBC
490 boundary_ids[d][0] = LaplaceBoundaryType::neumann;
491 }
492 else
493 {
494 AssertThrow(false, ExcNotImplemented());
495 }
496 }
497
498 // right neighbor or right boundary
499 if ((cell->at_boundary(2 * d + 1) == false) ||
500 cell->has_periodic_neighbor(2 * d + 1))
501 {
502 Assert(cell_extent[d][2] > 0.0, ExcInternalError());
503
504 boundary_ids[d][1] = LaplaceBoundaryType::internal_boundary;
505 }
506 else
507 {
508 const auto bid = cell->face(2 * d + 1)->boundary_id();
509 if (dirichlet_boundaries.find(bid) !=
510 dirichlet_boundaries.end() /*DBC*/)
511 {
512 // right DBC
513 boundary_ids[d][1] = LaplaceBoundaryType::dirichlet;
514 }
515 else if (neumann_boundaries.find(bid) !=
516 neumann_boundaries.end() /*NBC*/)
517 {
518 // right NBC
519 boundary_ids[d][1] = LaplaceBoundaryType::neumann;
520 }
521 else
522 {
523 AssertThrow(false, ExcNotImplemented());
524 }
525 }
526 }
527
528 return create_laplace_tensor_product_matrix<dim, Number>(
529 fe, quadrature, boundary_ids, cell_extent, n_overlap);
530 }
531
532 template <typename Number>
533 FullMatrix<Number>
534 create_1d_cell_mass_matrix(const FiniteElement<1> &fe,
535 const Number &h,
536 const std::pair<bool, bool> include_endpoints,
537 std::vector<unsigned int> numbering)
538 {
539 if (dynamic_cast<const FE_DGQ<1> *>(&fe) == nullptr &&
540 numbering.size() == 0)
541 {
542 Assert(
543 include_endpoints.first == true && include_endpoints.second == true,
544 ExcMessage(
545 "You tried to generate a 1D mass matrix with excluding boundary "
546 "dofs for a non-DGQ element without providing a numbering."));
547 }
548
549 if (numbering.size() == 0)
550 {
551 numbering.resize(fe.dofs_per_cell);
552 std::iota(numbering.begin(), numbering.end(), 0);
553 }
554
555 const unsigned int degree = fe.degree;
556 const unsigned int n_dofs_per_cell = fe.dofs_per_cell;
557 QGauss<1> quadrature(degree + 1);
558
559 FullMatrix<Number> cell_matrix(n_dofs_per_cell, n_dofs_per_cell);
560 cell_matrix = 0;
561
562 unsigned int start_dof = include_endpoints.first ? 0 : 1;
563 unsigned int end_dof =
564 include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1;
565 const unsigned int shift = include_endpoints.first ? 0 : 1;
566
567 for (unsigned int i = start_dof; i < end_dof; ++i)
568 for (unsigned int j = start_dof; j < end_dof; ++j)
569 for (unsigned int q = 0; q < quadrature.size(); ++q)
570 cell_matrix(i - shift, j - shift) +=
571 (fe.shape_value(numbering[i], quadrature.point(q)) *
572 fe.shape_value(numbering[j], quadrature.point(q))) *
573 (h * quadrature.weight(q));
574
575 return cell_matrix;
576 }
577
578 template <typename Number>
579 FullMatrix<Number>
580 create_1d_cell_laplace_matrix(const FiniteElement<1> &fe,
581 const Number &h,
582 const std::pair<bool, bool> include_endpoints,
583 std::vector<unsigned int> numbering)
584 {
585 if (dynamic_cast<const FE_DGQ<1> *>(&fe) == nullptr &&
586 numbering.size() == 0)
587 {
588 Assert(
589 include_endpoints.first == true && include_endpoints.second == true,
590 ExcMessage(
591 "You tried to generate a 1D derivative matrix with excluding boundary "
592 "dofs for a non-DGQ element without providing a numbering."));
593 }
594
595 if (numbering.size() == 0)
596 {
597 numbering.resize(fe.dofs_per_cell);
598 std::iota(numbering.begin(), numbering.end(), 0);
599 }
600
601 const unsigned int degree = fe.degree;
602 const unsigned int n_dofs_per_cell = fe.dofs_per_cell;
603 const Number &JxW = h;
604 QGauss<1> quadrature(degree + 1);
605
606 FullMatrix<Number> cell_matrix(n_dofs_per_cell, n_dofs_per_cell);
607 cell_matrix = 0;
608
609 unsigned int start_dof = include_endpoints.first ? 0 : 1;
610 unsigned int end_dof =
611 include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1;
612 const unsigned int shift = include_endpoints.first ? 0 : 1;
613
614 for (unsigned int i = start_dof; i < end_dof; ++i)
615 for (unsigned int j = start_dof; j < end_dof; ++j)
616 for (unsigned int q = 0; q < quadrature.size(); ++q)
617 cell_matrix(i - shift, j - shift) +=
618 (fe.shape_grad(numbering[i], quadrature.point(q)) / h *
619 fe.shape_grad(numbering[j], quadrature.point(q))) /
620 h * (h * quadrature.weight(q));
621
622 return cell_matrix;
623 }
624
625 template <typename Number>
626 FullMatrix<Number>
627 create_1D_discretization_matrix(FullMatrix<Number> &cell_matrix,
628 const unsigned int &n_cells,
629 const unsigned int &overlap,
630 const std::pair<bool, bool> include_endpoints)
631 {
632 const unsigned int n_dofs_per_cell = cell_matrix.n();
633
634 Assert(cell_matrix.m() == n_dofs_per_cell,
635 ExcMessage(
636 "The provided cell mass matrix must be a square matrix."));
637 AssertThrow(
638 n_cells <= 10,
639 ExcMessage(
640 "create_1D_discretization_matrix() returns a full matrix and is not meant to be used with a larger number of cells. "));
641 Assert(n_cells > 0,
642 ExcMessage("You are trying to get a mass matrix of zero cells."));
643 Assert(overlap < n_dofs_per_cell,
644 ExcMessage("The overlap must be smaller than the number of dofs."));
645
646 unsigned int n_total_dofs =
647 n_cells * n_dofs_per_cell - overlap * (n_cells - 1);
648
649 if (!include_endpoints.first)
650 n_total_dofs -= 1;
651 if (!include_endpoints.second)
652 n_total_dofs -= 1;
653
654 FullMatrix<Number> result_matrix(n_total_dofs, n_total_dofs);
655 result_matrix = 0;
656
657 const unsigned int left_shift = include_endpoints.first ? 0 : 1;
658
659 for (unsigned int cell = 0; cell < n_cells; ++cell)
660 {
661 const unsigned int dof_shift = cell * overlap + left_shift;
662
663 const unsigned int start_dof =
664 (cell == 0 && !include_endpoints.first) ? 1 : 0;
665
666 const unsigned int end_dof =
667 (cell == n_cells - 1 && !include_endpoints.second) ?
668 n_dofs_per_cell - 1 :
669 n_dofs_per_cell;
670 for (unsigned int i = start_dof; i < end_dof; ++i)
671 for (unsigned int j = start_dof; j < end_dof; ++j)
672 {
673 result_matrix(i + cell * n_dofs_per_cell - dof_shift,
674 j + cell * n_dofs_per_cell - dof_shift) +=
675 cell_matrix(i, j);
676 }
677 }
678 return result_matrix;
679 }
680
681
682
683 template <typename Number>
684 FullMatrix<Number>
685 create_1d_ghost_penalty_matrix(const FiniteElement<1> &fe,
686 const Number h,
687 std::vector<Number> coefficients)
688 {
689 Assert(dynamic_cast<const FE_Q<1> *>(&fe) != nullptr, ExcNotImplemented());
690 Assert(h > 0, ExcMessage("Provided element size h is negative"));
691
692 const unsigned int degree = fe.degree;
693 Assert(degree > 0,
694 ExcMessage("Provided element degree has to greater than 0"));
695
696
697 Assert(coefficients.size() == 0 || coefficients.size() == degree,
698 ExcMessage(
699 "Provided coefficients vector has to be empty or the same size "
700 "as the number of dofs"));
701
702 if (coefficients.size() == 0)
703 {
704 coefficients.resize(degree);
705
706 double inverse_factorial_square = 1.;
707 coefficients[0] = 1.;
708 for (unsigned int k = 2; k <= degree; ++k)
709 {
710 inverse_factorial_square /= (k * k);
711 coefficients[k - 1] = inverse_factorial_square;
712 }
713 }
714
715 std::vector<std::vector<Polynomials::Polynomial<double>>> polynomial_basis;
716
717 polynomial_basis.resize(degree + 1);
718
719 auto support_points = fe.get_unit_support_points();
720 std::sort(support_points.begin(),
721 support_points.end(),
722 [](const Point<1> &p, const Point<1> &q) -> bool {
723 return p(0) < q(0);
724 });
725
726 polynomial_basis[0] =
727 Polynomials::generate_complete_Lagrange_basis(support_points);
728
729 for (unsigned int k = 1; k < degree + 1; ++k)
730 {
731 polynomial_basis[k].reserve(degree + 1);
732 for (unsigned int i = 0; i < degree + 1; ++i)
733 polynomial_basis[k].push_back(
734 polynomial_basis[k - 1][i].derivative());
735 }
736
737
738 FullMatrix<Number> penalty_matrix =
739 create_1d_ghost_penalty_matrix(polynomial_basis[1]);
740 penalty_matrix *= coefficients[0];
741
742 for (unsigned int k = 2; k < degree + 1; ++k)
743 {
744 FullMatrix<Number> kth_matrix =
745 create_1d_ghost_penalty_matrix(polynomial_basis[k]);
746 penalty_matrix.add(coefficients[k - 1], kth_matrix);
747 }
748
749 penalty_matrix *= (1 / h);
750 return penalty_matrix;
751 }
752
753
754 template <typename Number>
755 FullMatrix<Number>
756 create_1d_ghost_penalty_matrix(
757 const std::vector<Polynomials::Polynomial<double>>
758 &polynomial_basis_derivative,
759 const unsigned int overlap)
760 {
761 const unsigned int n_dofs_per_cell = polynomial_basis_derivative.size();
762 const unsigned int n_total_dofs = 2 * n_dofs_per_cell - overlap;
763 const unsigned int shift = n_dofs_per_cell - overlap;
764
765 FullMatrix<Number> penalty_matrix(n_total_dofs, n_total_dofs);
766
767 std::vector<double> values_left(n_dofs_per_cell);
768 std::vector<double> values_right(n_dofs_per_cell);
769
770 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
771 {
772 values_left[i] = polynomial_basis_derivative[i].value(0);
773 values_right[i] = polynomial_basis_derivative[i].value(1);
774 }
775
776 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
777 for (unsigned int j = 0; j < n_dofs_per_cell; ++j)
778 penalty_matrix(i, j) += values_right[i] * values_right[j];
779
780
781 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
782 for (unsigned int j = 0; j < n_dofs_per_cell; ++j)
783 penalty_matrix(i + shift, j) -= values_left[i] * values_right[j];
784
785 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
786 for (unsigned int j = 0; j < n_dofs_per_cell; ++j)
787 penalty_matrix(i, j + shift) -= values_right[i] * values_left[j];
788
789 for (unsigned int i = 0; i < n_dofs_per_cell; ++i)
790 for (unsigned int j = 0; j < n_dofs_per_cell; ++j)
791 penalty_matrix(i + shift, j + shift) += values_left[i] * values_left[j];
792
793 return penalty_matrix;
794 }
795} // namespace TensorProductMatrixCreator
796
797
798
799DEAL_II_NAMESPACE_CLOSE
800
801#endif
push_back
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
DoFHandler::distribute_dofs
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
FEValuesBase::dof_indices
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEValuesBase::quadrature_point_indices
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
FEValuesBase::shape_grad
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
FEValuesBase::JxW
double JxW(const unsigned int q_point) const
FEValuesBase::shape_value
const double & shape_value(const unsigned int i, const unsigned int q_point) const
FEValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
FE_DGQ
Definition fe_dgq.h:110
FE_Q
Definition fe_q.h:552
FiniteElementData::degree
const unsigned int degree
Definition fe_data.h:450
FiniteElementData::n_dofs_per_cell
unsigned int n_dofs_per_cell() const
FiniteElementData::tensor_degree
unsigned int tensor_degree() const
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition fe_data.h:434
FiniteElement::shape_grad
virtual Tensor< 1, dim > shape_grad(const unsigned int i, const Point< dim > &p) const
FiniteElement::get_unit_support_points
const std::vector< Point< dim > > & get_unit_support_points() const
FiniteElement::shape_value
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
FullMatrix::add
void add(const number a, const FullMatrix< number2 > &A)
Point
Definition point.h:111
Triangulation::begin
cell_iterator begin(const unsigned int level=0) const
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
dof_handler.h
fe_values.h
fe_dgq.h
fe_q.h
fe_tools.h
grid_generator.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:73
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:132
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:79
mapping_q1.h
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
Polynomials::generate_complete_Lagrange_basis
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
Definition polynomial.cc:692
TensorProductMatrixCreator::internal::clear_row_and_column
void clear_row_and_column(const unsigned int n_dofs_1D_with_overlap, const unsigned int n, FullMatrix< Number > &matrix)
Definition tensor_product_matrix_creator.h:253
TensorProductMatrixCreator::internal::create_reference_mass_and_stiffness_matrices
std::tuple< FullMatrix< Number >, FullMatrix< Number >, bool > create_reference_mass_and_stiffness_matrices(const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature)
Definition tensor_product_matrix_creator.h:266
TensorProductMatrixCreator::create_laplace_tensor_product_matrix
std::pair< std::array< FullMatrix< Number >, dim >, std::array< FullMatrix< Number >, dim > > create_laplace_tensor_product_matrix(const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature, const ::ndarray< LaplaceBoundaryType, dim, 2 > &boundary_ids, const ::ndarray< double, dim, 3 > &cell_extent, const unsigned int n_overlap=1)
Definition tensor_product_matrix_creator.h:325
TensorProductMatrixCreator::create_1d_cell_mass_matrix
FullMatrix< Number > create_1d_cell_mass_matrix(const FiniteElement< 1 > &fe, const Number &h, const std::pair< bool, bool > include_endpoints={true, true}, std::vector< unsigned int > numbering=std::vector< unsigned int >())
Definition tensor_product_matrix_creator.h:534
TensorProductMatrixCreator::create_1D_discretization_matrix
FullMatrix< Number > create_1D_discretization_matrix(FullMatrix< Number > &cell_matrix, const unsigned int &n_cells, const unsigned int &overlap, const std::pair< bool, bool > include_endpoints={true, true})
Definition tensor_product_matrix_creator.h:627
TensorProductMatrixCreator::create_1d_cell_laplace_matrix
FullMatrix< Number > create_1d_cell_laplace_matrix(const FiniteElement< 1 > &fe, const Number &h, const std::pair< bool, bool > include_endpoints={true, true}, std::vector< unsigned int > numbering=std::vector< unsigned int >())
Definition tensor_product_matrix_creator.h:580
TensorProductMatrixCreator::create_1d_ghost_penalty_matrix
FullMatrix< Number > create_1d_ghost_penalty_matrix(const FiniteElement< 1 > &fe, const Number h, std::vector< Number > coefficients=std::vector< Number >())
Create a 1D ghost penalty matrix for a given finite element. Ghost penalty is used for stabilization ...
Definition tensor_product_matrix_creator.h:685
Utilities::invert_permutation
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1670
ndarray
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:105
quadrature_lib.h
Generated by doxygen 1.9.8