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