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