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