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