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