Reference documentation for deal.II version 9.2.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\}}\)
simple.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 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 
23 
26 
28 
31 
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 
41 
42 namespace MeshWorker
43 {
44  namespace Assembler
45  {
59  template <typename VectorType>
61  {
62  public:
69  void
70  initialize(AnyData &results);
71 
75  void
76  initialize(
78 
92  void
94 
102  template <class DOFINFO>
103  void
104  initialize_info(DOFINFO &info, bool face) const;
105 
112  template <class DOFINFO>
113  void
114  assemble(const DOFINFO &info);
115 
119  template <class DOFINFO>
120  void
121  assemble(const DOFINFO &info1, const DOFINFO &info2);
122 
123  protected:
128 
135  };
136 
137 
169  template <typename MatrixType>
171  {
172  public:
177  MatrixSimple(double threshold = 1.e-12);
178 
182  void
183  initialize(MatrixType &m);
184 
188  void
189  initialize(std::vector<MatrixType> &m);
190 
198  void
199  initialize(
201 
209  template <class DOFINFO>
210  void
211  initialize_info(DOFINFO &info, bool face) const;
212 
217  template <class DOFINFO>
218  void
219  assemble(const DOFINFO &info);
220 
225  template <class DOFINFO>
226  void
227  assemble(const DOFINFO &info1, const DOFINFO &info2);
228 
229  protected:
233  std::vector<SmartPointer<MatrixType, MatrixSimple<MatrixType>>> matrix;
234 
240  const double threshold;
241 
242  private:
247  void
248  assemble(const FullMatrix<double> & M,
249  const unsigned int index,
250  const std::vector<types::global_dof_index> &i1,
251  const std::vector<types::global_dof_index> &i2);
252 
259  };
260 
261 
272  template <typename MatrixType>
274  {
275  public:
280  MGMatrixSimple(double threshold = 1.e-12);
281 
285  void
287 
291  void
292  initialize(const MGConstrainedDoFs &mg_constrained_dofs);
293 
298  void
299  initialize_fluxes(MGLevelObject<MatrixType> &flux_up,
300  MGLevelObject<MatrixType> &flux_down);
301 
306  void
307  initialize_interfaces(MGLevelObject<MatrixType> &interface_in,
308  MGLevelObject<MatrixType> &interface_out);
316  template <class DOFINFO>
317  void
318  initialize_info(DOFINFO &info, bool face) const;
319 
323  template <class DOFINFO>
324  void
325  assemble(const DOFINFO &info);
326 
331  template <class DOFINFO>
332  void
333  assemble(const DOFINFO &info1, const DOFINFO &info2);
334 
335  private:
339  void
340  assemble(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 
348  void
349  assemble(MatrixType & G,
350  const FullMatrix<double> & M,
351  const std::vector<types::global_dof_index> &i1,
352  const std::vector<types::global_dof_index> &i2,
353  const unsigned int level);
354 
359  void
360  assemble_up(MatrixType & G,
361  const FullMatrix<double> & M,
362  const std::vector<types::global_dof_index> &i1,
363  const std::vector<types::global_dof_index> &i2,
364  const unsigned int level = numbers::invalid_unsigned_int);
369  void
370  assemble_down(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 
380  void
381  assemble_in(MatrixType & G,
382  const FullMatrix<double> & M,
383  const std::vector<types::global_dof_index> &i1,
384  const std::vector<types::global_dof_index> &i2,
385  const unsigned int level = numbers::invalid_unsigned_int);
386 
391  void
392  assemble_out(MatrixType & G,
393  const FullMatrix<double> & M,
394  const std::vector<types::global_dof_index> &i1,
395  const std::vector<types::global_dof_index> &i2,
396  const unsigned int level = numbers::invalid_unsigned_int);
397 
403 
410 
417 
424 
436 
442  const double threshold;
443  };
444 
445 
456  template <typename MatrixType, typename VectorType>
457  class SystemSimple : private MatrixSimple<MatrixType>,
458  private ResidualSimple<VectorType>
459  {
460  public:
464  SystemSimple(double threshold = 1.e-12);
465 
469  void
470  initialize(MatrixType &m, VectorType &rhs);
471 
479  void
480  initialize(
482 
490  template <class DOFINFO>
491  void
492  initialize_info(DOFINFO &info, bool face) const;
493 
497  template <class DOFINFO>
498  void
499  assemble(const DOFINFO &info);
500 
505  template <class DOFINFO>
506  void
507  assemble(const DOFINFO &info1, const DOFINFO &info2);
508 
509  private:
514  void
515  assemble(const FullMatrix<double> & M,
516  const Vector<double> & vector,
517  const unsigned int index,
518  const std::vector<types::global_dof_index> &indices);
519 
520  void
521  assemble(const FullMatrix<double> & M,
522  const Vector<double> & vector,
523  const unsigned int index,
524  const std::vector<types::global_dof_index> &i1,
525  const std::vector<types::global_dof_index> &i2);
526  };
527 
528 
529  //----------------------------------------------------------------------//
530 
531  template <typename VectorType>
532  inline void
534  {
535  residuals = results;
536  }
537 
538  template <typename VectorType>
539  inline void
542  {
543  constraints = &c;
544  }
545 
546 
547  template <typename MatrixType>
548  inline void
550  {}
551 
552 
553  template <typename VectorType>
554  template <class DOFINFO>
555  inline void
557  {
558  info.initialize_vectors(residuals.size());
559  }
560 
561 
562  template <typename VectorType>
563  template <class DOFINFO>
564  inline void
566  {
567  for (unsigned int k = 0; k < residuals.size(); ++k)
568  {
570  for (unsigned int i = 0; i != info.vector(k).n_blocks(); ++i)
571  {
572  const std::vector<types::global_dof_index> &ldi =
573  info.vector(k).n_blocks() == 1 ? info.indices :
574  info.indices_by_block[i];
575 
576  if (constraints != nullptr)
577  constraints->distribute_local_to_global(info.vector(k).block(i),
578  ldi,
579  *v);
580  else
581  v->add(ldi, info.vector(k).block(i));
582  }
583  }
584  }
585 
586  template <typename VectorType>
587  template <class DOFINFO>
588  inline void
590  const DOFINFO &info2)
591  {
592  assemble(info1);
593  assemble(info2);
594  }
595 
596 
597  //----------------------------------------------------------------------//
598 
599  template <typename MatrixType>
600  inline MatrixSimple<MatrixType>::MatrixSimple(double threshold)
601  : threshold(threshold)
602  {}
603 
604 
605  template <typename MatrixType>
606  inline void
608  {
609  matrix.resize(1);
610  matrix[0] = &m;
611  }
612 
613 
614  template <typename MatrixType>
615  inline void
616  MatrixSimple<MatrixType>::initialize(std::vector<MatrixType> &m)
617  {
618  matrix.resize(m.size());
619  for (unsigned int i = 0; i < m.size(); ++i)
620  matrix[i] = &m[i];
621  }
622 
623 
624  template <typename MatrixType>
625  inline void
628  {
629  constraints = &c;
630  }
631 
632 
633  template <typename MatrixType>
634  template <class DOFINFO>
635  inline void
636  MatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
637  {
638  Assert(matrix.size() != 0, ExcNotInitialized());
639 
640  const unsigned int n = info.indices_by_block.size();
641 
642  if (n == 0)
643  info.initialize_matrices(matrix.size(), face);
644  else
645  {
646  info.initialize_matrices(matrix.size() * n * n, face);
647  unsigned int k = 0;
648  for (unsigned int m = 0; m < matrix.size(); ++m)
649  for (unsigned int i = 0; i < n; ++i)
650  for (unsigned int j = 0; j < n; ++j, ++k)
651  {
652  info.matrix(k, false).row = i;
653  info.matrix(k, false).column = j;
654  if (face)
655  {
656  info.matrix(k, true).row = i;
657  info.matrix(k, true).column = j;
658  }
659  }
660  }
661  }
662 
663 
664 
665  template <typename MatrixType>
666  inline void
668  const FullMatrix<double> & M,
669  const unsigned int index,
670  const std::vector<types::global_dof_index> &i1,
671  const std::vector<types::global_dof_index> &i2)
672  {
673  AssertDimension(M.m(), i1.size());
674  AssertDimension(M.n(), i2.size());
675 
676  if (constraints == nullptr)
677  {
678  for (unsigned int j = 0; j < i1.size(); ++j)
679  for (unsigned int k = 0; k < i2.size(); ++k)
680  if (std::fabs(M(j, k)) >= threshold)
681  matrix[index]->add(i1[j], i2[k], M(j, k));
682  }
683  else
684  constraints->distribute_local_to_global(M, i1, i2, *matrix[index]);
685  }
686 
687 
688  template <typename MatrixType>
689  template <class DOFINFO>
690  inline void
692  {
693  Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
694  const unsigned int n = info.indices_by_block.size();
695 
696  if (n == 0)
697  for (unsigned int m = 0; m < matrix.size(); ++m)
698  assemble(info.matrix(m, false).matrix, m, info.indices, info.indices);
699  else
700  {
701  for (unsigned int m = 0; m < matrix.size(); ++m)
702  for (unsigned int k = 0; k < n * n; ++k)
703  {
704  assemble(
705  info.matrix(k + m * n * n, false).matrix,
706  m,
707  info.indices_by_block[info.matrix(k + m * n * n, false).row],
708  info.indices_by_block[info.matrix(k + m * n * n, false)
709  .column]);
710  }
711  }
712  }
713 
714 
715  template <typename MatrixType>
716  template <class DOFINFO>
717  inline void
719  const DOFINFO &info2)
720  {
721  Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
722  Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
723  AssertDimension(info1.indices_by_block.size(),
724  info2.indices_by_block.size());
725 
726  const unsigned int n = info1.indices_by_block.size();
727 
728  if (n == 0)
729  {
730  for (unsigned int m = 0; m < matrix.size(); ++m)
731  {
732  assemble(info1.matrix(m, false).matrix,
733  m,
734  info1.indices,
735  info1.indices);
736  assemble(info1.matrix(m, true).matrix,
737  m,
738  info1.indices,
739  info2.indices);
740  assemble(info2.matrix(m, false).matrix,
741  m,
742  info2.indices,
743  info2.indices);
744  assemble(info2.matrix(m, true).matrix,
745  m,
746  info2.indices,
747  info1.indices);
748  }
749  }
750  else
751  {
752  for (unsigned int m = 0; m < matrix.size(); ++m)
753  for (unsigned int k = 0; k < n * n; ++k)
754  {
755  const unsigned int row = info1.matrix(k + m * n * n, false).row;
756  const unsigned int column =
757  info1.matrix(k + m * n * n, false).column;
758 
759  assemble(info1.matrix(k + m * n * n, false).matrix,
760  m,
761  info1.indices_by_block[row],
762  info1.indices_by_block[column]);
763  assemble(info1.matrix(k + m * n * n, true).matrix,
764  m,
765  info1.indices_by_block[row],
766  info2.indices_by_block[column]);
767  assemble(info2.matrix(k + m * n * n, false).matrix,
768  m,
769  info2.indices_by_block[row],
770  info2.indices_by_block[column]);
771  assemble(info2.matrix(k + m * n * n, true).matrix,
772  m,
773  info2.indices_by_block[row],
774  info1.indices_by_block[column]);
775  }
776  }
777  }
778 
779 
780  //----------------------------------------------------------------------//
781 
782  template <typename MatrixType>
784  : threshold(threshold)
785  {}
786 
787 
788  template <typename MatrixType>
789  inline void
791  {
792  matrix = &m;
793  }
794 
795  template <typename MatrixType>
796  inline void
798  {
799  mg_constrained_dofs = &c;
800  }
801 
802 
803  template <typename MatrixType>
804  inline void
808  {
809  flux_up = &up;
810  flux_down = &down;
811  }
812 
813 
814  template <typename MatrixType>
815  inline void
819  {
820  interface_in = &in;
821  interface_out = &out;
822  }
823 
824 
825  template <typename MatrixType>
826  template <class DOFINFO>
827  inline void
828  MGMatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
829  {
830  const unsigned int n = info.indices_by_block.size();
831 
832  if (n == 0)
833  info.initialize_matrices(1, face);
834  else
835  {
836  info.initialize_matrices(n * n, face);
837  unsigned int k = 0;
838  for (unsigned int i = 0; i < n; ++i)
839  for (unsigned int j = 0; j < n; ++j, ++k)
840  {
841  info.matrix(k, false).row = i;
842  info.matrix(k, false).column = j;
843  if (face)
844  {
845  info.matrix(k, true).row = i;
846  info.matrix(k, true).column = j;
847  }
848  }
849  }
850  }
851 
852 
853  template <typename MatrixType>
854  inline void
856  MatrixType & G,
857  const FullMatrix<double> & M,
858  const std::vector<types::global_dof_index> &i1,
859  const std::vector<types::global_dof_index> &i2)
860  {
861  AssertDimension(M.m(), i1.size());
862  AssertDimension(M.n(), i2.size());
864  // TODO: Possibly remove this function all together
865 
866  for (unsigned int j = 0; j < i1.size(); ++j)
867  for (unsigned int k = 0; k < i2.size(); ++k)
868  if (std::fabs(M(j, k)) >= threshold)
869  G.add(i1[j], i2[k], M(j, k));
870  }
871 
872 
873  template <typename MatrixType>
874  inline void
876  MatrixType & G,
877  const FullMatrix<double> & M,
878  const std::vector<types::global_dof_index> &i1,
879  const std::vector<types::global_dof_index> &i2,
880  const unsigned int level)
881  {
882  AssertDimension(M.m(), i1.size());
883  AssertDimension(M.n(), i2.size());
884 
885  if (mg_constrained_dofs == nullptr)
886  {
887  for (unsigned int j = 0; j < i1.size(); ++j)
888  for (unsigned int k = 0; k < i2.size(); ++k)
889  if (std::fabs(M(j, k)) >= threshold)
890  G.add(i1[j], i2[k], M(j, k));
891  }
892  else
893  {
894  for (unsigned int j = 0; j < i1.size(); ++j)
895  for (unsigned int k = 0; k < i2.size(); ++k)
896  {
897  // Only enter the local values into the global matrix,
898  // if the value is larger than the threshold
899  if (std::fabs(M(j, k)) < threshold)
900  continue;
901 
902  // Do not enter, if either the row or the column
903  // corresponds to an index on the refinement edge. The
904  // level problems are solved with homogeneous
905  // Dirichlet boundary conditions, therefore we
906  // eliminate these rows and columns. The corresponding
907  // matrix entries are entered by assemble_in() and
908  // assemble_out().
909  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
910  mg_constrained_dofs->at_refinement_edge(level, i2[k]))
911  continue;
912 
913  // At the boundary, only enter the term on the
914  // diagonal, but not the coupling terms
915  if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
916  mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
917  (i1[j] != i2[k]))
918  continue;
919 
920  G.add(i1[j], i2[k], M(j, k));
921  }
922  }
923  }
924 
925 
926  template <typename MatrixType>
927  inline void
929  MatrixType & G,
930  const FullMatrix<double> & M,
931  const std::vector<types::global_dof_index> &i1,
932  const std::vector<types::global_dof_index> &i2,
933  const unsigned int level)
934  {
935  AssertDimension(M.n(), i1.size());
936  AssertDimension(M.m(), i2.size());
937 
938  if (mg_constrained_dofs == nullptr)
939  {
940  for (unsigned int j = 0; j < i1.size(); ++j)
941  for (unsigned int k = 0; k < i2.size(); ++k)
942  if (std::fabs(M(k, j)) >= threshold)
943  G.add(i1[j], i2[k], M(k, j));
944  }
945  else
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(k, j)) >= threshold)
950  if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
951  G.add(i1[j], i2[k], M(k, j));
952  }
953  }
954 
955  template <typename MatrixType>
956  inline void
958  MatrixType & G,
959  const FullMatrix<double> & M,
960  const std::vector<types::global_dof_index> &i1,
961  const std::vector<types::global_dof_index> &i2,
962  const unsigned int level)
963  {
964  AssertDimension(M.m(), i1.size());
965  AssertDimension(M.n(), i2.size());
966 
967  if (mg_constrained_dofs == nullptr)
968  {
969  for (unsigned int j = 0; j < i1.size(); ++j)
970  for (unsigned int k = 0; k < i2.size(); ++k)
971  if (std::fabs(M(j, k)) >= threshold)
972  G.add(i1[j], i2[k], M(j, k));
973  }
974  else
975  {
976  for (unsigned int j = 0; j < i1.size(); ++j)
977  for (unsigned int k = 0; k < i2.size(); ++k)
978  if (std::fabs(M(j, k)) >= threshold)
979  if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
980  G.add(i1[j], i2[k], M(j, k));
981  }
982  }
983 
984  template <typename MatrixType>
985  inline void
987  MatrixType & G,
988  const FullMatrix<double> & M,
989  const std::vector<types::global_dof_index> &i1,
990  const std::vector<types::global_dof_index> &i2,
991  const unsigned int level)
992  {
993  AssertDimension(M.m(), i1.size());
994  AssertDimension(M.n(), i2.size());
996 
997  for (unsigned int j = 0; j < i1.size(); ++j)
998  for (unsigned int k = 0; k < i2.size(); ++k)
999  if (std::fabs(M(j, k)) >= threshold)
1000  // Enter values into matrix only if j corresponds to a
1001  // degree of freedom on the refinement edge, k does
1002  // not, and both are not on the boundary. This is part
1003  // the difference between the complete matrix with no
1004  // boundary condition at the refinement edge and
1005  // the matrix assembled above by assemble().
1006 
1007  // Thus the logic is: enter the row if it is
1008  // constrained by hanging node constraints (actually,
1009  // the whole refinement edge), but not if it is
1010  // constrained by a boundary constraint.
1011  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1012  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1013  {
1014  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1015  !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
1016  (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1017  mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
1018  i1[j] == i2[k]))
1019  G.add(i1[j], i2[k], M(j, k));
1020  }
1021  }
1022 
1023 
1024  template <typename MatrixType>
1025  inline void
1027  MatrixType & G,
1028  const FullMatrix<double> & M,
1029  const std::vector<types::global_dof_index> &i1,
1030  const std::vector<types::global_dof_index> &i2,
1031  const unsigned int level)
1032  {
1033  AssertDimension(M.n(), i1.size());
1034  AssertDimension(M.m(), i2.size());
1036 
1037  for (unsigned int j = 0; j < i1.size(); ++j)
1038  for (unsigned int k = 0; k < i2.size(); ++k)
1039  if (std::fabs(M(k, j)) >= threshold)
1040  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1041  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1042  {
1043  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1044  !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
1045  (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1046  mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
1047  i1[j] == i2[k]))
1048  G.add(i1[j], i2[k], M(k, j));
1049  }
1050  }
1051 
1052 
1053  template <typename MatrixType>
1054  template <class DOFINFO>
1055  inline void
1057  {
1058  Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
1059  const unsigned int level = info.cell->level();
1060 
1061  if (info.indices_by_block.size() == 0)
1062  {
1063  assemble((*matrix)[level],
1064  info.matrix(0, false).matrix,
1065  info.indices,
1066  info.indices,
1067  level);
1068  if (mg_constrained_dofs != nullptr)
1069  {
1070  assemble_in((*interface_in)[level],
1071  info.matrix(0, false).matrix,
1072  info.indices,
1073  info.indices,
1074  level);
1075  assemble_out((*interface_out)[level],
1076  info.matrix(0, false).matrix,
1077  info.indices,
1078  info.indices,
1079  level);
1080  }
1081  }
1082  else
1083  for (unsigned int k = 0; k < info.n_matrices(); ++k)
1084  {
1085  const unsigned int row = info.matrix(k, false).row;
1086  const unsigned int column = info.matrix(k, false).column;
1087 
1088  assemble((*matrix)[level],
1089  info.matrix(k, false).matrix,
1090  info.indices_by_block[row],
1091  info.indices_by_block[column],
1092  level);
1093 
1094  if (mg_constrained_dofs != nullptr)
1095  {
1096  assemble_in((*interface_in)[level],
1097  info.matrix(k, false).matrix,
1098  info.indices_by_block[row],
1099  info.indices_by_block[column],
1100  level);
1101  assemble_out((*interface_out)[level],
1102  info.matrix(k, false).matrix,
1103  info.indices_by_block[column],
1104  info.indices_by_block[row],
1105  level);
1106  }
1107  }
1108  }
1109 
1110 
1111  template <typename MatrixType>
1112  template <class DOFINFO>
1113  inline void
1115  const DOFINFO &info2)
1116  {
1117  Assert(info1.level_cell, ExcMessage("Cell must access level dofs"));
1118  Assert(info2.level_cell, ExcMessage("Cell must access level dofs"));
1119  const unsigned int level1 = info1.cell->level();
1120  const unsigned int level2 = info2.cell->level();
1121 
1122  if (info1.indices_by_block.size() == 0)
1123  {
1124  if (level1 == level2)
1125  {
1126  assemble((*matrix)[level1],
1127  info1.matrix(0, false).matrix,
1128  info1.indices,
1129  info1.indices,
1130  level1);
1131  assemble((*matrix)[level1],
1132  info1.matrix(0, true).matrix,
1133  info1.indices,
1134  info2.indices,
1135  level1);
1136  assemble((*matrix)[level1],
1137  info2.matrix(0, false).matrix,
1138  info2.indices,
1139  info2.indices,
1140  level1);
1141  assemble((*matrix)[level1],
1142  info2.matrix(0, true).matrix,
1143  info2.indices,
1144  info1.indices,
1145  level1);
1146  }
1147  else
1148  {
1149  Assert(level1 > level2, ExcInternalError());
1150  // Do not add info2.M1,
1151  // which is done by
1152  // the coarser cell
1153  assemble((*matrix)[level1],
1154  info1.matrix(0, false).matrix,
1155  info1.indices,
1156  info1.indices,
1157  level1);
1158  if (level1 > 0)
1159  {
1160  assemble_up((*flux_up)[level1],
1161  info1.matrix(0, true).matrix,
1162  info2.indices,
1163  info1.indices,
1164  level1);
1165  assemble_down((*flux_down)[level1],
1166  info2.matrix(0, true).matrix,
1167  info2.indices,
1168  info1.indices,
1169  level1);
1170  }
1171  }
1172  }
1173  else
1174  for (unsigned int k = 0; k < info1.n_matrices(); ++k)
1175  {
1176  const unsigned int row = info1.matrix(k, false).row;
1177  const unsigned int column = info1.matrix(k, false).column;
1178 
1179  if (level1 == level2)
1180  {
1181  assemble((*matrix)[level1],
1182  info1.matrix(k, false).matrix,
1183  info1.indices_by_block[row],
1184  info1.indices_by_block[column],
1185  level1);
1186  assemble((*matrix)[level1],
1187  info1.matrix(k, true).matrix,
1188  info1.indices_by_block[row],
1189  info2.indices_by_block[column],
1190  level1);
1191  assemble((*matrix)[level1],
1192  info2.matrix(k, false).matrix,
1193  info2.indices_by_block[row],
1194  info2.indices_by_block[column],
1195  level1);
1196  assemble((*matrix)[level1],
1197  info2.matrix(k, true).matrix,
1198  info2.indices_by_block[row],
1199  info1.indices_by_block[column],
1200  level1);
1201  }
1202  else
1203  {
1204  Assert(level1 > level2, ExcInternalError());
1205  // Do not add info2.M1,
1206  // which is done by
1207  // the coarser cell
1208  assemble((*matrix)[level1],
1209  info1.matrix(k, false).matrix,
1210  info1.indices_by_block[row],
1211  info1.indices_by_block[column],
1212  level1);
1213  if (level1 > 0)
1214  {
1215  assemble_up((*flux_up)[level1],
1216  info1.matrix(k, true).matrix,
1217  info2.indices_by_block[column],
1218  info1.indices_by_block[row],
1219  level1);
1220  assemble_down((*flux_down)[level1],
1221  info2.matrix(k, true).matrix,
1222  info2.indices_by_block[row],
1223  info1.indices_by_block[column],
1224  level1);
1225  }
1226  }
1227  }
1228  }
1229 
1230  //----------------------------------------------------------------------//
1231 
1232  template <typename MatrixType, typename VectorType>
1234  : MatrixSimple<MatrixType>(t)
1235  {}
1236 
1237 
1238  template <typename MatrixType, typename VectorType>
1239  inline void
1241  VectorType &rhs)
1242  {
1243  AnyData data;
1244  VectorType *p = &rhs;
1245  data.add(p, "right hand side");
1246 
1249  }
1250 
1251  template <typename MatrixType, typename VectorType>
1252  inline void
1255  {
1257  }
1258 
1259 
1260  template <typename MatrixType, typename VectorType>
1261  template <class DOFINFO>
1262  inline void
1264  bool face) const
1265  {
1268  }
1269 
1270  template <typename MatrixType, typename VectorType>
1271  inline void
1273  const FullMatrix<double> & M,
1274  const Vector<double> & vector,
1275  const unsigned int index,
1276  const std::vector<types::global_dof_index> &indices)
1277  {
1278  AssertDimension(M.m(), indices.size());
1279  AssertDimension(M.n(), indices.size());
1280 
1282  VectorType *v = residuals.entry<VectorType *>(index);
1283 
1285  {
1286  for (unsigned int i = 0; i < indices.size(); ++i)
1287  (*v)(indices[i]) += vector(i);
1288 
1289  for (unsigned int j = 0; j < indices.size(); ++j)
1290  for (unsigned int k = 0; k < indices.size(); ++k)
1292  MatrixSimple<MatrixType>::matrix[index]->add(indices[j],
1293  indices[k],
1294  M(j, k));
1295  }
1296  else
1297  {
1298  ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1299  M,
1300  vector,
1301  indices,
1303  *v,
1304  true);
1305  }
1306  }
1307 
1308  template <typename MatrixType, typename VectorType>
1309  inline void
1311  const FullMatrix<double> & M,
1312  const Vector<double> & vector,
1313  const unsigned int index,
1314  const std::vector<types::global_dof_index> &i1,
1315  const std::vector<types::global_dof_index> &i2)
1316  {
1317  AssertDimension(M.m(), i1.size());
1318  AssertDimension(M.n(), i2.size());
1319 
1321  VectorType *v = residuals.entry<VectorType *>(index);
1322 
1324  {
1325  for (unsigned int j = 0; j < i1.size(); ++j)
1326  for (unsigned int k = 0; k < i2.size(); ++k)
1328  MatrixSimple<MatrixType>::matrix[index]->add(i1[j],
1329  i2[k],
1330  M(j, k));
1331  }
1332  else
1333  {
1334  ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1335  vector, i1, i2, *v, M, false);
1336  ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1337  M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
1338  }
1339  }
1340 
1341 
1342  template <typename MatrixType, typename VectorType>
1343  template <class DOFINFO>
1344  inline void
1346  {
1349  Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
1350  const unsigned int n = info.indices_by_block.size();
1351 
1352  if (n == 0)
1353  {
1354  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1355  ++m)
1356  assemble(info.matrix(m, false).matrix,
1357  info.vector(m).block(0),
1358  m,
1359  info.indices);
1360  }
1361  else
1362  {
1363  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1364  ++m)
1365  for (unsigned int k = 0; k < n * n; ++k)
1366  {
1367  const unsigned int row = info.matrix(k + m * n * n, false).row;
1368  const unsigned int column =
1369  info.matrix(k + m * n * n, false).column;
1370 
1371  if (row == column)
1372  assemble(info.matrix(k + m * n * n, false).matrix,
1373  info.vector(m).block(row),
1374  m,
1375  info.indices_by_block[row]);
1376  else
1377  assemble(info.matrix(k + m * n * n, false).matrix,
1378  info.vector(m).block(row),
1379  m,
1380  info.indices_by_block[row],
1381  info.indices_by_block[column]);
1382  }
1383  }
1384  }
1385 
1386 
1387  template <typename MatrixType, typename VectorType>
1388  template <class DOFINFO>
1389  inline void
1391  const DOFINFO &info2)
1392  {
1393  Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
1394  Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
1395  AssertDimension(info1.indices_by_block.size(),
1396  info2.indices_by_block.size());
1397 
1398  const unsigned int n = info1.indices_by_block.size();
1399 
1400  if (n == 0)
1401  {
1402  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1403  ++m)
1404  {
1405  assemble(info1.matrix(m, false).matrix,
1406  info1.vector(m).block(0),
1407  m,
1408  info1.indices);
1409  assemble(info1.matrix(m, true).matrix,
1410  info1.vector(m).block(0),
1411  m,
1412  info1.indices,
1413  info2.indices);
1414  assemble(info2.matrix(m, false).matrix,
1415  info2.vector(m).block(0),
1416  m,
1417  info2.indices);
1418  assemble(info2.matrix(m, true).matrix,
1419  info2.vector(m).block(0),
1420  m,
1421  info2.indices,
1422  info1.indices);
1423  }
1424  }
1425  else
1426  {
1427  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1428  ++m)
1429  for (unsigned int k = 0; k < n * n; ++k)
1430  {
1431  const unsigned int row = info1.matrix(k + m * n * n, false).row;
1432  const unsigned int column =
1433  info1.matrix(k + m * n * n, false).column;
1434 
1435  if (row == column)
1436  {
1437  assemble(info1.matrix(k + m * n * n, false).matrix,
1438  info1.vector(m).block(row),
1439  m,
1440  info1.indices_by_block[row]);
1441  assemble(info2.matrix(k + m * n * n, false).matrix,
1442  info2.vector(m).block(row),
1443  m,
1444  info2.indices_by_block[row]);
1445  }
1446  else
1447  {
1448  assemble(info1.matrix(k + m * n * n, false).matrix,
1449  info1.vector(m).block(row),
1450  m,
1451  info1.indices_by_block[row],
1452  info1.indices_by_block[column]);
1453  assemble(info2.matrix(k + m * n * n, false).matrix,
1454  info2.vector(m).block(row),
1455  m,
1456  info2.indices_by_block[row],
1457  info2.indices_by_block[column]);
1458  }
1459  assemble(info1.matrix(k + m * n * n, true).matrix,
1460  info1.vector(m).block(row),
1461  m,
1462  info1.indices_by_block[row],
1463  info2.indices_by_block[column]);
1464  assemble(info2.matrix(k + m * n * n, true).matrix,
1465  info2.vector(m).block(row),
1466  m,
1467  info2.indices_by_block[row],
1468  info1.indices_by_block[column]);
1469  }
1470  }
1471  }
1472  } // namespace Assembler
1473 } // namespace MeshWorker
1474 
1476 
1477 #endif
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > flux_down
Definition: simple.h:416
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:191
type entry(const std::string &name)
Access to stored data object by name.
Definition: any_data.h:351
void initialize_interfaces(MGLevelObject< MatrixType > &interface_in, MGLevelObject< MatrixType > &interface_out)
Definition: simple.h:816
SystemSimple(double threshold=1.e-12)
Definition: simple.h:1233
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
std::vector< SmartPointer< MatrixType, MatrixSimple< MatrixType > > > matrix
Definition: simple.h:233
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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:1026
void initialize(MGLevelObject< MatrixType > &m)
Definition: simple.h:790
static ::ExceptionBase & ExcNotInitialized()
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > matrix
Definition: simple.h:402
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > interface_out
Definition: simple.h:430
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:928
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:556
void initialize(MatrixType &m, VectorType &rhs)
Definition: simple.h:1240
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:986
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1419
void assemble(const DOFINFO &info)
Definition: simple.h:565
void initialize_fluxes(MGLevelObject< MatrixType > &flux_up, MGLevelObject< MatrixType > &flux_down)
Definition: simple.h:805
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
unsigned int level
Definition: grid_out.cc:4355
void initialize(AnyData &results)
Definition: simple.h:533
MGMatrixSimple(double threshold=1.e-12)
Definition: simple.h:783
SmartPointer< const AffineConstraints< typename VectorType::value_type >, ResidualSimple< VectorType > > constraints
Definition: simple.h:134
Expression fabs(const Expression &x)
void initialize_local_blocks(const BlockIndices &)
Definition: simple.h:549
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:957
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > interface_in
Definition: simple.h:423
void add(type entry, const std::string &name)
Add a new data object.
Definition: any_data.h:432
unsigned int size() const
Number of stored data objects.
Definition: any_data.h:223
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > flux_up
Definition: simple.h:409
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
void assemble(const DOFINFO &info)
Definition: simple.h:1345
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:828
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:1263
void assemble(const DOFINFO &info)
Definition: simple.h:1056
void initialize(MatrixType &m)
Definition: simple.h:607
SmartPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixSimple< MatrixType > > constraints
Definition: simple.h:258
MatrixSimple(double threshold=1.e-12)
Definition: simple.h:600
#define DEAL_II_DEPRECATED
Definition: config.h:98
void assemble(const DOFINFO &info)
Definition: simple.h:691
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:636
static ::ExceptionBase & ExcInternalError()
SmartPointer< const MGConstrainedDoFs, MGMatrixSimple< MatrixType > > mg_constrained_dofs
Definition: simple.h:435