Reference documentation for deal.II version 9.6.0
\(\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
assembler.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2010 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_mesh_worker_assembler_h
17#define dealii_mesh_worker_assembler_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/algorithms/any_data.h>
22
23#include <deal.II/base/mg_level_object.h>
24#include <deal.II/base/smartpointer.h>
25
26#include <deal.II/lac/block_vector.h>
27
28#include <deal.II/meshworker/dof_info.h>
29#include <deal.II/meshworker/functional.h>
30#include <deal.II/meshworker/simple.h>
31
32#include <deal.II/multigrid/mg_constrained_dofs.h>
33
34
35DEAL_II_NAMESPACE_OPEN
36
37namespace MeshWorker
38{
87 namespace Assembler
88 {
109 template <typename VectorType>
110 class ResidualLocalBlocksToGlobalBlocks
111 {
112 public:
116 void
117 initialize(const BlockInfo *block_info, AnyData &residuals);
118
122 void
123 initialize(
124 const AffineConstraints<typename VectorType::value_type> &constraints);
125
133 template <class DOFINFO>
134 void
135 initialize_info(DOFINFO &info, bool face) const;
136
137
141 template <class DOFINFO>
142 void
143 assemble(const DOFINFO &info);
144
148 template <class DOFINFO>
149 void
150 assemble(const DOFINFO &info1, const DOFINFO &info2);
151
152 private:
156 void
157 assemble(VectorType &global,
158 const BlockVector<double> &local,
159 const std::vector<types::global_dof_index> &dof);
160
164 AnyData residuals;
165
169 SmartPointer<const BlockInfo,
170 ResidualLocalBlocksToGlobalBlocks<VectorType>>
171 block_info;
172
176 SmartPointer<const AffineConstraints<typename VectorType::value_type>,
177 ResidualLocalBlocksToGlobalBlocks<VectorType>>
178 constraints;
179 };
180
181
208 template <typename MatrixType, typename number = double>
209 class MatrixLocalBlocksToGlobalBlocks
210 {
211 public:
216 MatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
217
222 void
223 initialize(const BlockInfo *block_info,
224 MatrixBlockVector<MatrixType> &matrices);
225
229 void
230 initialize(
231 const AffineConstraints<typename MatrixType::value_type> &constraints);
232
240 template <class DOFINFO>
241 void
242 initialize_info(DOFINFO &info, bool face) const;
243
244
248 template <class DOFINFO>
249 void
250 assemble(const DOFINFO &info);
251
255 template <class DOFINFO>
256 void
257 assemble(const DOFINFO &info1, const DOFINFO &info2);
258
259 private:
263 void
264 assemble(MatrixBlock<MatrixType> &global,
265 const FullMatrix<number> &local,
266 const unsigned int block_row,
267 const unsigned int block_col,
268 const std::vector<types::global_dof_index> &dof1,
269 const std::vector<types::global_dof_index> &dof2);
270
274 SmartPointer<MatrixBlockVector<MatrixType>,
275 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>>
276 matrices;
277
281 SmartPointer<const BlockInfo,
282 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>>
283 block_info;
288 SmartPointer<const AffineConstraints<typename MatrixType::value_type>,
289 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>>
290 constraints;
291
297 const double threshold;
298 };
299
323 template <typename MatrixType, typename number = double>
324 class MGMatrixLocalBlocksToGlobalBlocks
325 {
326 public:
327 using MatrixPtrVector = MGMatrixBlockVector<MatrixType>;
328 using MatrixPtrVectorPtr =
329 SmartPointer<MatrixPtrVector,
330 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>>;
331
336 MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
337
342 void
343 initialize(const BlockInfo *block_info, MatrixPtrVector &matrices);
344
348 void
349 initialize(const MGConstrainedDoFs &mg_constrained_dofs);
350
356 void
357 initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down);
358
364 void
365 initialize_interfaces(MatrixPtrVector &interface_in,
366 MatrixPtrVector &interface_out);
374 template <class DOFINFO>
375 void
376 initialize_info(DOFINFO &info, bool face) const;
377
378
382 template <class DOFINFO>
383 void
384 assemble(const DOFINFO &info);
385
389 template <class DOFINFO>
390 void
391 assemble(const DOFINFO &info1, const DOFINFO &info2);
392
393 private:
397 void
398 assemble(MatrixType &global,
399 const FullMatrix<number> &local,
400 const unsigned int block_row,
401 const unsigned int block_col,
402 const std::vector<types::global_dof_index> &dof1,
403 const std::vector<types::global_dof_index> &dof2,
404 const unsigned int level1,
405 const unsigned int level2,
406 bool transpose = false);
407
411 void
412 assemble_fluxes(MatrixType &global,
413 const FullMatrix<number> &local,
414 const unsigned int block_row,
415 const unsigned int block_col,
416 const std::vector<types::global_dof_index> &dof1,
417 const std::vector<types::global_dof_index> &dof2,
418 const unsigned int level1,
419 const unsigned int level2);
420
424 void
425 assemble_up(MatrixType &global,
426 const FullMatrix<number> &local,
427 const unsigned int block_row,
428 const unsigned int block_col,
429 const std::vector<types::global_dof_index> &dof1,
430 const std::vector<types::global_dof_index> &dof2,
431 const unsigned int level1,
432 const unsigned int level2);
433
437 void
438 assemble_down(MatrixType &global,
439 const FullMatrix<number> &local,
440 const unsigned int block_row,
441 const unsigned int block_col,
442 const std::vector<types::global_dof_index> &dof1,
443 const std::vector<types::global_dof_index> &dof2,
444 const unsigned int level1,
445 const unsigned int level2);
446
450 void
451 assemble_in(MatrixType &global,
452 const FullMatrix<number> &local,
453 const unsigned int block_row,
454 const unsigned int block_col,
455 const std::vector<types::global_dof_index> &dof1,
456 const std::vector<types::global_dof_index> &dof2,
457 const unsigned int level1,
458 const unsigned int level2);
459
463 void
464 assemble_out(MatrixType &global,
465 const FullMatrix<number> &local,
466 const unsigned int block_row,
467 const unsigned int block_col,
468 const std::vector<types::global_dof_index> &dof1,
469 const std::vector<types::global_dof_index> &dof2,
470 const unsigned int level1,
471 const unsigned int level2);
472
476 MatrixPtrVectorPtr matrices;
477
482 MatrixPtrVectorPtr flux_down;
483
488 MatrixPtrVectorPtr flux_up;
489
494 MatrixPtrVectorPtr interface_out;
495
500 MatrixPtrVectorPtr interface_in;
501
505 SmartPointer<const BlockInfo,
506 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>>
507 block_info;
508
512 SmartPointer<const MGConstrainedDoFs,
513 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>>
514 mg_constrained_dofs;
515
516
522 const double threshold;
523 };
524
525 //----------------------------------------------------------------------//
526
527 template <typename VectorType>
528 inline void
529 ResidualLocalBlocksToGlobalBlocks<VectorType>::initialize(
530 const BlockInfo *b,
531 AnyData &m)
532 {
533 block_info = b;
534 residuals = m;
535 }
536
537 template <typename VectorType>
538 inline void
544
545
546 template <typename VectorType>
547 template <class DOFINFO>
548 inline void
549 ResidualLocalBlocksToGlobalBlocks<VectorType>::initialize_info(
550 DOFINFO &info,
551 bool) const
552 {
553 info.initialize_vectors(residuals.size());
554 }
555
556 template <typename VectorType>
557 inline void
558 ResidualLocalBlocksToGlobalBlocks<VectorType>::assemble(
559 VectorType &global,
560 const BlockVector<double> &local,
561 const std::vector<types::global_dof_index> &dof)
562 {
563 if (constraints == 0)
564 {
565 for (unsigned int b = 0; b < local.n_blocks(); ++b)
566 for (unsigned int j = 0; j < local.block(b).size(); ++j)
567 {
568 // The coordinates of
569 // the current entry in
570 // DoFHandler
571 // numbering, which
572 // differs from the
573 // block-wise local
574 // numbering we use in
575 // our local vectors
576 const unsigned int jcell =
577 this->block_info->local().local_to_global(b, j);
578 global(dof[jcell]) += local.block(b)(j);
579 }
580 }
581 else
582 constraints->distribute_local_to_global(local, dof, global);
583 }
584
585
586 template <typename VectorType>
587 template <class DOFINFO>
588 inline void
589 ResidualLocalBlocksToGlobalBlocks<VectorType>::assemble(const DOFINFO &info)
590 {
591 for (unsigned int i = 0; i < residuals.size(); ++i)
592 assemble(*(residuals.entry<VectorType>(i)),
593 info.vector(i),
594 info.indices);
595 }
596
597
598 template <typename VectorType>
599 template <class DOFINFO>
600 inline void
601 ResidualLocalBlocksToGlobalBlocks<VectorType>::assemble(
602 const DOFINFO &info1,
603 const DOFINFO &info2)
604 {
605 for (unsigned int i = 0; i < residuals.size(); ++i)
606 {
607 assemble(*(residuals.entry<VectorType>(i)),
608 info1.vector(i),
609 info1.indices);
610 assemble(*(residuals.entry<VectorType>(i)),
611 info2.vector(i),
612 info2.indices);
613 }
614 }
615
616
617 //----------------------------------------------------------------------//
618
619 template <typename MatrixType, typename number>
620 inline MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::
621 MatrixLocalBlocksToGlobalBlocks(double threshold)
622 : threshold(threshold)
623 {}
624
625
626 template <typename MatrixType, typename number>
627 inline void
628 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize(
629 const BlockInfo *b,
630 MatrixBlockVector<MatrixType> &m)
631 {
632 block_info = b;
633 matrices = &m;
634 }
635
636
637
638 template <typename MatrixType, typename number>
639 inline void
645
646
647
648 template <typename MatrixType, typename number>
649 template <class DOFINFO>
650 inline void
651 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize_info(
652 DOFINFO &info,
653 bool face) const
654 {
655 info.initialize_matrices(*matrices, face);
656 }
657
658
659
660 template <typename MatrixType, typename number>
661 inline void
662 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(
663 MatrixBlock<MatrixType> &global,
664 const FullMatrix<number> &local,
665 const unsigned int block_row,
666 const unsigned int block_col,
667 const std::vector<types::global_dof_index> &dof1,
668 const std::vector<types::global_dof_index> &dof2)
669 {
670 if (constraints == nullptr)
671 {
672 for (unsigned int j = 0; j < local.n_rows(); ++j)
673 for (unsigned int k = 0; k < local.n_cols(); ++k)
674 if (std::fabs(local(j, k)) >= threshold)
675 {
676 // The coordinates of
677 // the current entry in
678 // DoFHandler
679 // numbering, which
680 // differs from the
681 // block-wise local
682 // numbering we use in
683 // our local matrices
684 const unsigned int jcell =
685 this->block_info->local().local_to_global(block_row, j);
686 const unsigned int kcell =
687 this->block_info->local().local_to_global(block_col, k);
688
689 global.add(dof1[jcell], dof2[kcell], local(j, k));
690 }
691 }
692 else
693 {
694 const BlockIndices &bi = this->block_info->local();
695 std::vector<types::global_dof_index> sliced_row_indices(
696 bi.block_size(block_row));
697 for (unsigned int i = 0; i < sliced_row_indices.size(); ++i)
698 sliced_row_indices[i] = dof1[bi.block_start(block_row) + i];
699
700 std::vector<types::global_dof_index> sliced_col_indices(
701 bi.block_size(block_col));
702 for (unsigned int i = 0; i < sliced_col_indices.size(); ++i)
703 sliced_col_indices[i] = dof2[bi.block_start(block_col) + i];
704
705 constraints->distribute_local_to_global(local,
706 sliced_row_indices,
707 sliced_col_indices,
708 global);
709 }
710 }
711
712
713 template <typename MatrixType, typename number>
714 template <class DOFINFO>
715 inline void
716 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(
717 const DOFINFO &info)
718 {
719 for (unsigned int i = 0; i < matrices->size(); ++i)
720 {
721 // Row and column index of
722 // the block we are dealing with
723 const types::global_dof_index row = matrices->block(i).row;
724 const types::global_dof_index col = matrices->block(i).column;
725
726 assemble(matrices->block(i),
727 info.matrix(i, false).matrix,
728 row,
729 col,
730 info.indices,
731 info.indices);
732 }
733 }
734
735
736 template <typename MatrixType, typename number>
737 template <class DOFINFO>
738 inline void
739 MatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(
740 const DOFINFO &info1,
741 const DOFINFO &info2)
742 {
743 for (unsigned int i = 0; i < matrices->size(); ++i)
744 {
745 // Row and column index of
746 // the block we are dealing with
747 const types::global_dof_index row = matrices->block(i).row;
748 const types::global_dof_index col = matrices->block(i).column;
749
750 assemble(matrices->block(i),
751 info1.matrix(i, false).matrix,
752 row,
753 col,
754 info1.indices,
755 info1.indices);
756 assemble(matrices->block(i),
757 info1.matrix(i, true).matrix,
758 row,
759 col,
760 info1.indices,
761 info2.indices);
762 assemble(matrices->block(i),
763 info2.matrix(i, false).matrix,
764 row,
765 col,
766 info2.indices,
767 info2.indices);
768 assemble(matrices->block(i),
769 info2.matrix(i, true).matrix,
770 row,
771 col,
772 info2.indices,
773 info1.indices);
774 }
775 }
776
777
778 // ----------------------------------------------------------------------//
779
780 template <typename MatrixType, typename number>
781 inline MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::
782 MGMatrixLocalBlocksToGlobalBlocks(double threshold)
783 : threshold(threshold)
784 {}
785
786
787 template <typename MatrixType, typename number>
788 inline void
789 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize(
790 const BlockInfo *b,
791 MatrixPtrVector &m)
792 {
793 block_info = b;
794 AssertDimension(block_info->local().size(), block_info->global().size());
795 matrices = &m;
796 }
797
798
799 template <typename MatrixType, typename number>
800 inline void
801 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize(
802 const MGConstrainedDoFs &mg_c)
803 {
804 mg_constrained_dofs = &mg_c;
805 }
806
807
808 template <typename MatrixType, typename number>
809 template <class DOFINFO>
810 inline void
811 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize_info(
812 DOFINFO &info,
813 bool face) const
814 {
815 info.initialize_matrices(*matrices, face);
816 }
817
818
819
820 template <typename MatrixType, typename number>
821 inline void
822 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::initialize_edge_flux(
823 MatrixPtrVector &up,
824 MatrixPtrVector &down)
825 {
826 flux_up = up;
827 flux_down = down;
828 }
829
830
831 template <typename MatrixType, typename number>
832 inline void
833 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::
834 initialize_interfaces(MatrixPtrVector &in, MatrixPtrVector &out)
835 {
836 interface_in = in;
837 interface_out = out;
838 }
839
840
841 template <typename MatrixType, typename number>
842 inline void
843 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(
844 MatrixType &global,
845 const FullMatrix<number> &local,
846 const unsigned int block_row,
847 const unsigned int block_col,
848 const std::vector<types::global_dof_index> &dof1,
849 const std::vector<types::global_dof_index> &dof2,
850 const unsigned int level1,
851 const unsigned int level2,
852 bool transpose)
853 {
854 for (unsigned int j = 0; j < local.n_rows(); ++j)
855 for (unsigned int k = 0; k < local.n_cols(); ++k)
856 if (std::fabs(local(j, k)) >= threshold)
857 {
858 // The coordinates of
859 // the current entry in
860 // DoFHandler
861 // numbering, which
862 // differs from the
863 // block-wise local
864 // numbering we use in
865 // our local matrices
866 const unsigned int jcell =
867 this->block_info->local().local_to_global(block_row, j);
868 const unsigned int kcell =
869 this->block_info->local().local_to_global(block_col, k);
870
871 // The global dof
872 // indices to assemble
873 // in. Since we may
874 // have face matrices
875 // coupling two
876 // different cells, we
877 // provide two sets of
878 // dof indices.
879 const unsigned int jglobal = this->block_info->level(level1)
880 .global_to_local(dof1[jcell])
881 .second;
882 const unsigned int kglobal = this->block_info->level(level2)
883 .global_to_local(dof2[kcell])
884 .second;
885
886 if (mg_constrained_dofs == 0)
887 {
888 if (transpose)
889 global.add(kglobal, jglobal, local(j, k));
890 else
891 global.add(jglobal, kglobal, local(j, k));
892 }
893 else
894 {
895 if (!mg_constrained_dofs->at_refinement_edge(level1,
896 jglobal) &&
897 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
898 {
899 if (mg_constrained_dofs->set_boundary_values())
900 {
901 if ((!mg_constrained_dofs->is_boundary_index(
902 level1, jglobal) &&
903 !mg_constrained_dofs->is_boundary_index(
904 level2, kglobal)) ||
905 (mg_constrained_dofs->is_boundary_index(
906 level1, jglobal) &&
907 mg_constrained_dofs->is_boundary_index(
908 level2, kglobal) &&
909 jglobal == kglobal))
910 {
911 if (transpose)
912 global.add(kglobal, jglobal, local(j, k));
913 else
914 global.add(jglobal, kglobal, local(j, k));
915 }
916 }
917 else
918 {
919 if (transpose)
920 global.add(kglobal, jglobal, local(j, k));
921 else
922 global.add(jglobal, kglobal, local(j, k));
923 }
924 }
925 }
926 }
927 }
928
929
930 template <typename MatrixType, typename number>
931 inline void
932 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_fluxes(
933 MatrixType &global,
934 const FullMatrix<number> &local,
935 const unsigned int block_row,
936 const unsigned int block_col,
937 const std::vector<types::global_dof_index> &dof1,
938 const std::vector<types::global_dof_index> &dof2,
939 const unsigned int level1,
940 const unsigned int level2)
941 {
942 for (unsigned int j = 0; j < local.n_rows(); ++j)
943 for (unsigned int k = 0; k < local.n_cols(); ++k)
944 if (std::fabs(local(j, k)) >= threshold)
945 {
946 // The coordinates of
947 // the current entry in
948 // DoFHandler
949 // numbering, which
950 // differs from the
951 // block-wise local
952 // numbering we use in
953 // our local matrices
954 const unsigned int jcell =
955 this->block_info->local().local_to_global(block_row, j);
956 const unsigned int kcell =
957 this->block_info->local().local_to_global(block_col, k);
958
959 // The global dof
960 // indices to assemble
961 // in. Since we may
962 // have face matrices
963 // coupling two
964 // different cells, we
965 // provide two sets of
966 // dof indices.
967 const unsigned int jglobal = this->block_info->level(level1)
968 .global_to_local(dof1[jcell])
969 .second;
970 const unsigned int kglobal = this->block_info->level(level2)
971 .global_to_local(dof2[kcell])
972 .second;
973
974 if (mg_constrained_dofs == 0)
975 global.add(jglobal, kglobal, local(j, k));
976 else
977 {
978 if (!mg_constrained_dofs->non_refinement_edge_index(
979 level1, jglobal) &&
980 !mg_constrained_dofs->non_refinement_edge_index(level2,
981 kglobal))
982 {
983 if (!mg_constrained_dofs->at_refinement_edge(level1,
984 jglobal) &&
985 !mg_constrained_dofs->at_refinement_edge(level2,
986 kglobal))
987 global.add(jglobal, kglobal, local(j, k));
988 }
989 }
990 }
991 }
992
993 template <typename MatrixType, typename number>
994 inline void
995 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_up(
996 MatrixType &global,
997 const FullMatrix<number> &local,
998 const unsigned int block_row,
999 const unsigned int block_col,
1000 const std::vector<types::global_dof_index> &dof1,
1001 const std::vector<types::global_dof_index> &dof2,
1002 const unsigned int level1,
1003 const unsigned int level2)
1004 {
1005 for (unsigned int j = 0; j < local.n_rows(); ++j)
1006 for (unsigned int k = 0; k < local.n_cols(); ++k)
1007 if (std::fabs(local(j, k)) >= threshold)
1008 {
1009 // The coordinates of
1010 // the current entry in
1011 // DoFHandler
1012 // numbering, which
1013 // differs from the
1014 // block-wise local
1015 // numbering we use in
1016 // our local matrices
1017 const unsigned int jcell =
1018 this->block_info->local().local_to_global(block_row, j);
1019 const unsigned int kcell =
1020 this->block_info->local().local_to_global(block_col, k);
1021
1022 // The global dof
1023 // indices to assemble
1024 // in. Since we may
1025 // have face matrices
1026 // coupling two
1027 // different cells, we
1028 // provide two sets of
1029 // dof indices.
1030 const unsigned int jglobal = this->block_info->level(level1)
1031 .global_to_local(dof1[jcell])
1032 .second;
1033 const unsigned int kglobal = this->block_info->level(level2)
1034 .global_to_local(dof2[kcell])
1035 .second;
1036
1037 if (mg_constrained_dofs == 0)
1038 global.add(jglobal, kglobal, local(j, k));
1039 else
1040 {
1041 if (!mg_constrained_dofs->non_refinement_edge_index(
1042 level1, jglobal) &&
1043 !mg_constrained_dofs->non_refinement_edge_index(level2,
1044 kglobal))
1045 {
1046 if (!mg_constrained_dofs->at_refinement_edge(level1,
1047 jglobal) &&
1048 !mg_constrained_dofs->at_refinement_edge(level2,
1049 kglobal))
1050 global.add(jglobal, kglobal, local(j, k));
1051 }
1052 }
1053 }
1054 }
1055
1056 template <typename MatrixType, typename number>
1057 inline void
1058 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_down(
1059 MatrixType &global,
1060 const FullMatrix<number> &local,
1061 const unsigned int block_row,
1062 const unsigned int block_col,
1063 const std::vector<types::global_dof_index> &dof1,
1064 const std::vector<types::global_dof_index> &dof2,
1065 const unsigned int level1,
1066 const unsigned int level2)
1067 {
1068 for (unsigned int j = 0; j < local.n_rows(); ++j)
1069 for (unsigned int k = 0; k < local.n_cols(); ++k)
1070 if (std::fabs(local(k, j)) >= threshold)
1071 {
1072 // The coordinates of
1073 // the current entry in
1074 // DoFHandler
1075 // numbering, which
1076 // differs from the
1077 // block-wise local
1078 // numbering we use in
1079 // our local matrices
1080 const unsigned int jcell =
1081 this->block_info->local().local_to_global(block_row, j);
1082 const unsigned int kcell =
1083 this->block_info->local().local_to_global(block_col, k);
1084
1085 // The global dof
1086 // indices to assemble
1087 // in. Since we may
1088 // have face matrices
1089 // coupling two
1090 // different cells, we
1091 // provide two sets of
1092 // dof indices.
1093 const unsigned int jglobal = this->block_info->level(level1)
1094 .global_to_local(dof1[jcell])
1095 .second;
1096 const unsigned int kglobal = this->block_info->level(level2)
1097 .global_to_local(dof2[kcell])
1098 .second;
1099
1100 if (mg_constrained_dofs == 0)
1101 global.add(jglobal, kglobal, local(k, j));
1102 else
1103 {
1104 if (!mg_constrained_dofs->non_refinement_edge_index(
1105 level1, jglobal) &&
1106 !mg_constrained_dofs->non_refinement_edge_index(level2,
1107 kglobal))
1108 {
1109 if (!mg_constrained_dofs->at_refinement_edge(level1,
1110 jglobal) &&
1111 !mg_constrained_dofs->at_refinement_edge(level2,
1112 kglobal))
1113 global.add(jglobal, kglobal, local(k, j));
1114 }
1115 }
1116 }
1117 }
1118
1119 template <typename MatrixType, typename number>
1120 inline void
1121 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_in(
1122 MatrixType &global,
1123 const FullMatrix<number> &local,
1124 const unsigned int block_row,
1125 const unsigned int block_col,
1126 const std::vector<types::global_dof_index> &dof1,
1127 const std::vector<types::global_dof_index> &dof2,
1128 const unsigned int level1,
1129 const unsigned int level2)
1130 {
1131 // AssertDimension(local.n(), dof1.size());
1132 // AssertDimension(local.m(), dof2.size());
1133
1134 for (unsigned int j = 0; j < local.n_rows(); ++j)
1135 for (unsigned int k = 0; k < local.n_cols(); ++k)
1136 if (std::fabs(local(j, k)) >= threshold)
1137 {
1138 // The coordinates of
1139 // the current entry in
1140 // DoFHandler
1141 // numbering, which
1142 // differs from the
1143 // block-wise local
1144 // numbering we use in
1145 // our local matrices
1146 const unsigned int jcell =
1147 this->block_info->local().local_to_global(block_row, j);
1148 const unsigned int kcell =
1149 this->block_info->local().local_to_global(block_col, k);
1150
1151 // The global dof
1152 // indices to assemble
1153 // in. Since we may
1154 // have face matrices
1155 // coupling two
1156 // different cells, we
1157 // provide two sets of
1158 // dof indices.
1159 const unsigned int jglobal = this->block_info->level(level1)
1160 .global_to_local(dof1[jcell])
1161 .second;
1162 const unsigned int kglobal = this->block_info->level(level2)
1163 .global_to_local(dof2[kcell])
1164 .second;
1165
1166 if (mg_constrained_dofs == 0)
1167 global.add(jglobal, kglobal, local(j, k));
1168 else
1169 {
1170 if (mg_constrained_dofs->at_refinement_edge(level1,
1171 jglobal) &&
1172 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1173 {
1174 if (mg_constrained_dofs->set_boundary_values())
1175 {
1176 if ((!mg_constrained_dofs->is_boundary_index(
1177 level1, jglobal) &&
1178 !mg_constrained_dofs->is_boundary_index(
1179 level2, kglobal)) ||
1180 (mg_constrained_dofs->is_boundary_index(
1181 level1, jglobal) &&
1182 mg_constrained_dofs->is_boundary_index(
1183 level2, kglobal) &&
1184 jglobal == kglobal))
1185 global.add(jglobal, kglobal, local(j, k));
1186 }
1187 else
1188 global.add(jglobal, kglobal, local(j, k));
1189 }
1190 }
1191 }
1192 }
1193
1194 template <typename MatrixType, typename number>
1195 inline void
1196 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble_out(
1197 MatrixType &global,
1198 const FullMatrix<number> &local,
1199 const unsigned int block_row,
1200 const unsigned int block_col,
1201 const std::vector<types::global_dof_index> &dof1,
1202 const std::vector<types::global_dof_index> &dof2,
1203 const unsigned int level1,
1204 const unsigned int level2)
1205 {
1206 // AssertDimension(local.n(), dof1.size());
1207 // AssertDimension(local.m(), dof2.size());
1208
1209 for (unsigned int j = 0; j < local.n_rows(); ++j)
1210 for (unsigned int k = 0; k < local.n_cols(); ++k)
1211 if (std::fabs(local(k, j)) >= threshold)
1212 {
1213 // The coordinates of
1214 // the current entry in
1215 // DoFHandler
1216 // numbering, which
1217 // differs from the
1218 // block-wise local
1219 // numbering we use in
1220 // our local matrices
1221 const unsigned int jcell =
1222 this->block_info->local().local_to_global(block_row, j);
1223 const unsigned int kcell =
1224 this->block_info->local().local_to_global(block_col, k);
1225
1226 // The global dof
1227 // indices to assemble
1228 // in. Since we may
1229 // have face matrices
1230 // coupling two
1231 // different cells, we
1232 // provide two sets of
1233 // dof indices.
1234 const unsigned int jglobal = this->block_info->level(level1)
1235 .global_to_local(dof1[jcell])
1236 .second;
1237 const unsigned int kglobal = this->block_info->level(level2)
1238 .global_to_local(dof2[kcell])
1239 .second;
1240
1241 if (mg_constrained_dofs == 0)
1242 global.add(jglobal, kglobal, local(k, j));
1243 else
1244 {
1245 if (mg_constrained_dofs->at_refinement_edge(level1,
1246 jglobal) &&
1247 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1248 {
1249 if (mg_constrained_dofs->set_boundary_values())
1250 {
1251 if ((!mg_constrained_dofs->is_boundary_index(
1252 level1, jglobal) &&
1253 !mg_constrained_dofs->is_boundary_index(
1254 level2, kglobal)) ||
1255 (mg_constrained_dofs->is_boundary_index(
1256 level1, jglobal) &&
1257 mg_constrained_dofs->is_boundary_index(
1258 level2, kglobal) &&
1259 jglobal == kglobal))
1260 global.add(jglobal, kglobal, local(k, j));
1261 }
1262 else
1263 global.add(jglobal, kglobal, local(k, j));
1264 }
1265 }
1266 }
1267 }
1268
1269
1270 template <typename MatrixType, typename number>
1271 template <class DOFINFO>
1272 inline void
1273 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(
1274 const DOFINFO &info)
1275 {
1276 const unsigned int level = info.cell->level();
1277
1278 for (unsigned int i = 0; i < matrices->size(); ++i)
1279 {
1280 // Row and column index of
1281 // the block we are dealing with
1282 const unsigned int row = matrices->block(i)[level].row;
1283 const unsigned int col = matrices->block(i)[level].column;
1284
1285 assemble(matrices->block(i)[level].matrix,
1286 info.matrix(i, false).matrix,
1287 row,
1288 col,
1289 info.indices,
1290 info.indices,
1291 level,
1292 level);
1293 if (mg_constrained_dofs != 0)
1294 {
1295 if (interface_in != 0)
1296 assemble_in(interface_in->block(i)[level],
1297 info.matrix(i, false).matrix,
1298 row,
1299 col,
1300 info.indices,
1301 info.indices,
1302 level,
1303 level);
1304 if (interface_out != 0)
1305 assemble_in(interface_out->block(i)[level],
1306 info.matrix(i, false).matrix,
1307 row,
1308 col,
1309 info.indices,
1310 info.indices,
1311 level,
1312 level);
1313
1314 assemble_in(matrices->block_in(i)[level],
1315 info.matrix(i, false).matrix,
1316 row,
1317 col,
1318 info.indices,
1319 info.indices,
1320 level,
1321 level);
1322 assemble_out(matrices->block_out(i)[level],
1323 info.matrix(i, false).matrix,
1324 row,
1325 col,
1326 info.indices,
1327 info.indices,
1328 level,
1329 level);
1330 }
1331 }
1332 }
1333
1334
1335 template <typename MatrixType, typename number>
1336 template <class DOFINFO>
1337 inline void
1338 MGMatrixLocalBlocksToGlobalBlocks<MatrixType, number>::assemble(
1339 const DOFINFO &info1,
1340 const DOFINFO &info2)
1341 {
1342 const unsigned int level1 = info1.cell->level();
1343 const unsigned int level2 = info2.cell->level();
1344
1345 for (unsigned int i = 0; i < matrices->size(); ++i)
1346 {
1347 MGLevelObject<MatrixBlock<MatrixType>> &o = matrices->block(i);
1348
1349 // Row and column index of
1350 // the block we are dealing with
1351 const unsigned int row = o[level1].row;
1352 const unsigned int col = o[level1].column;
1353
1354 if (level1 == level2)
1355 {
1356 if (mg_constrained_dofs == 0)
1357 {
1358 assemble(o[level1].matrix,
1359 info1.matrix(i, false).matrix,
1360 row,
1361 col,
1362 info1.indices,
1363 info1.indices,
1364 level1,
1365 level1);
1366 assemble(o[level1].matrix,
1367 info1.matrix(i, true).matrix,
1368 row,
1369 col,
1370 info1.indices,
1371 info2.indices,
1372 level1,
1373 level2);
1374 assemble(o[level1].matrix,
1375 info2.matrix(i, false).matrix,
1376 row,
1377 col,
1378 info2.indices,
1379 info2.indices,
1380 level2,
1381 level2);
1382 assemble(o[level1].matrix,
1383 info2.matrix(i, true).matrix,
1384 row,
1385 col,
1386 info2.indices,
1387 info1.indices,
1388 level2,
1389 level1);
1390 }
1391 else
1392 {
1393 assemble_fluxes(o[level1].matrix,
1394 info1.matrix(i, false).matrix,
1395 row,
1396 col,
1397 info1.indices,
1398 info1.indices,
1399 level1,
1400 level1);
1401 assemble_fluxes(o[level1].matrix,
1402 info1.matrix(i, true).matrix,
1403 row,
1404 col,
1405 info1.indices,
1406 info2.indices,
1407 level1,
1408 level2);
1409 assemble_fluxes(o[level1].matrix,
1410 info2.matrix(i, false).matrix,
1411 row,
1412 col,
1413 info2.indices,
1414 info2.indices,
1415 level2,
1416 level2);
1417 assemble_fluxes(o[level1].matrix,
1418 info2.matrix(i, true).matrix,
1419 row,
1420 col,
1421 info2.indices,
1422 info1.indices,
1423 level2,
1424 level1);
1425 }
1426 }
1427 else
1428 {
1429 Assert(level1 > level2, ExcNotImplemented());
1430 if (flux_up->size() != 0)
1431 {
1432 // Do not add M22,
1433 // which is done by
1434 // the coarser cell
1435 assemble_fluxes(o[level1].matrix,
1436 info1.matrix(i, false).matrix,
1437 row,
1438 col,
1439 info1.indices,
1440 info1.indices,
1441 level1,
1442 level1);
1443 assemble_up(flux_up->block(i)[level1].matrix,
1444 info1.matrix(i, true).matrix,
1445 row,
1446 col,
1447 info1.indices,
1448 info2.indices,
1449 level1,
1450 level2);
1451 assemble_down(flux_down->block(i)[level1].matrix,
1452 info2.matrix(i, true).matrix,
1453 row,
1454 col,
1455 info2.indices,
1456 info1.indices,
1457 level2,
1458 level1);
1459 }
1460 }
1461 }
1462 }
1463 } // namespace Assembler
1464} // namespace MeshWorker
1465
1466DEAL_II_NAMESPACE_CLOSE
1467
1468#endif
any_data.h
block_vector.h
BlockIndices::block_size
size_type block_size(const unsigned int i) const
Definition block_indices.h:371
BlockIndices::block_start
size_type block_start(const unsigned int i) const
Definition block_indices.h:396
BlockInfo
A small class collecting the different BlockIndices involved in global, multilevel and local computat...
Definition block_info.h:95
BlockVectorBase::n_blocks
unsigned int n_blocks() const
BlockVectorBase::block
BlockType & block(const unsigned int i)
MatrixBlock::add
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
Definition matrix_block.h:674
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::assemble_down
void assemble_down(MatrixType &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition assembler.h:1058
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::initialize_info
void initialize_info(DOFINFO &info, bool face) const
Definition assembler.h:811
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::assemble_in
void assemble_in(MatrixType &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition assembler.h:1121
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::MatrixPtrVector
MGMatrixBlockVector< MatrixType > MatrixPtrVector
Definition assembler.h:327
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::mg_constrained_dofs
SmartPointer< const MGConstrainedDoFs, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > mg_constrained_dofs
Definition assembler.h:514
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::initialize
void initialize(const BlockInfo *block_info, MatrixPtrVector &matrices)
Definition assembler.h:789
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::assemble_up
void assemble_up(MatrixType &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition assembler.h:995
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::assemble_out
void assemble_out(MatrixType &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition assembler.h:1196
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::initialize_interfaces
void initialize_interfaces(MatrixPtrVector &interface_in, MatrixPtrVector &interface_out)
Definition assembler.h:834
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::assemble_fluxes
void assemble_fluxes(MatrixType &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition assembler.h:932
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::block_info
SmartPointer< const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
Definition assembler.h:507
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks::initialize_edge_flux
void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down)
Definition assembler.h:822
MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks::matrices
SmartPointer< MatrixBlockVector< MatrixType >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > matrices
Definition assembler.h:276
MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks::block_info
SmartPointer< const BlockInfo, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
Definition assembler.h:283
MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks::constraints
SmartPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > constraints
Definition assembler.h:290
MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks::initialize_info
void initialize_info(DOFINFO &info, bool face) const
Definition assembler.h:651
MeshWorker::Assembler::MatrixLocalBlocksToGlobalBlocks::initialize
void initialize(const BlockInfo *block_info, MatrixBlockVector< MatrixType > &matrices)
Definition assembler.h:628
MeshWorker::Assembler::ResidualLocalBlocksToGlobalBlocks::initialize_info
void initialize_info(DOFINFO &info, bool face) const
Definition assembler.h:549
MeshWorker::Assembler::ResidualLocalBlocksToGlobalBlocks::initialize
void initialize(const BlockInfo *block_info, AnyData &residuals)
Definition assembler.h:529
MeshWorker::Assembler::ResidualLocalBlocksToGlobalBlocks::block_info
SmartPointer< const BlockInfo, ResidualLocalBlocksToGlobalBlocks< VectorType > > block_info
Definition assembler.h:171
MeshWorker::Assembler::ResidualLocalBlocksToGlobalBlocks::constraints
SmartPointer< const AffineConstraints< typename VectorType::value_type >, ResidualLocalBlocksToGlobalBlocks< VectorType > > constraints
Definition assembler.h:178
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:610
level
unsigned int level
Definition grid_out.cc:4626
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1638
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
dof_info.h
mg_constrained_dofs.h
mg_level_object.h
Generated by doxygen 1.11.0