Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00:02+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
simple.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2011 - 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_simple_h
17#define dealii_mesh_worker_simple_h
18
19#include <deal.II/base/config.h>
20
22
25
27
30
32
33/*
34 * The header containing the classes MeshWorker::Assembler::MatrixSimple,
35 * MeshWorker::Assembler::MGMatrixSimple, MeshWorker::Assembler::ResidualSimple,
36 * and MeshWorker::Assembler::SystemSimple.
37 */
38
40
41namespace MeshWorker
42{
43 namespace Assembler
44 {
57 template <typename VectorType>
59 {
60 public:
67 void
68 initialize(AnyData &results);
69
73 void
76
84 template <class DOFINFO>
85 void
86 initialize_info(DOFINFO &info, bool face) const;
87
94 template <class DOFINFO>
95 void
96 assemble(const DOFINFO &info);
97
101 template <class DOFINFO>
102 void
103 assemble(const DOFINFO &info1, const DOFINFO &info2);
104
105 protected:
110
117 };
118
119
150 template <typename MatrixType>
152 {
153 public:
158 MatrixSimple(double threshold = 1.e-12);
159
163 void
164 initialize(MatrixType &m);
165
169 void
170 initialize(std::vector<MatrixType> &m);
171
179 void
182
190 template <class DOFINFO>
191 void
192 initialize_info(DOFINFO &info, bool face) const;
193
198 template <class DOFINFO>
199 void
200 assemble(const DOFINFO &info);
201
206 template <class DOFINFO>
207 void
208 assemble(const DOFINFO &info1, const DOFINFO &info2);
209
210 protected:
214 std::vector<SmartPointer<MatrixType, MatrixSimple<MatrixType>>> matrix;
215
221 const double threshold;
222
223 private:
228 void
230 const unsigned int index,
231 const std::vector<types::global_dof_index> &i1,
232 const std::vector<types::global_dof_index> &i2);
233
240 };
241
242
252 template <typename MatrixType>
254 {
255 public:
260 MGMatrixSimple(double threshold = 1.e-12);
261
265 void
267
271 void
273
278 void
281
286 void
296 template <class DOFINFO>
297 void
298 initialize_info(DOFINFO &info, bool face) const;
299
303 template <class DOFINFO>
304 void
305 assemble(const DOFINFO &info);
306
311 template <class DOFINFO>
312 void
313 assemble(const DOFINFO &info1, const DOFINFO &info2);
314
315 private:
319 void
320 assemble(MatrixType &G,
321 const FullMatrix<double> &M,
322 const std::vector<types::global_dof_index> &i1,
323 const std::vector<types::global_dof_index> &i2);
324
328 void
329 assemble(MatrixType &G,
330 const FullMatrix<double> &M,
331 const std::vector<types::global_dof_index> &i1,
332 const std::vector<types::global_dof_index> &i2,
333 const unsigned int level);
334
339 void
340 assemble_up(MatrixType &G,
341 const FullMatrix<double> &M,
342 const std::vector<types::global_dof_index> &i1,
343 const std::vector<types::global_dof_index> &i2,
344 const unsigned int level = numbers::invalid_unsigned_int);
349 void
350 assemble_down(MatrixType &G,
351 const FullMatrix<double> &M,
352 const std::vector<types::global_dof_index> &i1,
353 const std::vector<types::global_dof_index> &i2,
354 const unsigned int level = numbers::invalid_unsigned_int);
355
360 void
361 assemble_in(MatrixType &G,
362 const FullMatrix<double> &M,
363 const std::vector<types::global_dof_index> &i1,
364 const std::vector<types::global_dof_index> &i2,
365 const unsigned int level = numbers::invalid_unsigned_int);
366
371 void
372 assemble_out(MatrixType &G,
373 const FullMatrix<double> &M,
374 const std::vector<types::global_dof_index> &i1,
375 const std::vector<types::global_dof_index> &i2,
376 const unsigned int level = numbers::invalid_unsigned_int);
377
383
390
397
404
416
422 const double threshold;
423 };
424
425
435 template <typename MatrixType, typename VectorType>
436 class SystemSimple : private MatrixSimple<MatrixType>,
437 private ResidualSimple<VectorType>
438 {
439 public:
443 SystemSimple(double threshold = 1.e-12);
444
448 void
449 initialize(MatrixType &m, VectorType &rhs);
450
458 void
461
469 template <class DOFINFO>
470 void
471 initialize_info(DOFINFO &info, bool face) const;
472
476 template <class DOFINFO>
477 void
478 assemble(const DOFINFO &info);
479
484 template <class DOFINFO>
485 void
486 assemble(const DOFINFO &info1, const DOFINFO &info2);
487
488 private:
493 void
495 const Vector<double> &vector,
496 const unsigned int index,
497 const std::vector<types::global_dof_index> &indices);
498
499 void
501 const Vector<double> &vector,
502 const unsigned int index,
503 const std::vector<types::global_dof_index> &i1,
504 const std::vector<types::global_dof_index> &i2);
505 };
506
507
508 //----------------------------------------------------------------------//
509
510 template <typename VectorType>
511 inline void
513 {
514 residuals = results;
515 }
516
517
518
519 template <typename VectorType>
520 inline void
526
527
528
529 template <typename VectorType>
530 template <class DOFINFO>
531 inline void
533 {
534 info.initialize_vectors(residuals.size());
535 }
536
537
538
539 template <typename VectorType>
540 template <class DOFINFO>
541 inline void
543 {
544 for (unsigned int k = 0; k < residuals.size(); ++k)
545 {
546 VectorType *v = residuals.entry<VectorType *>(k);
547 for (unsigned int i = 0; i != info.vector(k).n_blocks(); ++i)
548 {
549 const std::vector<types::global_dof_index> &ldi =
550 info.vector(k).n_blocks() == 1 ? info.indices :
551 info.indices_by_block[i];
552
553 if (constraints != nullptr)
554 constraints->distribute_local_to_global(info.vector(k).block(i),
555 ldi,
556 *v);
557 else
558 v->add(ldi, info.vector(k).block(i));
559 }
560 }
561 }
562
563 template <typename VectorType>
564 template <class DOFINFO>
565 inline void
567 const DOFINFO &info2)
568 {
569 assemble(info1);
570 assemble(info2);
571 }
572
573
574 //----------------------------------------------------------------------//
575
576 template <typename MatrixType>
578 : threshold(threshold)
579 {}
580
581
582 template <typename MatrixType>
583 inline void
585 {
586 matrix.resize(1);
587 matrix[0] = &m;
588 }
589
590
591 template <typename MatrixType>
592 inline void
593 MatrixSimple<MatrixType>::initialize(std::vector<MatrixType> &m)
594 {
595 matrix.resize(m.size());
596 for (unsigned int i = 0; i < m.size(); ++i)
597 matrix[i] = &m[i];
598 }
599
600
601 template <typename MatrixType>
602 inline void
608
609
610 template <typename MatrixType>
611 template <class DOFINFO>
612 inline void
613 MatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
614 {
615 Assert(matrix.size() != 0, ExcNotInitialized());
616
617 const unsigned int n = info.indices_by_block.size();
618
619 if (n == 0)
620 info.initialize_matrices(matrix.size(), face);
621 else
622 {
623 info.initialize_matrices(matrix.size() * n * n, face);
624 unsigned int k = 0;
625 for (unsigned int m = 0; m < matrix.size(); ++m)
626 for (unsigned int i = 0; i < n; ++i)
627 for (unsigned int j = 0; j < n; ++j, ++k)
628 {
629 info.matrix(k, false).row = i;
630 info.matrix(k, false).column = j;
631 if (face)
632 {
633 info.matrix(k, true).row = i;
634 info.matrix(k, true).column = j;
635 }
636 }
637 }
638 }
639
640
641
642 template <typename MatrixType>
643 inline void
645 const FullMatrix<double> &M,
646 const unsigned int index,
647 const std::vector<types::global_dof_index> &i1,
648 const std::vector<types::global_dof_index> &i2)
649 {
650 AssertDimension(M.m(), i1.size());
651 AssertDimension(M.n(), i2.size());
652
653 if (constraints == nullptr)
654 {
655 for (unsigned int j = 0; j < i1.size(); ++j)
656 for (unsigned int k = 0; k < i2.size(); ++k)
657 if (std::fabs(M(j, k)) >= threshold)
658 matrix[index]->add(i1[j], i2[k], M(j, k));
659 }
660 else
661 constraints->distribute_local_to_global(M, i1, i2, *matrix[index]);
662 }
663
664
665 template <typename MatrixType>
666 template <class DOFINFO>
667 inline void
669 {
670 Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
671 const unsigned int n = info.indices_by_block.size();
672
673 if (n == 0)
674 for (unsigned int m = 0; m < matrix.size(); ++m)
675 assemble(info.matrix(m, false).matrix, m, info.indices, info.indices);
676 else
677 {
678 for (unsigned int m = 0; m < matrix.size(); ++m)
679 for (unsigned int k = 0; k < n * n; ++k)
680 {
681 assemble(
682 info.matrix(k + m * n * n, false).matrix,
683 m,
684 info.indices_by_block[info.matrix(k + m * n * n, false).row],
685 info.indices_by_block[info.matrix(k + m * n * n, false)
686 .column]);
687 }
688 }
689 }
690
691
692 template <typename MatrixType>
693 template <class DOFINFO>
694 inline void
696 const DOFINFO &info2)
697 {
698 Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
699 Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
700 AssertDimension(info1.indices_by_block.size(),
701 info2.indices_by_block.size());
702
703 const unsigned int n = info1.indices_by_block.size();
704
705 if (n == 0)
706 {
707 for (unsigned int m = 0; m < matrix.size(); ++m)
708 {
709 assemble(info1.matrix(m, false).matrix,
710 m,
711 info1.indices,
712 info1.indices);
713 assemble(info1.matrix(m, true).matrix,
714 m,
715 info1.indices,
716 info2.indices);
717 assemble(info2.matrix(m, false).matrix,
718 m,
719 info2.indices,
720 info2.indices);
721 assemble(info2.matrix(m, true).matrix,
722 m,
723 info2.indices,
724 info1.indices);
725 }
726 }
727 else
728 {
729 for (unsigned int m = 0; m < matrix.size(); ++m)
730 for (unsigned int k = 0; k < n * n; ++k)
731 {
732 const unsigned int row = info1.matrix(k + m * n * n, false).row;
733 const unsigned int column =
734 info1.matrix(k + m * n * n, false).column;
735
736 assemble(info1.matrix(k + m * n * n, false).matrix,
737 m,
738 info1.indices_by_block[row],
739 info1.indices_by_block[column]);
740 assemble(info1.matrix(k + m * n * n, true).matrix,
741 m,
742 info1.indices_by_block[row],
743 info2.indices_by_block[column]);
744 assemble(info2.matrix(k + m * n * n, false).matrix,
745 m,
746 info2.indices_by_block[row],
747 info2.indices_by_block[column]);
748 assemble(info2.matrix(k + m * n * n, true).matrix,
749 m,
750 info2.indices_by_block[row],
751 info1.indices_by_block[column]);
752 }
753 }
754 }
755
756
757 //----------------------------------------------------------------------//
758
759 template <typename MatrixType>
761 : threshold(threshold)
762 {}
763
764
765 template <typename MatrixType>
766 inline void
771
772 template <typename MatrixType>
773 inline void
775 {
776 mg_constrained_dofs = &c;
777 }
778
779
780 template <typename MatrixType>
781 inline void
785 {
786 flux_up = &up;
787 flux_down = &down;
788 }
789
790
791 template <typename MatrixType>
792 inline void
796 {
797 interface_in = &in;
798 interface_out = &out;
799 }
800
801
802 template <typename MatrixType>
803 template <class DOFINFO>
804 inline void
805 MGMatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
806 {
807 const unsigned int n = info.indices_by_block.size();
808
809 if (n == 0)
810 info.initialize_matrices(1, face);
811 else
812 {
813 info.initialize_matrices(n * n, face);
814 unsigned int k = 0;
815 for (unsigned int i = 0; i < n; ++i)
816 for (unsigned int j = 0; j < n; ++j, ++k)
817 {
818 info.matrix(k, false).row = i;
819 info.matrix(k, false).column = j;
820 if (face)
821 {
822 info.matrix(k, true).row = i;
823 info.matrix(k, true).column = j;
824 }
825 }
826 }
827 }
828
829
830 template <typename MatrixType>
831 inline void
833 MatrixType &G,
834 const FullMatrix<double> &M,
835 const std::vector<types::global_dof_index> &i1,
836 const std::vector<types::global_dof_index> &i2)
837 {
838 AssertDimension(M.m(), i1.size());
839 AssertDimension(M.n(), i2.size());
840 Assert(mg_constrained_dofs == 0, ExcInternalError());
841 // TODO: Possibly remove this function all together
842
843 for (unsigned int j = 0; j < i1.size(); ++j)
844 for (unsigned int k = 0; k < i2.size(); ++k)
845 if (std::fabs(M(j, k)) >= threshold)
846 G.add(i1[j], i2[k], M(j, k));
847 }
848
849
850 template <typename MatrixType>
851 inline void
853 MatrixType &G,
854 const FullMatrix<double> &M,
855 const std::vector<types::global_dof_index> &i1,
856 const std::vector<types::global_dof_index> &i2,
857 const unsigned int level)
858 {
859 AssertDimension(M.m(), i1.size());
860 AssertDimension(M.n(), i2.size());
861
862 if (mg_constrained_dofs == nullptr)
863 {
864 for (unsigned int j = 0; j < i1.size(); ++j)
865 for (unsigned int k = 0; k < i2.size(); ++k)
866 if (std::fabs(M(j, k)) >= threshold)
867 G.add(i1[j], i2[k], M(j, k));
868 }
869 else
870 {
871 for (unsigned int j = 0; j < i1.size(); ++j)
872 for (unsigned int k = 0; k < i2.size(); ++k)
873 {
874 // Only enter the local values into the global matrix,
875 // if the value is larger than the threshold
876 if (std::fabs(M(j, k)) < threshold)
877 continue;
878
879 // Do not enter, if either the row or the column
880 // corresponds to an index on the refinement edge. The
881 // level problems are solved with homogeneous
882 // Dirichlet boundary conditions, therefore we
883 // eliminate these rows and columns. The corresponding
884 // matrix entries are entered by assemble_in() and
885 // assemble_out().
886 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
887 mg_constrained_dofs->at_refinement_edge(level, i2[k]))
888 continue;
889
890 // At the boundary, only enter the term on the
891 // diagonal, but not the coupling terms
892 if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
893 mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
894 (i1[j] != i2[k]))
895 continue;
896
897 G.add(i1[j], i2[k], M(j, k));
898 }
899 }
900 }
901
902
903 template <typename MatrixType>
904 inline void
906 MatrixType &G,
907 const FullMatrix<double> &M,
908 const std::vector<types::global_dof_index> &i1,
909 const std::vector<types::global_dof_index> &i2,
910 const unsigned int level)
911 {
912 AssertDimension(M.n(), i1.size());
913 AssertDimension(M.m(), i2.size());
914
915 if (mg_constrained_dofs == nullptr)
916 {
917 for (unsigned int j = 0; j < i1.size(); ++j)
918 for (unsigned int k = 0; k < i2.size(); ++k)
919 if (std::fabs(M(k, j)) >= threshold)
920 G.add(i1[j], i2[k], M(k, j));
921 }
922 else
923 {
924 for (unsigned int j = 0; j < i1.size(); ++j)
925 for (unsigned int k = 0; k < i2.size(); ++k)
926 if (std::fabs(M(k, j)) >= threshold)
927 if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
928 G.add(i1[j], i2[k], M(k, j));
929 }
930 }
931
932 template <typename MatrixType>
933 inline void
935 MatrixType &G,
936 const FullMatrix<double> &M,
937 const std::vector<types::global_dof_index> &i1,
938 const std::vector<types::global_dof_index> &i2,
939 const unsigned int level)
940 {
941 AssertDimension(M.m(), i1.size());
942 AssertDimension(M.n(), i2.size());
943
944 if (mg_constrained_dofs == nullptr)
945 {
946 for (unsigned int j = 0; j < i1.size(); ++j)
947 for (unsigned int k = 0; k < i2.size(); ++k)
948 if (std::fabs(M(j, k)) >= threshold)
949 G.add(i1[j], i2[k], M(j, k));
950 }
951 else
952 {
953 for (unsigned int j = 0; j < i1.size(); ++j)
954 for (unsigned int k = 0; k < i2.size(); ++k)
955 if (std::fabs(M(j, k)) >= threshold)
956 if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
957 G.add(i1[j], i2[k], M(j, k));
958 }
959 }
960
961 template <typename MatrixType>
962 inline void
964 MatrixType &G,
965 const FullMatrix<double> &M,
966 const std::vector<types::global_dof_index> &i1,
967 const std::vector<types::global_dof_index> &i2,
968 const unsigned int level)
969 {
970 AssertDimension(M.m(), i1.size());
971 AssertDimension(M.n(), i2.size());
972 Assert(mg_constrained_dofs != nullptr, ExcInternalError());
973
974 for (unsigned int j = 0; j < i1.size(); ++j)
975 for (unsigned int k = 0; k < i2.size(); ++k)
976 if (std::fabs(M(j, k)) >= threshold)
977 // Enter values into matrix only if j corresponds to a
978 // degree of freedom on the refinement edge, k does
979 // not, and both are not on the boundary. This is part
980 // the difference between the complete matrix with no
981 // boundary condition at the refinement edge and
982 // the matrix assembled above by assemble().
983
984 // Thus the logic is: enter the row if it is
985 // constrained by hanging node constraints (actually,
986 // the whole refinement edge), but not if it is
987 // constrained by a boundary constraint.
988 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
989 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
990 {
991 if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
992 !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
993 (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
994 mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
995 i1[j] == i2[k]))
996 G.add(i1[j], i2[k], M(j, k));
997 }
998 }
999
1000
1001 template <typename MatrixType>
1002 inline void
1004 MatrixType &G,
1005 const FullMatrix<double> &M,
1006 const std::vector<types::global_dof_index> &i1,
1007 const std::vector<types::global_dof_index> &i2,
1008 const unsigned int level)
1009 {
1010 AssertDimension(M.n(), i1.size());
1011 AssertDimension(M.m(), i2.size());
1012 Assert(mg_constrained_dofs != nullptr, ExcInternalError());
1013
1014 for (unsigned int j = 0; j < i1.size(); ++j)
1015 for (unsigned int k = 0; k < i2.size(); ++k)
1016 if (std::fabs(M(k, j)) >= threshold)
1017 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1018 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1019 {
1020 if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1021 !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
1022 (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1023 mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
1024 i1[j] == i2[k]))
1025 G.add(i1[j], i2[k], M(k, j));
1026 }
1027 }
1028
1029
1030 template <typename MatrixType>
1031 template <class DOFINFO>
1032 inline void
1034 {
1035 Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
1036 const unsigned int level = info.cell->level();
1037
1038 if (info.indices_by_block.empty())
1039 {
1040 assemble((*matrix)[level],
1041 info.matrix(0, false).matrix,
1042 info.indices,
1043 info.indices,
1044 level);
1045 if (mg_constrained_dofs != nullptr)
1046 {
1047 assemble_in((*interface_in)[level],
1048 info.matrix(0, false).matrix,
1049 info.indices,
1050 info.indices,
1051 level);
1052 assemble_out((*interface_out)[level],
1053 info.matrix(0, false).matrix,
1054 info.indices,
1055 info.indices,
1056 level);
1057 }
1058 }
1059 else
1060 for (unsigned int k = 0; k < info.n_matrices(); ++k)
1061 {
1062 const unsigned int row = info.matrix(k, false).row;
1063 const unsigned int column = info.matrix(k, false).column;
1064
1065 assemble((*matrix)[level],
1066 info.matrix(k, false).matrix,
1067 info.indices_by_block[row],
1068 info.indices_by_block[column],
1069 level);
1070
1071 if (mg_constrained_dofs != nullptr)
1072 {
1073 assemble_in((*interface_in)[level],
1074 info.matrix(k, false).matrix,
1075 info.indices_by_block[row],
1076 info.indices_by_block[column],
1077 level);
1078 assemble_out((*interface_out)[level],
1079 info.matrix(k, false).matrix,
1080 info.indices_by_block[column],
1081 info.indices_by_block[row],
1082 level);
1083 }
1084 }
1085 }
1086
1087
1088 template <typename MatrixType>
1089 template <class DOFINFO>
1090 inline void
1092 const DOFINFO &info2)
1093 {
1094 Assert(info1.level_cell, ExcMessage("Cell must access level dofs"));
1095 Assert(info2.level_cell, ExcMessage("Cell must access level dofs"));
1096 const unsigned int level1 = info1.cell->level();
1097 const unsigned int level2 = info2.cell->level();
1098
1099 if (info1.indices_by_block.empty())
1100 {
1101 if (level1 == level2)
1102 {
1103 assemble((*matrix)[level1],
1104 info1.matrix(0, false).matrix,
1105 info1.indices,
1106 info1.indices,
1107 level1);
1108 assemble((*matrix)[level1],
1109 info1.matrix(0, true).matrix,
1110 info1.indices,
1111 info2.indices,
1112 level1);
1113 assemble((*matrix)[level1],
1114 info2.matrix(0, false).matrix,
1115 info2.indices,
1116 info2.indices,
1117 level1);
1118 assemble((*matrix)[level1],
1119 info2.matrix(0, true).matrix,
1120 info2.indices,
1121 info1.indices,
1122 level1);
1123 }
1124 else
1125 {
1126 Assert(level1 > level2, ExcInternalError());
1127 // Do not add info2.M1,
1128 // which is done by
1129 // the coarser cell
1130 assemble((*matrix)[level1],
1131 info1.matrix(0, false).matrix,
1132 info1.indices,
1133 info1.indices,
1134 level1);
1135 if (level1 > 0)
1136 {
1137 assemble_up((*flux_up)[level1],
1138 info1.matrix(0, true).matrix,
1139 info2.indices,
1140 info1.indices,
1141 level1);
1142 assemble_down((*flux_down)[level1],
1143 info2.matrix(0, true).matrix,
1144 info2.indices,
1145 info1.indices,
1146 level1);
1147 }
1148 }
1149 }
1150 else
1151 for (unsigned int k = 0; k < info1.n_matrices(); ++k)
1152 {
1153 const unsigned int row = info1.matrix(k, false).row;
1154 const unsigned int column = info1.matrix(k, false).column;
1155
1156 if (level1 == level2)
1157 {
1158 assemble((*matrix)[level1],
1159 info1.matrix(k, false).matrix,
1160 info1.indices_by_block[row],
1161 info1.indices_by_block[column],
1162 level1);
1163 assemble((*matrix)[level1],
1164 info1.matrix(k, true).matrix,
1165 info1.indices_by_block[row],
1166 info2.indices_by_block[column],
1167 level1);
1168 assemble((*matrix)[level1],
1169 info2.matrix(k, false).matrix,
1170 info2.indices_by_block[row],
1171 info2.indices_by_block[column],
1172 level1);
1173 assemble((*matrix)[level1],
1174 info2.matrix(k, true).matrix,
1175 info2.indices_by_block[row],
1176 info1.indices_by_block[column],
1177 level1);
1178 }
1179 else
1180 {
1181 Assert(level1 > level2, ExcInternalError());
1182 // Do not add info2.M1,
1183 // which is done by
1184 // the coarser cell
1185 assemble((*matrix)[level1],
1186 info1.matrix(k, false).matrix,
1187 info1.indices_by_block[row],
1188 info1.indices_by_block[column],
1189 level1);
1190 if (level1 > 0)
1191 {
1192 assemble_up((*flux_up)[level1],
1193 info1.matrix(k, true).matrix,
1194 info2.indices_by_block[column],
1195 info1.indices_by_block[row],
1196 level1);
1197 assemble_down((*flux_down)[level1],
1198 info2.matrix(k, true).matrix,
1199 info2.indices_by_block[row],
1200 info1.indices_by_block[column],
1201 level1);
1202 }
1203 }
1204 }
1205 }
1206
1207 //----------------------------------------------------------------------//
1208
1209 template <typename MatrixType, typename VectorType>
1213
1214
1215 template <typename MatrixType, typename VectorType>
1216 inline void
1218 VectorType &rhs)
1219 {
1220 AnyData data;
1221 VectorType *p = &rhs;
1222 data.add(p, "right hand side");
1223
1226 }
1227
1228 template <typename MatrixType, typename VectorType>
1229 inline void
1235
1236
1237 template <typename MatrixType, typename VectorType>
1238 template <class DOFINFO>
1239 inline void
1246
1247 template <typename MatrixType, typename VectorType>
1248 inline void
1250 const FullMatrix<double> &M,
1251 const Vector<double> &vector,
1252 const unsigned int index,
1253 const std::vector<types::global_dof_index> &indices)
1254 {
1255 AssertDimension(M.m(), indices.size());
1256 AssertDimension(M.n(), indices.size());
1257
1259 VectorType *v = residuals.entry<VectorType *>(index);
1260
1262 {
1263 for (unsigned int i = 0; i < indices.size(); ++i)
1264 (*v)(indices[i]) += vector(i);
1265
1266 for (unsigned int j = 0; j < indices.size(); ++j)
1267 for (unsigned int k = 0; k < indices.size(); ++k)
1268 if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
1269 MatrixSimple<MatrixType>::matrix[index]->add(indices[j],
1270 indices[k],
1271 M(j, k));
1272 }
1273 else
1274 {
1275 ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1276 M,
1277 vector,
1278 indices,
1280 *v,
1281 true);
1282 }
1283 }
1284
1285 template <typename MatrixType, typename VectorType>
1286 inline void
1288 const FullMatrix<double> &M,
1289 const Vector<double> &vector,
1290 const unsigned int index,
1291 const std::vector<types::global_dof_index> &i1,
1292 const std::vector<types::global_dof_index> &i2)
1293 {
1294 AssertDimension(M.m(), i1.size());
1295 AssertDimension(M.n(), i2.size());
1296
1298 VectorType *v = residuals.entry<VectorType *>(index);
1299
1301 {
1302 for (unsigned int j = 0; j < i1.size(); ++j)
1303 for (unsigned int k = 0; k < i2.size(); ++k)
1304 if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
1305 MatrixSimple<MatrixType>::matrix[index]->add(i1[j],
1306 i2[k],
1307 M(j, k));
1308 }
1309 else
1310 {
1311 ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1312 vector, i1, i2, *v, M, false);
1313 ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1314 M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
1315 }
1316 }
1317
1318
1319 template <typename MatrixType, typename VectorType>
1320 template <class DOFINFO>
1321 inline void
1323 {
1326 Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
1327 const unsigned int n = info.indices_by_block.size();
1328
1329 if (n == 0)
1330 {
1331 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1332 ++m)
1333 assemble(info.matrix(m, false).matrix,
1334 info.vector(m).block(0),
1335 m,
1336 info.indices);
1337 }
1338 else
1339 {
1340 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1341 ++m)
1342 for (unsigned int k = 0; k < n * n; ++k)
1343 {
1344 const unsigned int row = info.matrix(k + m * n * n, false).row;
1345 const unsigned int column =
1346 info.matrix(k + m * n * n, false).column;
1347
1348 if (row == column)
1349 assemble(info.matrix(k + m * n * n, false).matrix,
1350 info.vector(m).block(row),
1351 m,
1352 info.indices_by_block[row]);
1353 else
1354 assemble(info.matrix(k + m * n * n, false).matrix,
1355 info.vector(m).block(row),
1356 m,
1357 info.indices_by_block[row],
1358 info.indices_by_block[column]);
1359 }
1360 }
1361 }
1362
1363
1364 template <typename MatrixType, typename VectorType>
1365 template <class DOFINFO>
1366 inline void
1368 const DOFINFO &info2)
1369 {
1370 Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
1371 Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
1372 AssertDimension(info1.indices_by_block.size(),
1373 info2.indices_by_block.size());
1374
1375 const unsigned int n = info1.indices_by_block.size();
1376
1377 if (n == 0)
1378 {
1379 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1380 ++m)
1381 {
1382 assemble(info1.matrix(m, false).matrix,
1383 info1.vector(m).block(0),
1384 m,
1385 info1.indices);
1386 assemble(info1.matrix(m, true).matrix,
1387 info1.vector(m).block(0),
1388 m,
1389 info1.indices,
1390 info2.indices);
1391 assemble(info2.matrix(m, false).matrix,
1392 info2.vector(m).block(0),
1393 m,
1394 info2.indices);
1395 assemble(info2.matrix(m, true).matrix,
1396 info2.vector(m).block(0),
1397 m,
1398 info2.indices,
1399 info1.indices);
1400 }
1401 }
1402 else
1403 {
1404 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1405 ++m)
1406 for (unsigned int k = 0; k < n * n; ++k)
1407 {
1408 const unsigned int row = info1.matrix(k + m * n * n, false).row;
1409 const unsigned int column =
1410 info1.matrix(k + m * n * n, false).column;
1411
1412 if (row == column)
1413 {
1414 assemble(info1.matrix(k + m * n * n, false).matrix,
1415 info1.vector(m).block(row),
1416 m,
1417 info1.indices_by_block[row]);
1418 assemble(info2.matrix(k + m * n * n, false).matrix,
1419 info2.vector(m).block(row),
1420 m,
1421 info2.indices_by_block[row]);
1422 }
1423 else
1424 {
1425 assemble(info1.matrix(k + m * n * n, false).matrix,
1426 info1.vector(m).block(row),
1427 m,
1428 info1.indices_by_block[row],
1429 info1.indices_by_block[column]);
1430 assemble(info2.matrix(k + m * n * n, false).matrix,
1431 info2.vector(m).block(row),
1432 m,
1433 info2.indices_by_block[row],
1434 info2.indices_by_block[column]);
1435 }
1436 assemble(info1.matrix(k + m * n * n, true).matrix,
1437 info1.vector(m).block(row),
1438 m,
1439 info1.indices_by_block[row],
1440 info2.indices_by_block[column]);
1441 assemble(info2.matrix(k + m * n * n, true).matrix,
1442 info2.vector(m).block(row),
1443 m,
1444 info2.indices_by_block[row],
1445 info1.indices_by_block[column]);
1446 }
1447 }
1448 }
1449 } // namespace Assembler
1450} // namespace MeshWorker
1451
1453
1454#endif
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
type entry(const std::string &name)
Access to stored data object by name.
Definition any_data.h:349
void add(type entry, const std::string &name)
Add a new data object.
Definition any_data.h:430
size_type n() const
size_type m() const
void assemble_out(MatrixType &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
Definition simple.h:1003
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > flux_down
Definition simple.h:396
void initialize_fluxes(MGLevelObject< MatrixType > &flux_up, MGLevelObject< MatrixType > &flux_down)
Definition simple.h:782
void initialize_interfaces(MGLevelObject< MatrixType > &interface_in, MGLevelObject< MatrixType > &interface_out)
Definition simple.h:793
void assemble(const DOFINFO &info)
Definition simple.h:1033
void initialize_info(DOFINFO &info, bool face) const
Definition simple.h:805
void assemble_in(MatrixType &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
Definition simple.h:963
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > interface_out
Definition simple.h:410
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > matrix
Definition simple.h:382
MGMatrixSimple(double threshold=1.e-12)
Definition simple.h:760
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > interface_in
Definition simple.h:403
void initialize(MGLevelObject< MatrixType > &m)
Definition simple.h:767
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > flux_up
Definition simple.h:389
void assemble_down(MatrixType &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
Definition simple.h:934
void assemble_up(MatrixType &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
Definition simple.h:905
SmartPointer< const MGConstrainedDoFs, MGMatrixSimple< MatrixType > > mg_constrained_dofs
Definition simple.h:415
SmartPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixSimple< MatrixType > > constraints
Definition simple.h:239
void initialize(MatrixType &m)
Definition simple.h:584
MatrixSimple(double threshold=1.e-12)
Definition simple.h:577
void initialize_info(DOFINFO &info, bool face) const
Definition simple.h:613
void assemble(const DOFINFO &info)
Definition simple.h:668
std::vector< SmartPointer< MatrixType, MatrixSimple< MatrixType > > > matrix
Definition simple.h:214
void assemble(const DOFINFO &info)
Definition simple.h:542
void initialize_info(DOFINFO &info, bool face) const
Definition simple.h:532
SmartPointer< const AffineConstraints< typename VectorType::value_type >, ResidualSimple< VectorType > > constraints
Definition simple.h:116
void initialize(AnyData &results)
Definition simple.h:512
SystemSimple(double threshold=1.e-12)
Definition simple.h:1210
void assemble(const DOFINFO &info)
Definition simple.h:1322
void initialize(MatrixType &m, VectorType &rhs)
Definition simple.h:1217
void initialize_info(DOFINFO &info, bool face) const
Definition simple.h:1240
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
unsigned int level
Definition grid_out.cc:4626
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
static const unsigned int invalid_unsigned_int
Definition types.h:220