Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
simple.h
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/algorithms/any_data.h>
21 
22 #include <deal.II/base/mg_level_object.h>
23 #include <deal.II/base/smartpointer.h>
24 
25 #include <deal.II/lac/block_vector.h>
26 
27 #include <deal.II/meshworker/dof_info.h>
28 #include <deal.II/meshworker/functional.h>
29 
30 #include <deal.II/multigrid/mg_constrained_dofs.h>
31 
32 /*
33  * The header containing the classes MeshWorker::Assembler::MatrixSimple,
34  * MeshWorker::Assembler::MGMatrixSimple, MeshWorker::Assembler::ResidualSimple,
35  * and MeshWorker::Assembler::SystemSimple.
36  */
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 namespace MeshWorker
41 {
42  namespace Assembler
43  {
57  template <typename VectorType>
59  {
60  public:
67  void
68  initialize(AnyData &results);
69 
73  void
74  initialize(
76 
89  DEAL_II_DEPRECATED
90  void
92 
100  template <class DOFINFO>
101  void
102  initialize_info(DOFINFO &info, bool face) const;
103 
110  template <class DOFINFO>
111  void
112  assemble(const DOFINFO &info);
113 
117  template <class DOFINFO>
118  void
119  assemble(const DOFINFO &info1, const DOFINFO &info2);
120 
121  protected:
126 
133  };
134 
135 
167  template <typename MatrixType>
169  {
170  public:
175  MatrixSimple(double threshold = 1.e-12);
176 
180  void
181  initialize(MatrixType &m);
182 
186  void
187  initialize(std::vector<MatrixType> &m);
188 
196  void
197  initialize(
199 
207  template <class DOFINFO>
208  void
209  initialize_info(DOFINFO &info, bool face) const;
210 
215  template <class DOFINFO>
216  void
217  assemble(const DOFINFO &info);
218 
223  template <class DOFINFO>
224  void
225  assemble(const DOFINFO &info1, const DOFINFO &info2);
226 
227  protected:
231  std::vector<SmartPointer<MatrixType, MatrixSimple<MatrixType>>> matrix;
232 
238  const double threshold;
239 
240  private:
245  void
246  assemble(const FullMatrix<double> & M,
247  const unsigned int index,
248  const std::vector<types::global_dof_index> &i1,
249  const std::vector<types::global_dof_index> &i2);
250 
257  };
258 
259 
270  template <typename MatrixType>
272  {
273  public:
278  MGMatrixSimple(double threshold = 1.e-12);
279 
283  void
285 
289  void
291 
296  void
299 
304  void
314  template <class DOFINFO>
315  void
316  initialize_info(DOFINFO &info, bool face) const;
317 
321  template <class DOFINFO>
322  void
323  assemble(const DOFINFO &info);
324 
329  template <class DOFINFO>
330  void
331  assemble(const DOFINFO &info1, const DOFINFO &info2);
332 
333  private:
337  void
338  assemble(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 
346  void
347  assemble(MatrixType & G,
348  const FullMatrix<double> & M,
349  const std::vector<types::global_dof_index> &i1,
350  const std::vector<types::global_dof_index> &i2,
351  const unsigned int level);
352 
357  void
358  assemble_up(MatrixType & G,
359  const FullMatrix<double> & M,
360  const std::vector<types::global_dof_index> &i1,
361  const std::vector<types::global_dof_index> &i2,
362  const unsigned int level = numbers::invalid_unsigned_int);
367  void
368  assemble_down(MatrixType & G,
369  const FullMatrix<double> & M,
370  const std::vector<types::global_dof_index> &i1,
371  const std::vector<types::global_dof_index> &i2,
372  const unsigned int level = numbers::invalid_unsigned_int);
373 
378  void
379  assemble_in(MatrixType & G,
380  const FullMatrix<double> & M,
381  const std::vector<types::global_dof_index> &i1,
382  const std::vector<types::global_dof_index> &i2,
383  const unsigned int level = numbers::invalid_unsigned_int);
384 
389  void
390  assemble_out(MatrixType & G,
391  const FullMatrix<double> & M,
392  const std::vector<types::global_dof_index> &i1,
393  const std::vector<types::global_dof_index> &i2,
394  const unsigned int level = numbers::invalid_unsigned_int);
395 
401 
408 
415 
422 
434 
440  const double threshold;
441  };
442 
443 
454  template <typename MatrixType, typename VectorType>
455  class SystemSimple : private MatrixSimple<MatrixType>,
456  private ResidualSimple<VectorType>
457  {
458  public:
462  SystemSimple(double threshold = 1.e-12);
463 
467  void
468  initialize(MatrixType &m, VectorType &rhs);
469 
477  void
478  initialize(
480 
488  template <class DOFINFO>
489  void
490  initialize_info(DOFINFO &info, bool face) const;
491 
495  template <class DOFINFO>
496  void
497  assemble(const DOFINFO &info);
498 
503  template <class DOFINFO>
504  void
505  assemble(const DOFINFO &info1, const DOFINFO &info2);
506 
507  private:
512  void
513  assemble(const FullMatrix<double> & M,
514  const Vector<double> & vector,
515  const unsigned int index,
516  const std::vector<types::global_dof_index> &indices);
517 
518  void
519  assemble(const FullMatrix<double> & M,
520  const Vector<double> & vector,
521  const unsigned int index,
522  const std::vector<types::global_dof_index> &i1,
523  const std::vector<types::global_dof_index> &i2);
524  };
525 
526 
527  //----------------------------------------------------------------------//
528 
529  template <typename VectorType>
530  inline void
532  {
533  residuals = results;
534  }
535 
536  template <typename VectorType>
537  inline void
540  {
541  constraints = &c;
542  }
543 
544 
545  template <typename MatrixType>
546  inline void
548  {}
549 
550 
551  template <typename VectorType>
552  template <class DOFINFO>
553  inline void
555  {
556  info.initialize_vectors(residuals.size());
557  }
558 
559 
560  template <typename VectorType>
561  template <class DOFINFO>
562  inline void
564  {
565  for (unsigned int k = 0; k < residuals.size(); ++k)
566  {
567  VectorType *v = residuals.entry<VectorType *>(k);
568  for (unsigned int i = 0; i != info.vector(k).n_blocks(); ++i)
569  {
570  const std::vector<types::global_dof_index> &ldi =
571  info.vector(k).n_blocks() == 1 ? info.indices :
572  info.indices_by_block[i];
573 
574  if (constraints != nullptr)
575  constraints->distribute_local_to_global(info.vector(k).block(i),
576  ldi,
577  *v);
578  else
579  v->add(ldi, info.vector(k).block(i));
580  }
581  }
582  }
583 
584  template <typename VectorType>
585  template <class DOFINFO>
586  inline void
588  const DOFINFO &info2)
589  {
590  assemble(info1);
591  assemble(info2);
592  }
593 
594 
595  //----------------------------------------------------------------------//
596 
597  template <typename MatrixType>
598  inline MatrixSimple<MatrixType>::MatrixSimple(double threshold)
599  : threshold(threshold)
600  {}
601 
602 
603  template <typename MatrixType>
604  inline void
606  {
607  matrix.resize(1);
608  matrix[0] = &m;
609  }
610 
611 
612  template <typename MatrixType>
613  inline void
614  MatrixSimple<MatrixType>::initialize(std::vector<MatrixType> &m)
615  {
616  matrix.resize(m.size());
617  for (unsigned int i = 0; i < m.size(); ++i)
618  matrix[i] = &m[i];
619  }
620 
621 
622  template <typename MatrixType>
623  inline void
626  {
627  constraints = &c;
628  }
629 
630 
631  template <typename MatrixType>
632  template <class DOFINFO>
633  inline void
634  MatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
635  {
636  Assert(matrix.size() != 0, ExcNotInitialized());
637 
638  const unsigned int n = info.indices_by_block.size();
639 
640  if (n == 0)
641  info.initialize_matrices(matrix.size(), face);
642  else
643  {
644  info.initialize_matrices(matrix.size() * n * n, face);
645  unsigned int k = 0;
646  for (unsigned int m = 0; m < matrix.size(); ++m)
647  for (unsigned int i = 0; i < n; ++i)
648  for (unsigned int j = 0; j < n; ++j, ++k)
649  {
650  info.matrix(k, false).row = i;
651  info.matrix(k, false).column = j;
652  if (face)
653  {
654  info.matrix(k, true).row = i;
655  info.matrix(k, true).column = j;
656  }
657  }
658  }
659  }
660 
661 
662 
663  template <typename MatrixType>
664  inline void
666  const FullMatrix<double> & M,
667  const unsigned int index,
668  const std::vector<types::global_dof_index> &i1,
669  const std::vector<types::global_dof_index> &i2)
670  {
671  AssertDimension(M.m(), i1.size());
672  AssertDimension(M.n(), i2.size());
673 
674  if (constraints == nullptr)
675  {
676  for (unsigned int j = 0; j < i1.size(); ++j)
677  for (unsigned int k = 0; k < i2.size(); ++k)
678  if (std::fabs(M(j, k)) >= threshold)
679  matrix[index]->add(i1[j], i2[k], M(j, k));
680  }
681  else
682  constraints->distribute_local_to_global(M, i1, i2, *matrix[index]);
683  }
684 
685 
686  template <typename MatrixType>
687  template <class DOFINFO>
688  inline void
690  {
691  Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
692  const unsigned int n = info.indices_by_block.size();
693 
694  if (n == 0)
695  for (unsigned int m = 0; m < matrix.size(); ++m)
696  assemble(info.matrix(m, false).matrix, m, info.indices, info.indices);
697  else
698  {
699  for (unsigned int m = 0; m < matrix.size(); ++m)
700  for (unsigned int k = 0; k < n * n; ++k)
701  {
702  assemble(
703  info.matrix(k + m * n * n, false).matrix,
704  m,
705  info.indices_by_block[info.matrix(k + m * n * n, false).row],
706  info.indices_by_block[info.matrix(k + m * n * n, false)
707  .column]);
708  }
709  }
710  }
711 
712 
713  template <typename MatrixType>
714  template <class DOFINFO>
715  inline void
717  const DOFINFO &info2)
718  {
719  Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
720  Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
721  AssertDimension(info1.indices_by_block.size(),
722  info2.indices_by_block.size());
723 
724  const unsigned int n = info1.indices_by_block.size();
725 
726  if (n == 0)
727  {
728  for (unsigned int m = 0; m < matrix.size(); ++m)
729  {
730  assemble(info1.matrix(m, false).matrix,
731  m,
732  info1.indices,
733  info1.indices);
734  assemble(info1.matrix(m, true).matrix,
735  m,
736  info1.indices,
737  info2.indices);
738  assemble(info2.matrix(m, false).matrix,
739  m,
740  info2.indices,
741  info2.indices);
742  assemble(info2.matrix(m, true).matrix,
743  m,
744  info2.indices,
745  info1.indices);
746  }
747  }
748  else
749  {
750  for (unsigned int m = 0; m < matrix.size(); ++m)
751  for (unsigned int k = 0; k < n * n; ++k)
752  {
753  const unsigned int row = info1.matrix(k + m * n * n, false).row;
754  const unsigned int column =
755  info1.matrix(k + m * n * n, false).column;
756 
757  assemble(info1.matrix(k + m * n * n, false).matrix,
758  m,
759  info1.indices_by_block[row],
760  info1.indices_by_block[column]);
761  assemble(info1.matrix(k + m * n * n, true).matrix,
762  m,
763  info1.indices_by_block[row],
764  info2.indices_by_block[column]);
765  assemble(info2.matrix(k + m * n * n, false).matrix,
766  m,
767  info2.indices_by_block[row],
768  info2.indices_by_block[column]);
769  assemble(info2.matrix(k + m * n * n, true).matrix,
770  m,
771  info2.indices_by_block[row],
772  info1.indices_by_block[column]);
773  }
774  }
775  }
776 
777 
778  //----------------------------------------------------------------------//
779 
780  template <typename MatrixType>
782  : threshold(threshold)
783  {}
784 
785 
786  template <typename MatrixType>
787  inline void
789  {
790  matrix = &m;
791  }
792 
793  template <typename MatrixType>
794  inline void
796  {
797  mg_constrained_dofs = &c;
798  }
799 
800 
801  template <typename MatrixType>
802  inline void
806  {
807  flux_up = &up;
808  flux_down = &down;
809  }
810 
811 
812  template <typename MatrixType>
813  inline void
817  {
818  interface_in = &in;
819  interface_out = &out;
820  }
821 
822 
823  template <typename MatrixType>
824  template <class DOFINFO>
825  inline void
826  MGMatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
827  {
828  const unsigned int n = info.indices_by_block.size();
829 
830  if (n == 0)
831  info.initialize_matrices(1, face);
832  else
833  {
834  info.initialize_matrices(n * n, face);
835  unsigned int k = 0;
836  for (unsigned int i = 0; i < n; ++i)
837  for (unsigned int j = 0; j < n; ++j, ++k)
838  {
839  info.matrix(k, false).row = i;
840  info.matrix(k, false).column = j;
841  if (face)
842  {
843  info.matrix(k, true).row = i;
844  info.matrix(k, true).column = j;
845  }
846  }
847  }
848  }
849 
850 
851  template <typename MatrixType>
852  inline void
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  {
859  AssertDimension(M.m(), i1.size());
860  AssertDimension(M.n(), i2.size());
861  Assert(mg_constrained_dofs == 0, ExcInternalError());
862  // TODO: Possibly remove this function all together
863 
864  for (unsigned int j = 0; j < i1.size(); ++j)
865  for (unsigned int k = 0; k < i2.size(); ++k)
866  if (std::fabs(M(j, k)) >= threshold)
867  G.add(i1[j], i2[k], M(j, k));
868  }
869 
870 
871  template <typename MatrixType>
872  inline void
874  MatrixType & G,
875  const FullMatrix<double> & M,
876  const std::vector<types::global_dof_index> &i1,
877  const std::vector<types::global_dof_index> &i2,
878  const unsigned int level)
879  {
880  AssertDimension(M.m(), i1.size());
881  AssertDimension(M.n(), i2.size());
882 
883  if (mg_constrained_dofs == nullptr)
884  {
885  for (unsigned int j = 0; j < i1.size(); ++j)
886  for (unsigned int k = 0; k < i2.size(); ++k)
887  if (std::fabs(M(j, k)) >= threshold)
888  G.add(i1[j], i2[k], M(j, k));
889  }
890  else
891  {
892  for (unsigned int j = 0; j < i1.size(); ++j)
893  for (unsigned int k = 0; k < i2.size(); ++k)
894  {
895  // Only enter the local values into the global matrix,
896  // if the value is larger than the threshold
897  if (std::fabs(M(j, k)) < threshold)
898  continue;
899 
900  // Do not enter, if either the row or the column
901  // corresponds to an index on the refinement edge. The
902  // level problems are solved with homogeneous
903  // Dirichlet boundary conditions, therefore we
904  // eliminate these rows and columns. The corresponding
905  // matrix entries are entered by assemble_in() and
906  // assemble_out().
907  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
908  mg_constrained_dofs->at_refinement_edge(level, i2[k]))
909  continue;
910 
911  // At the boundary, only enter the term on the
912  // diagonal, but not the coupling terms
913  if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
914  mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
915  (i1[j] != i2[k]))
916  continue;
917 
918  G.add(i1[j], i2[k], M(j, k));
919  }
920  }
921  }
922 
923 
924  template <typename MatrixType>
925  inline void
927  MatrixType & G,
928  const FullMatrix<double> & M,
929  const std::vector<types::global_dof_index> &i1,
930  const std::vector<types::global_dof_index> &i2,
931  const unsigned int level)
932  {
933  AssertDimension(M.n(), i1.size());
934  AssertDimension(M.m(), i2.size());
935 
936  if (mg_constrained_dofs == nullptr)
937  {
938  for (unsigned int j = 0; j < i1.size(); ++j)
939  for (unsigned int k = 0; k < i2.size(); ++k)
940  if (std::fabs(M(k, j)) >= threshold)
941  G.add(i1[j], i2[k], M(k, j));
942  }
943  else
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(k, j)) >= threshold)
948  if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
949  G.add(i1[j], i2[k], M(k, j));
950  }
951  }
952 
953  template <typename MatrixType>
954  inline void
956  MatrixType & G,
957  const FullMatrix<double> & M,
958  const std::vector<types::global_dof_index> &i1,
959  const std::vector<types::global_dof_index> &i2,
960  const unsigned int level)
961  {
962  AssertDimension(M.m(), i1.size());
963  AssertDimension(M.n(), i2.size());
964 
965  if (mg_constrained_dofs == nullptr)
966  {
967  for (unsigned int j = 0; j < i1.size(); ++j)
968  for (unsigned int k = 0; k < i2.size(); ++k)
969  if (std::fabs(M(j, k)) >= threshold)
970  G.add(i1[j], i2[k], M(j, k));
971  }
972  else
973  {
974  for (unsigned int j = 0; j < i1.size(); ++j)
975  for (unsigned int k = 0; k < i2.size(); ++k)
976  if (std::fabs(M(j, k)) >= threshold)
977  if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
978  G.add(i1[j], i2[k], M(j, k));
979  }
980  }
981 
982  template <typename MatrixType>
983  inline void
985  MatrixType & G,
986  const FullMatrix<double> & M,
987  const std::vector<types::global_dof_index> &i1,
988  const std::vector<types::global_dof_index> &i2,
989  const unsigned int level)
990  {
991  AssertDimension(M.m(), i1.size());
992  AssertDimension(M.n(), i2.size());
993  Assert(mg_constrained_dofs != nullptr, ExcInternalError());
994 
995  for (unsigned int j = 0; j < i1.size(); ++j)
996  for (unsigned int k = 0; k < i2.size(); ++k)
997  if (std::fabs(M(j, k)) >= threshold)
998  // Enter values into matrix only if j corresponds to a
999  // degree of freedom on the refinement edge, k does
1000  // not, and both are not on the boundary. This is part
1001  // the difference between the complete matrix with no
1002  // boundary condition at the refinement edge and
1003  // the matrix assembled above by assemble().
1004 
1005  // Thus the logic is: enter the row if it is
1006  // constrained by hanging node constraints (actually,
1007  // the whole refinement edge), but not if it is
1008  // constrained by a boundary constraint.
1009  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1010  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1011  {
1012  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1013  !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
1014  (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1015  mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
1016  i1[j] == i2[k]))
1017  G.add(i1[j], i2[k], M(j, k));
1018  }
1019  }
1020 
1021 
1022  template <typename MatrixType>
1023  inline void
1025  MatrixType & G,
1026  const FullMatrix<double> & M,
1027  const std::vector<types::global_dof_index> &i1,
1028  const std::vector<types::global_dof_index> &i2,
1029  const unsigned int level)
1030  {
1031  AssertDimension(M.n(), i1.size());
1032  AssertDimension(M.m(), i2.size());
1033  Assert(mg_constrained_dofs != nullptr, ExcInternalError());
1034 
1035  for (unsigned int j = 0; j < i1.size(); ++j)
1036  for (unsigned int k = 0; k < i2.size(); ++k)
1037  if (std::fabs(M(k, j)) >= threshold)
1038  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1039  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1040  {
1041  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1042  !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
1043  (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1044  mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
1045  i1[j] == i2[k]))
1046  G.add(i1[j], i2[k], M(k, j));
1047  }
1048  }
1049 
1050 
1051  template <typename MatrixType>
1052  template <class DOFINFO>
1053  inline void
1055  {
1056  Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
1057  const unsigned int level = info.cell->level();
1058 
1059  if (info.indices_by_block.size() == 0)
1060  {
1061  assemble((*matrix)[level],
1062  info.matrix(0, false).matrix,
1063  info.indices,
1064  info.indices,
1065  level);
1066  if (mg_constrained_dofs != nullptr)
1067  {
1068  assemble_in((*interface_in)[level],
1069  info.matrix(0, false).matrix,
1070  info.indices,
1071  info.indices,
1072  level);
1073  assemble_out((*interface_out)[level],
1074  info.matrix(0, false).matrix,
1075  info.indices,
1076  info.indices,
1077  level);
1078  }
1079  }
1080  else
1081  for (unsigned int k = 0; k < info.n_matrices(); ++k)
1082  {
1083  const unsigned int row = info.matrix(k, false).row;
1084  const unsigned int column = info.matrix(k, false).column;
1085 
1086  assemble((*matrix)[level],
1087  info.matrix(k, false).matrix,
1088  info.indices_by_block[row],
1089  info.indices_by_block[column],
1090  level);
1091 
1092  if (mg_constrained_dofs != nullptr)
1093  {
1094  assemble_in((*interface_in)[level],
1095  info.matrix(k, false).matrix,
1096  info.indices_by_block[row],
1097  info.indices_by_block[column],
1098  level);
1099  assemble_out((*interface_out)[level],
1100  info.matrix(k, false).matrix,
1101  info.indices_by_block[column],
1102  info.indices_by_block[row],
1103  level);
1104  }
1105  }
1106  }
1107 
1108 
1109  template <typename MatrixType>
1110  template <class DOFINFO>
1111  inline void
1113  const DOFINFO &info2)
1114  {
1115  Assert(info1.level_cell, ExcMessage("Cell must access level dofs"));
1116  Assert(info2.level_cell, ExcMessage("Cell must access level dofs"));
1117  const unsigned int level1 = info1.cell->level();
1118  const unsigned int level2 = info2.cell->level();
1119 
1120  if (info1.indices_by_block.size() == 0)
1121  {
1122  if (level1 == level2)
1123  {
1124  assemble((*matrix)[level1],
1125  info1.matrix(0, false).matrix,
1126  info1.indices,
1127  info1.indices,
1128  level1);
1129  assemble((*matrix)[level1],
1130  info1.matrix(0, true).matrix,
1131  info1.indices,
1132  info2.indices,
1133  level1);
1134  assemble((*matrix)[level1],
1135  info2.matrix(0, false).matrix,
1136  info2.indices,
1137  info2.indices,
1138  level1);
1139  assemble((*matrix)[level1],
1140  info2.matrix(0, true).matrix,
1141  info2.indices,
1142  info1.indices,
1143  level1);
1144  }
1145  else
1146  {
1147  Assert(level1 > level2, ExcInternalError());
1148  // Do not add info2.M1,
1149  // which is done by
1150  // the coarser cell
1151  assemble((*matrix)[level1],
1152  info1.matrix(0, false).matrix,
1153  info1.indices,
1154  info1.indices,
1155  level1);
1156  if (level1 > 0)
1157  {
1158  assemble_up((*flux_up)[level1],
1159  info1.matrix(0, true).matrix,
1160  info2.indices,
1161  info1.indices,
1162  level1);
1163  assemble_down((*flux_down)[level1],
1164  info2.matrix(0, true).matrix,
1165  info2.indices,
1166  info1.indices,
1167  level1);
1168  }
1169  }
1170  }
1171  else
1172  for (unsigned int k = 0; k < info1.n_matrices(); ++k)
1173  {
1174  const unsigned int row = info1.matrix(k, false).row;
1175  const unsigned int column = info1.matrix(k, false).column;
1176 
1177  if (level1 == level2)
1178  {
1179  assemble((*matrix)[level1],
1180  info1.matrix(k, false).matrix,
1181  info1.indices_by_block[row],
1182  info1.indices_by_block[column],
1183  level1);
1184  assemble((*matrix)[level1],
1185  info1.matrix(k, true).matrix,
1186  info1.indices_by_block[row],
1187  info2.indices_by_block[column],
1188  level1);
1189  assemble((*matrix)[level1],
1190  info2.matrix(k, false).matrix,
1191  info2.indices_by_block[row],
1192  info2.indices_by_block[column],
1193  level1);
1194  assemble((*matrix)[level1],
1195  info2.matrix(k, true).matrix,
1196  info2.indices_by_block[row],
1197  info1.indices_by_block[column],
1198  level1);
1199  }
1200  else
1201  {
1202  Assert(level1 > level2, ExcInternalError());
1203  // Do not add info2.M1,
1204  // which is done by
1205  // the coarser cell
1206  assemble((*matrix)[level1],
1207  info1.matrix(k, false).matrix,
1208  info1.indices_by_block[row],
1209  info1.indices_by_block[column],
1210  level1);
1211  if (level1 > 0)
1212  {
1213  assemble_up((*flux_up)[level1],
1214  info1.matrix(k, true).matrix,
1215  info2.indices_by_block[column],
1216  info1.indices_by_block[row],
1217  level1);
1218  assemble_down((*flux_down)[level1],
1219  info2.matrix(k, true).matrix,
1220  info2.indices_by_block[row],
1221  info1.indices_by_block[column],
1222  level1);
1223  }
1224  }
1225  }
1226  }
1227 
1228  //----------------------------------------------------------------------//
1229 
1230  template <typename MatrixType, typename VectorType>
1232  : MatrixSimple<MatrixType>(t)
1233  {}
1234 
1235 
1236  template <typename MatrixType, typename VectorType>
1237  inline void
1239  VectorType &rhs)
1240  {
1241  AnyData data;
1242  VectorType *p = &rhs;
1243  data.add(p, "right hand side");
1244 
1247  }
1248 
1249  template <typename MatrixType, typename VectorType>
1250  inline void
1253  {
1255  }
1256 
1257 
1258  template <typename MatrixType, typename VectorType>
1259  template <class DOFINFO>
1260  inline void
1262  bool face) const
1263  {
1266  }
1267 
1268  template <typename MatrixType, typename VectorType>
1269  inline void
1271  const FullMatrix<double> & M,
1272  const Vector<double> & vector,
1273  const unsigned int index,
1274  const std::vector<types::global_dof_index> &indices)
1275  {
1276  AssertDimension(M.m(), indices.size());
1277  AssertDimension(M.n(), indices.size());
1278 
1280  VectorType *v = residuals.entry<VectorType *>(index);
1281 
1283  {
1284  for (unsigned int i = 0; i < indices.size(); ++i)
1285  (*v)(indices[i]) += vector(i);
1286 
1287  for (unsigned int j = 0; j < indices.size(); ++j)
1288  for (unsigned int k = 0; k < indices.size(); ++k)
1289  if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
1290  MatrixSimple<MatrixType>::matrix[index]->add(indices[j],
1291  indices[k],
1292  M(j, k));
1293  }
1294  else
1295  {
1296  ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1297  M,
1298  vector,
1299  indices,
1301  *v,
1302  true);
1303  }
1304  }
1305 
1306  template <typename MatrixType, typename VectorType>
1307  inline void
1309  const FullMatrix<double> & M,
1310  const Vector<double> & vector,
1311  const unsigned int index,
1312  const std::vector<types::global_dof_index> &i1,
1313  const std::vector<types::global_dof_index> &i2)
1314  {
1315  AssertDimension(M.m(), i1.size());
1316  AssertDimension(M.n(), i2.size());
1317 
1319  VectorType *v = residuals.entry<VectorType *>(index);
1320 
1322  {
1323  for (unsigned int j = 0; j < i1.size(); ++j)
1324  for (unsigned int k = 0; k < i2.size(); ++k)
1325  if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
1326  MatrixSimple<MatrixType>::matrix[index]->add(i1[j],
1327  i2[k],
1328  M(j, k));
1329  }
1330  else
1331  {
1332  ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1333  vector, i1, i2, *v, M, false);
1334  ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1335  M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
1336  }
1337  }
1338 
1339 
1340  template <typename MatrixType, typename VectorType>
1341  template <class DOFINFO>
1342  inline void
1344  {
1347  Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
1348  const unsigned int n = info.indices_by_block.size();
1349 
1350  if (n == 0)
1351  {
1352  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1353  ++m)
1354  assemble(info.matrix(m, false).matrix,
1355  info.vector(m).block(0),
1356  m,
1357  info.indices);
1358  }
1359  else
1360  {
1361  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1362  ++m)
1363  for (unsigned int k = 0; k < n * n; ++k)
1364  {
1365  const unsigned int row = info.matrix(k + m * n * n, false).row;
1366  const unsigned int column =
1367  info.matrix(k + m * n * n, false).column;
1368 
1369  if (row == column)
1370  assemble(info.matrix(k + m * n * n, false).matrix,
1371  info.vector(m).block(row),
1372  m,
1373  info.indices_by_block[row]);
1374  else
1375  assemble(info.matrix(k + m * n * n, false).matrix,
1376  info.vector(m).block(row),
1377  m,
1378  info.indices_by_block[row],
1379  info.indices_by_block[column]);
1380  }
1381  }
1382  }
1383 
1384 
1385  template <typename MatrixType, typename VectorType>
1386  template <class DOFINFO>
1387  inline void
1389  const DOFINFO &info2)
1390  {
1391  Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
1392  Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
1393  AssertDimension(info1.indices_by_block.size(),
1394  info2.indices_by_block.size());
1395 
1396  const unsigned int n = info1.indices_by_block.size();
1397 
1398  if (n == 0)
1399  {
1400  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1401  ++m)
1402  {
1403  assemble(info1.matrix(m, false).matrix,
1404  info1.vector(m).block(0),
1405  m,
1406  info1.indices);
1407  assemble(info1.matrix(m, true).matrix,
1408  info1.vector(m).block(0),
1409  m,
1410  info1.indices,
1411  info2.indices);
1412  assemble(info2.matrix(m, false).matrix,
1413  info2.vector(m).block(0),
1414  m,
1415  info2.indices);
1416  assemble(info2.matrix(m, true).matrix,
1417  info2.vector(m).block(0),
1418  m,
1419  info2.indices,
1420  info1.indices);
1421  }
1422  }
1423  else
1424  {
1425  for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1426  ++m)
1427  for (unsigned int k = 0; k < n * n; ++k)
1428  {
1429  const unsigned int row = info1.matrix(k + m * n * n, false).row;
1430  const unsigned int column =
1431  info1.matrix(k + m * n * n, false).column;
1432 
1433  if (row == column)
1434  {
1435  assemble(info1.matrix(k + m * n * n, false).matrix,
1436  info1.vector(m).block(row),
1437  m,
1438  info1.indices_by_block[row]);
1439  assemble(info2.matrix(k + m * n * n, false).matrix,
1440  info2.vector(m).block(row),
1441  m,
1442  info2.indices_by_block[row]);
1443  }
1444  else
1445  {
1446  assemble(info1.matrix(k + m * n * n, false).matrix,
1447  info1.vector(m).block(row),
1448  m,
1449  info1.indices_by_block[row],
1450  info1.indices_by_block[column]);
1451  assemble(info2.matrix(k + m * n * n, false).matrix,
1452  info2.vector(m).block(row),
1453  m,
1454  info2.indices_by_block[row],
1455  info2.indices_by_block[column]);
1456  }
1457  assemble(info1.matrix(k + m * n * n, true).matrix,
1458  info1.vector(m).block(row),
1459  m,
1460  info1.indices_by_block[row],
1461  info2.indices_by_block[column]);
1462  assemble(info2.matrix(k + m * n * n, true).matrix,
1463  info2.vector(m).block(row),
1464  m,
1465  info2.indices_by_block[row],
1466  info1.indices_by_block[column]);
1467  }
1468  }
1469  }
1470  } // namespace Assembler
1471 } // namespace MeshWorker
1472 
1473 DEAL_II_NAMESPACE_CLOSE
1474 
1475 #endif
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > flux_down
Definition: simple.h:414
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void initialize_interfaces(MGLevelObject< MatrixType > &interface_in, MGLevelObject< MatrixType > &interface_out)
Definition: simple.h:814
SystemSimple(double threshold=1.e-12)
Definition: simple.h:1231
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::vector< SmartPointer< MatrixType, MatrixSimple< MatrixType > > > matrix
Definition: simple.h:231
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:1024
void initialize(MGLevelObject< MatrixType > &m)
Definition: simple.h:788
static ::ExceptionBase & ExcNotInitialized()
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > matrix
Definition: simple.h:400
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > interface_out
Definition: simple.h:428
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:926
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:554
void initialize(MatrixType &m, VectorType &rhs)
Definition: simple.h:1238
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:984
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
void assemble(const DOFINFO &info)
Definition: simple.h:563
void initialize_fluxes(MGLevelObject< MatrixType > &flux_up, MGLevelObject< MatrixType > &flux_down)
Definition: simple.h:803
void initialize(AnyData &results)
Definition: simple.h:531
MGMatrixSimple(double threshold=1.e-12)
Definition: simple.h:781
SmartPointer< const AffineConstraints< typename VectorType::value_type >, ResidualSimple< VectorType > > constraints
Definition: simple.h:132
void initialize_local_blocks(const BlockIndices &)
Definition: simple.h:547
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:955
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > interface_in
Definition: simple.h:421
void add(type entry, const std::string &name)
Add a new data object.
Definition: any_data.h:432
SmartPointer< MGLevelObject< MatrixType >, MGMatrixSimple< MatrixType > > flux_up
Definition: simple.h:407
void assemble(const DOFINFO &info)
Definition: simple.h:1343
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:826
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:1261
void assemble(const DOFINFO &info)
Definition: simple.h:1054
void initialize(MatrixType &m)
Definition: simple.h:605
SmartPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixSimple< MatrixType > > constraints
Definition: simple.h:256
MatrixSimple(double threshold=1.e-12)
Definition: simple.h:598
void assemble(const DOFINFO &info)
Definition: simple.h:689
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:634
static ::ExceptionBase & ExcInternalError()
SmartPointer< const MGConstrainedDoFs, MGMatrixSimple< MatrixType > > mg_constrained_dofs
Definition: simple.h:433