Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
simple.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 // @f$Id: simple.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #ifndef __deal2__mesh_worker_simple_h
19 #define __deal2__mesh_worker_simple_h
20 
21 #include <deal.II/base/named_data.h>
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/base/mg_level_object.h>
24 #include <deal.II/lac/block_vector.h>
25 #include <deal.II/meshworker/dof_info.h>
26 #include <deal.II/meshworker/functional.h>
27 #include <deal.II/multigrid/mg_constrained_dofs.h>
28 
41 
42 namespace MeshWorker
43 {
44  namespace Assembler
45  {
59  template <class VECTOR>
61  {
62  public:
63  void initialize(NamedData<VECTOR *> &results);
64 
68  void initialize(const ConstraintMatrix &constraints);
69 
86 
94  template <class DOFINFO>
95  void initialize_info(DOFINFO &info, bool face) const;
96 
104  template <class DOFINFO>
105  void assemble(const DOFINFO &info);
106 
110  template<class DOFINFO>
111  void assemble(const DOFINFO &info1,
112  const DOFINFO &info2);
113  private:
122  };
123 
124 
151  template <class MATRIX>
153  {
154  public:
159  MatrixSimple(double threshold = 1.e-12);
160 
164  void initialize(MATRIX &m);
173 
190 
198  template <class DOFINFO>
199  void initialize_info(DOFINFO &info, bool face) const;
200 
204  template<class DOFINFO>
205  void assemble(const DOFINFO &info);
206 
211  template<class DOFINFO>
212  void assemble(const DOFINFO &info1,
213  const DOFINFO &info2);
214  private:
219  void assemble(const FullMatrix<double> &M,
220  const std::vector<types::global_dof_index> &i1,
221  const std::vector<types::global_dof_index> &i2);
222 
233 
239  const double threshold;
240 
241  };
242 
243 
254  template <class MATRIX>
256  {
257  public:
262  MGMatrixSimple(double threshold = 1.e-12);
263 
268 
273 
289 
296 
311  template <class DOFINFO>
312  void initialize_info(DOFINFO &info, bool face) const;
313 
317  template<class DOFINFO>
318  void assemble(const DOFINFO &info);
319 
324  template<class DOFINFO>
325  void assemble(const DOFINFO &info1,
326  const DOFINFO &info2);
327  private:
331  void assemble(MATRIX &G,
332  const FullMatrix<double> &M,
333  const std::vector<types::global_dof_index> &i1,
334  const std::vector<types::global_dof_index> &i2);
335 
339  void assemble(MATRIX &G,
340  const FullMatrix<double> &M,
341  const std::vector<types::global_dof_index> &i1,
342  const std::vector<types::global_dof_index> &i2,
343  const unsigned int level);
344 
349  void assemble_up(MATRIX &G,
350  const FullMatrix<double> &M,
351  const std::vector<types::global_dof_index> &i1,
352  const std::vector<types::global_dof_index> &i2,
353  const unsigned int level = numbers::invalid_unsigned_int);
358  void assemble_down(MATRIX &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);
363 
368  void assemble_in(MATRIX &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 assemble_out(MATRIX &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 
394 
400 
406 
416 
422  const double threshold;
423 
424  };
425 
426 
437  template <class MATRIX, class VECTOR>
438  class SystemSimple :
439  private MatrixSimple<MATRIX>,
440  private ResidualSimple<VECTOR>
441  {
442  public:
448  SystemSimple(double threshold = 1.e-12);
449 
454  void initialize(MATRIX &m, VECTOR &rhs);
455 
464 
478  template <class DOFINFO>
479  void initialize_info(DOFINFO &info, bool face) const;
480 
486  template<class DOFINFO>
487  void assemble(const DOFINFO &info);
488 
495  template<class DOFINFO>
496  void assemble(const DOFINFO &info1,
497  const DOFINFO &info2);
498  };
499 
500 
501 //----------------------------------------------------------------------//
502 
503  template <class VECTOR>
504  inline void
506  {
507  residuals = results;
508  }
509 
510  template <class VECTOR>
511  inline void
513  {
514  constraints = &c;
515  }
516 
517 
518  template <class MATRIX>
519  inline void
521  {}
522 
523 
524  template <class VECTOR>
525  template <class DOFINFO>
526  inline void
527  ResidualSimple<VECTOR>::initialize_info(DOFINFO &info, bool) const
528  {
529  info.initialize_vectors(residuals.size());
530  }
531 
532 
533  template <class VECTOR>
534  template <class DOFINFO>
535  inline void
536  ResidualSimple<VECTOR>::assemble(const DOFINFO &info)
537  {
538  for (unsigned int k=0; k<residuals.size(); ++k)
539  {
540  if (constraints == 0)
541  {
542  for (unsigned int i=0; i<info.vector(k).block(0).size(); ++i)
543  (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
544  }
545  else
546  {
547  if (info.indices_by_block.size() == 0)
548  constraints->distribute_local_to_global(info.vector(k).block(0), info.indices, (*residuals(k)));
549  else
550  for (unsigned int i=0; i != info.vector(k).n_blocks(); ++i)
551  constraints->distribute_local_to_global(info.vector(k).block(i), info.indices_by_block[i], (*residuals(k)));
552  }
553  }
554  }
555 
556  template <class VECTOR>
557  template <class DOFINFO>
558  inline void
559  ResidualSimple<VECTOR>::assemble(const DOFINFO &info1,
560  const DOFINFO &info2)
561  {
562  for (unsigned int k=0; k<residuals.size(); ++k)
563  {
564  if (constraints == 0)
565  {
566  for (unsigned int i=0; i<info1.vector(k).block(0).size(); ++i)
567  (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
568  for (unsigned int i=0; i<info2.vector(k).block(0).size(); ++i)
569  (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
570  }
571  else
572  {
573  if (info1.indices_by_block.size() == 0 && info2.indices_by_block.size() == 0)
574  {
575  constraints->distribute_local_to_global
576  (info1.vector(k).block(0), info1.indices, (*residuals(k)));
577  constraints->distribute_local_to_global
578  (info2.vector(k).block(0), info2.indices, (*residuals(k)));
579  }
580  else if (info1.indices_by_block.size() != 0 && info2.indices_by_block.size() != 0)
581  {
582  for (unsigned int i=0; i<info1.vector(k).n_blocks(); ++i)
583  {
584  constraints->distribute_local_to_global
585  (info1.vector(k).block(i), info1.indices_by_block[i], (*residuals(k)));
586  constraints->distribute_local_to_global
587  (info2.vector(k).block(i), info2.indices_by_block[i], (*residuals(k)));
588  }
589  }
590  else
591  {
592  Assert(false, ExcNotImplemented());
593  }
594  }
595  }
596  }
597 
598 
599 //----------------------------------------------------------------------//
600 
601  template <class MATRIX>
602  inline
604  :
605  threshold(threshold)
606  {}
607 
608 
609  template <class MATRIX>
610  inline void
612  {
613  matrix = &m;
614  }
615 
616 
617  template <class MATRIX>
618  inline void
620  {
621  constraints = &c;
622  }
623 
624 
625  template <class MATRIX>
626  inline void
628  {}
629 
630 
631  template <class MATRIX >
632  template <class DOFINFO>
633  inline void
634  MatrixSimple<MATRIX>::initialize_info(DOFINFO &info, bool face) const
635  {
636  const unsigned int n = info.indices_by_block.size();
637 
638  if (n == 0)
639  info.initialize_matrices(1, face);
640  else
641  {
642  info.initialize_matrices(n*n, face);
643  unsigned int k=0;
644  for (unsigned int i=0; i<n; ++i)
645  for (unsigned int j=0; j<n; ++j,++k)
646  {
647  info.matrix(k,false).row = i;
648  info.matrix(k,false).column = j;
649  if (face)
650  {
651  info.matrix(k,true).row = i;
652  info.matrix(k,true).column = j;
653  }
654  }
655  }
656  }
657 
658 
659 
660  template <class MATRIX>
661  inline void
663  const std::vector<types::global_dof_index> &i1,
664  const std::vector<types::global_dof_index> &i2)
665  {
666  AssertDimension(M.m(), i1.size());
667  AssertDimension(M.n(), i2.size());
668 
669  if (constraints == 0)
670  {
671  for (unsigned int j=0; j<i1.size(); ++j)
672  for (unsigned int k=0; k<i2.size(); ++k)
673  if (std::fabs(M(j,k)) >= threshold)
674  matrix->add(i1[j], i2[k], M(j,k));
675  }
676  else
677  constraints->distribute_local_to_global(M, i1, i2, *matrix);
678  }
679 
680 
681  template <class MATRIX>
682  template <class DOFINFO>
683  inline void
684  MatrixSimple<MATRIX>::assemble(const DOFINFO &info)
685  {
686  Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
687 
688  if (info.indices_by_block.size() == 0)
689  assemble(info.matrix(0,false).matrix, info.indices, info.indices);
690  else
691  {
692  for (unsigned int k=0; k<info.n_matrices(); ++k)
693  {
694  assemble(info.matrix(k,false).matrix,
695  info.indices_by_block[info.matrix(k,false).row],
696  info.indices_by_block[info.matrix(k,false).column]);
697  }
698  }
699  }
700 
701 
702  template <class MATRIX>
703  template <class DOFINFO>
704  inline void
705  MatrixSimple<MATRIX>::assemble(const DOFINFO &info1, const DOFINFO &info2)
706  {
707  Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
708  Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
709 
710  if (info1.indices_by_block.size() == 0 && info2.indices_by_block.size() == 0)
711  {
712  assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices);
713  assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices);
714  assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices);
715  assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices);
716  }
717  else if (info1.indices_by_block.size() != 0 && info2.indices_by_block.size() != 0)
718  for (unsigned int k=0; k<info1.n_matrices(); ++k)
719  {
720  const unsigned int row = info1.matrix(k,false).row;
721  const unsigned int column = info1.matrix(k,false).column;
722 
723  assemble(info1.matrix(k,false).matrix,
724  info1.indices_by_block[row], info1.indices_by_block[column]);
725  assemble(info1.matrix(k,true).matrix,
726  info1.indices_by_block[row], info2.indices_by_block[column]);
727  assemble(info2.matrix(k,false).matrix,
728  info2.indices_by_block[row], info2.indices_by_block[column]);
729  assemble(info2.matrix(k,true).matrix,
730  info2.indices_by_block[row], info1.indices_by_block[column]);
731  }
732  else
733  {
734  Assert(false, ExcNotImplemented());
735  }
736  }
737 
738 
739 //----------------------------------------------------------------------//
740 
741  template <class MATRIX>
742  inline
744  :
745  threshold(threshold)
746  {}
747 
748 
749  template <class MATRIX>
750  inline void
752  {
753  matrix = &m;
754  }
755 
756  template <class MATRIX>
757  inline void
759  {
760  mg_constrained_dofs = &c;
761  }
762 
763  template <class MATRIX>
764  inline void
766  {}
767 
768 
769  template <class MATRIX>
770  inline void
773  {
774  flux_up = &up;
775  flux_down = &down;
776  }
777 
778 
779  template <class MATRIX>
780  inline void
783  {
784  interface_in = &in;
785  interface_out = &out;
786  }
787 
788 
789  template <class MATRIX >
790  template <class DOFINFO>
791  inline void
792  MGMatrixSimple<MATRIX>::initialize_info(DOFINFO &info, bool face) const
793  {
794  const unsigned int n = info.indices_by_block.size();
795 
796  if (n == 0)
797  info.initialize_matrices(1, face);
798  else
799  {
800  info.initialize_matrices(n*n, face);
801  unsigned int k=0;
802  for (unsigned int i=0; i<n; ++i)
803  for (unsigned int j=0; j<n; ++j,++k)
804  {
805  info.matrix(k,false).row = i;
806  info.matrix(k,false).column = j;
807  if (face)
808  {
809  info.matrix(k,true).row = i;
810  info.matrix(k,true).column = j;
811  }
812  }
813  }
814  }
815 
816 
817  template <class MATRIX>
818  inline void
820  MATRIX &G,
821  const FullMatrix<double> &M,
822  const std::vector<types::global_dof_index> &i1,
823  const std::vector<types::global_dof_index> &i2)
824  {
825  AssertDimension(M.m(), i1.size());
826  AssertDimension(M.n(), i2.size());
827 
828  if (mg_constrained_dofs == 0)
829  {
830  for (unsigned int j=0; j<i1.size(); ++j)
831  for (unsigned int k=0; k<i2.size(); ++k)
832  if (std::fabs(M(j,k)) >= threshold)
833  G.add(i1[j], i2[k], M(j,k));
834  }
835  else
836  {
837  for (unsigned int j=0; j<i1.size(); ++j)
838  for (unsigned int k=0; k<i2.size(); ++k)
839  if (std::fabs(M(j,k)) >= threshold)
840  {
841  if (!mg_constrained_dofs->continuity_across_refinement_edges())
842  G.add(i1[j], i2[k], M(j,k));
843  }
844  }
845  }
846 
847 
848  template <class MATRIX>
849  inline void
851  MATRIX &G,
852  const FullMatrix<double> &M,
853  const std::vector<types::global_dof_index> &i1,
854  const std::vector<types::global_dof_index> &i2,
855  const unsigned int level)
856  {
857  AssertDimension(M.m(), i1.size());
858  AssertDimension(M.n(), i2.size());
859 
860  if (mg_constrained_dofs == 0)
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  else
868  {
869  for (unsigned int j=0; j<i1.size(); ++j)
870  for (unsigned int k=0; k<i2.size(); ++k)
871  if (std::fabs(M(j,k)) >= threshold)
872  if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
873  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
874  {
875  if (mg_constrained_dofs->set_boundary_values())
876  {
877  // At the boundary, only enter the term on the
878  // diagonal, but not the coupling terms
879  if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
880  !mg_constrained_dofs->is_boundary_index(level, i2[k]))
881  ||
882  (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
883  mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
884  i1[j] == i2[k]))
885  G.add(i1[j], i2[k], M(j,k));
886  }
887  else
888  G.add(i1[j], i2[k], M(j,k));
889  }
890  }
891  }
892 
893 
894  template <class MATRIX>
895  inline void
897  MATRIX &G,
898  const FullMatrix<double> &M,
899  const std::vector<types::global_dof_index> &i1,
900  const std::vector<types::global_dof_index> &i2,
901  const unsigned int level)
902  {
903  AssertDimension(M.n(), i1.size());
904  AssertDimension(M.m(), i2.size());
905 
906  if (mg_constrained_dofs == 0)
907  {
908  for (unsigned int j=0; j<i1.size(); ++j)
909  for (unsigned int k=0; k<i2.size(); ++k)
910  if (std::fabs(M(k,j)) >= threshold)
911  G.add(i1[j], i2[k], M(k,j));
912  }
913  else
914  {
915  for (unsigned int j=0; j<i1.size(); ++j)
916  for (unsigned int k=0; k<i2.size(); ++k)
917  if (std::fabs(M(k,j)) >= threshold)
918  if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
919  G.add(i1[j], i2[k], M(k,j));
920  }
921  }
922 
923  template <class MATRIX>
924  inline void
926  MATRIX &G,
927  const FullMatrix<double> &M,
928  const std::vector<types::global_dof_index> &i1,
929  const std::vector<types::global_dof_index> &i2,
930  const unsigned int level)
931  {
932  AssertDimension(M.m(), i1.size());
933  AssertDimension(M.n(), i2.size());
934 
935  if (mg_constrained_dofs == 0)
936  {
937  for (unsigned int j=0; j<i1.size(); ++j)
938  for (unsigned int k=0; k<i2.size(); ++k)
939  if (std::fabs(M(j,k)) >= threshold)
940  G.add(i1[j], i2[k], M(j,k));
941  }
942  else
943  {
944  for (unsigned int j=0; j<i1.size(); ++j)
945  for (unsigned int k=0; k<i2.size(); ++k)
946  if (std::fabs(M(j,k)) >= threshold)
947  if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
948  G.add(i1[j], i2[k], M(j,k));
949  }
950  }
951 
952  template <class MATRIX>
953  inline void
955  MATRIX &G,
956  const FullMatrix<double> &M,
957  const std::vector<types::global_dof_index> &i1,
958  const std::vector<types::global_dof_index> &i2,
959  const unsigned int level)
960  {
961  AssertDimension(M.m(), i1.size());
962  AssertDimension(M.n(), i2.size());
963  Assert(mg_constrained_dofs != 0, ExcInternalError());
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  // Enter values into matrix only if j corresponds to a
969  // degree of freedom on the refinemenent edge, k does
970  // not, and both are not on the boundary. This is part
971  // the difference between the complete matrix with no
972  // boundary condition at the refinement edge and and
973  // the matrix assembled above by assemble().
974 
975  // Thus the logic is: enter the row if it is
976  // constrained by hanging node constraints (actually,
977  // the whole refinement edge), but not if it is
978  // constrained by a boundary constraint.
979  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
980  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
981  {
982  if (mg_constrained_dofs->set_boundary_values())
983  {
984  if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
985  !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
986  ||
987  (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
988  mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
989  i1[j] == i2[k]))
990  G.add(i1[j], i2[k], M(j,k));
991  }
992  else
993  G.add(i1[j], i2[k], M(j,k));
994  }
995  }
996 
997 
998  template <class MATRIX>
999  inline void
1001  MATRIX &G,
1002  const FullMatrix<double> &M,
1003  const std::vector<types::global_dof_index> &i1,
1004  const std::vector<types::global_dof_index> &i2,
1005  const unsigned int level)
1006  {
1007  AssertDimension(M.n(), i1.size());
1008  AssertDimension(M.m(), i2.size());
1009  Assert(mg_constrained_dofs != 0, ExcInternalError());
1010 
1011  for (unsigned int j=0; j<i1.size(); ++j)
1012  for (unsigned int k=0; k<i2.size(); ++k)
1013  if (std::fabs(M(k,j)) >= threshold)
1014  if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1015  !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1016  {
1017  if (mg_constrained_dofs->set_boundary_values())
1018  {
1019  if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
1020  !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
1021  ||
1022  (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
1023  mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
1024  i1[j] == i2[k]))
1025  G.add(i1[j], i2[k], M(k,j));
1026  }
1027  else
1028  G.add(i1[j], i2[k], M(k,j));
1029  }
1030  }
1031 
1032 
1033  template <class MATRIX>
1034  template <class DOFINFO>
1035  inline void
1037  {
1038  Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
1039  const unsigned int level = info.cell->level();
1040 
1041  if (info.indices_by_block.size() == 0)
1042  {
1043  assemble((*matrix)[level], info.matrix(0,false).matrix,
1044  info.indices, info.indices, level);
1045  if (mg_constrained_dofs != 0)
1046  {
1047  assemble_in((*interface_in)[level], info.matrix(0,false).matrix,
1048  info.indices, info.indices, level);
1049  assemble_out((*interface_out)[level],info.matrix(0,false).matrix,
1050  info.indices, info.indices, level);
1051  }
1052  }
1053  else
1054  for (unsigned int k=0; k<info.n_matrices(); ++k)
1055  {
1056  const unsigned int row = info.matrix(k,false).row;
1057  const unsigned int column = info.matrix(k,false).column;
1058 
1059  assemble((*matrix)[level], info.matrix(k,false).matrix,
1060  info.indices_by_block[row], info.indices_by_block[column], level);
1061 
1062  if (mg_constrained_dofs != 0)
1063  {
1064  assemble_in((*interface_in)[level], info.matrix(k,false).matrix,
1065  info.indices_by_block[row], info.indices_by_block[column], level);
1066  assemble_out((*interface_out)[level],info.matrix(k,false).matrix,
1067  info.indices_by_block[column], info.indices_by_block[row], level);
1068  }
1069  }
1070  }
1071 
1072 
1073  template <class MATRIX>
1074  template <class DOFINFO>
1075  inline void
1076  MGMatrixSimple<MATRIX>::assemble(const DOFINFO &info1,
1077  const DOFINFO &info2)
1078  {
1079  Assert(info1.level_cell, ExcMessage("Cell must access level dofs"));
1080  Assert(info2.level_cell, ExcMessage("Cell must access level dofs"));
1081  const unsigned int level1 = info1.cell->level();
1082  const unsigned int level2 = info2.cell->level();
1083 
1084  if (info1.indices_by_block.size() == 0)
1085  {
1086  if (level1 == level2)
1087  {
1088  if (mg_constrained_dofs == 0)
1089  {
1090  assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
1091  assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices);
1092  assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
1093  assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
1094  }
1095  else
1096  {
1097  assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
1098  assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
1099  assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
1100  assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
1101  }
1102  }
1103  else
1104  {
1105  Assert(level1 > level2, ExcInternalError());
1106  // Do not add info2.M1,
1107  // which is done by
1108  // the coarser cell
1109  assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
1110  if (level1>0)
1111  {
1112  assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
1113  assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
1114  }
1115  }
1116  }
1117  else
1118  for (unsigned int k=0; k<info1.n_matrices(); ++k)
1119  {
1120  const unsigned int row = info1.matrix(k,false).row;
1121  const unsigned int column = info1.matrix(k,false).column;
1122 
1123  if (level1 == level2)
1124  {
1125  if (mg_constrained_dofs == 0)
1126  {
1127  assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
1128  assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column]);
1129  assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column]);
1130  assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column]);
1131  }
1132  else
1133  {
1134  assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
1135  assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
1136  assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
1137  assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
1138  }
1139  }
1140  else
1141  {
1142  Assert(level1 > level2, ExcInternalError());
1143  // Do not add info2.M1,
1144  // which is done by
1145  // the coarser cell
1146  assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
1147  if (level1>0)
1148  {
1149  assemble_up((*flux_up)[level1],info1.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
1150  assemble_down((*flux_down)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
1151  }
1152  }
1153  }
1154  }
1155 
1156 //----------------------------------------------------------------------//
1157 
1158  template <class MATRIX, class VECTOR>
1160  :
1161  MatrixSimple<MATRIX>(t)
1162  {}
1163 
1164 
1165  template <class MATRIX, class VECTOR>
1166  inline void
1168  {
1169  NamedData<VECTOR *> data;
1170  VECTOR *p = &rhs;
1171  data.add(p, "right hand side");
1172 
1175  }
1176 
1177  template <class MATRIX, class VECTOR>
1178  inline void
1180  {
1183  }
1184 
1185 
1186  template <class MATRIX, class VECTOR>
1187  template <class DOFINFO>
1188  inline void
1190  bool face) const
1191  {
1194  }
1195 
1196 
1197  template <class MATRIX, class VECTOR>
1198  template <class DOFINFO>
1199  inline void
1201  {
1204  }
1205 
1206 
1207  template <class MATRIX, class VECTOR>
1208  template <class DOFINFO>
1209  inline void
1211  const DOFINFO &info2)
1212  {
1213  MatrixSimple<MATRIX>::assemble(info1, info2);
1214  ResidualSimple<VECTOR>::assemble(info1, info2);
1215  }
1216  }
1217 }
1218 
1219 DEAL_II_NAMESPACE_CLOSE
1220 
1221 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void assemble(const DOFINFO &info)
Definition: simple.h:1200
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:634
void assemble(const DOFINFO &info)
Definition: simple.h:1036
void add(DATA &v, const std::string &name)
Definition: named_data.h:253
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:1189
size_type m() const
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > flux_up
Definition: simple.h:393
::ExceptionBase & ExcMessage(std::string arg1)
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:55
MatrixSimple(double threshold=1.e-12)
Definition: simple.h:603
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > interface_in
Definition: simple.h:405
void initialize_local_blocks(const BlockIndices &)
Definition: simple.h:627
size_type n() const
void initialize_local_blocks(const BlockIndices &)
Definition: simple.h:765
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > matrix
Definition: simple.h:387
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:527
void initialize_fluxes(MGLevelObject< MATRIX > &flux_up, MGLevelObject< MATRIX > &flux_down)
Definition: simple.h:771
void assemble_down(MATRIX &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:925
MGMatrixSimple(double threshold=1.e-12)
Definition: simple.h:743
SmartPointer< const MGConstrainedDoFs, MGMatrixSimple< MATRIX > > mg_constrained_dofs
Definition: simple.h:415
#define Assert(cond, exc)
Definition: exceptions.h:299
SmartPointer< const ConstraintMatrix, MatrixSimple< MATRIX > > constraints
Definition: simple.h:232
void assemble_out(MATRIX &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:1000
SystemSimple(double threshold=1.e-12)
Definition: simple.h:1159
void initialize_local_blocks(const BlockIndices &)
Definition: simple.h:520
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > flux_down
Definition: simple.h:399
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > interface_out
Definition: simple.h:411
NamedData< SmartPointer< VECTOR, ResidualSimple< VECTOR > > > residuals
Definition: simple.h:117
void initialize(MATRIX &m, VECTOR &rhs)
Definition: simple.h:1167
void assemble(const DOFINFO &info)
Definition: simple.h:536
void initialize_info(DOFINFO &info, bool face) const
Definition: simple.h:792
void assemble_up(MATRIX &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:896
void assemble_in(MATRIX &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:954
SmartPointer< const ConstraintMatrix, ResidualSimple< VECTOR > > constraints
Definition: simple.h:121
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcInternalError()
SmartPointer< MATRIX, MatrixSimple< MATRIX > > matrix
Definition: simple.h:227
void assemble(const DOFINFO &info)
Definition: simple.h:684
void initialize(MGLevelObject< MATRIX > &m)
Definition: simple.h:751
void initialize_interfaces(MGLevelObject< MATRIX > &interface_in, MGLevelObject< MATRIX > &interface_out)
Definition: simple.h:781