Reference documentation for deal.II version 8.1.0
assembler.h
1 // ---------------------------------------------------------------------
2 // @f$Id: assembler.h 30036 2013-07-18 16:55:32Z maier @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_assembler_h
19 #define __deal2__mesh_worker_assembler_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>
28 #include <deal.II/multigrid/mg_constrained_dofs.h>
29 
30 
32 
33 namespace MeshWorker
34 {
86  namespace Assembler
87  {
109  template <class VECTOR>
111  {
112  public:
118  void initialize(const BlockInfo *block_info,
137  template <class DOFINFO>
138  void initialize_info(DOFINFO &info, bool face) const;
139 
140 
145  template<class DOFINFO>
146  void assemble(const DOFINFO &info);
147 
152  template<class DOFINFO>
153  void assemble(const DOFINFO &info1,
154  const DOFINFO &info2);
155  private:
160  void assemble(VECTOR &global,
161  const BlockVector<double> &local,
162  const std::vector<types::global_dof_index> &dof);
163 
169  NamedData<SmartPointer<VECTOR,
171 
175  SmartPointer<const BlockInfo,
182  };
183 
184 
212  template <class MATRIX, typename number = double>
214  {
215  public:
223  MatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
224 
231  void initialize(const BlockInfo *block_info,
232  MatrixBlockVector<MATRIX> &matrices);
233 
251  template <class DOFINFO>
252  void initialize_info(DOFINFO &info, bool face) const;
253 
254 
259  template<class DOFINFO>
260  void assemble(const DOFINFO &info);
261 
266  template<class DOFINFO>
267  void assemble(const DOFINFO &info1,
268  const DOFINFO &info2);
269 
270  private:
275  void assemble(
276  MatrixBlock<MATRIX> &global,
277  const FullMatrix<number> &local,
278  const unsigned int block_row,
279  const unsigned int block_col,
280  const std::vector<types::global_dof_index> &dof1,
281  const std::vector<types::global_dof_index> &dof2);
282 
290 
294  SmartPointer<const BlockInfo,
301 
311  const double threshold;
312 
313  };
314 
339  template <class MATRIX, typename number = double>
341  {
342  public:
346 
354  MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
355 
362  void initialize(const BlockInfo *block_info,
363  MatrixPtrVector &matrices);
364 
369  void initialize(const MGConstrainedDoFs &mg_constrained_dofs);
370 
382  void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down);
383 
395  void initialize_interfaces (MatrixPtrVector &interface_in, MatrixPtrVector &interface_out);
409  template <class DOFINFO>
410  void initialize_info(DOFINFO &info, bool face) const;
411 
412 
417  template<class DOFINFO>
418  void assemble(const DOFINFO &info);
419 
424  template<class DOFINFO>
425  void assemble(const DOFINFO &info1,
426  const DOFINFO &info2);
427 
428  private:
433  void assemble(
434  MATRIX &global,
435  const FullMatrix<number> &local,
436  const unsigned int block_row,
437  const unsigned int block_col,
438  const std::vector<types::global_dof_index> &dof1,
439  const std::vector<types::global_dof_index> &dof2,
440  const unsigned int level1,
441  const unsigned int level2,
442  bool transpose = false);
443 
448  void assemble_fluxes(
449  MATRIX &global,
450  const FullMatrix<number> &local,
451  const unsigned int block_row,
452  const unsigned int block_col,
453  const std::vector<types::global_dof_index> &dof1,
454  const std::vector<types::global_dof_index> &dof2,
455  const unsigned int level1,
456  const unsigned int level2);
457 
462  void assemble_up(
463  MATRIX &global,
464  const FullMatrix<number> &local,
465  const unsigned int block_row,
466  const unsigned int block_col,
467  const std::vector<types::global_dof_index> &dof1,
468  const std::vector<types::global_dof_index> &dof2,
469  const unsigned int level1,
470  const unsigned int level2);
471 
476  void assemble_down(
477  MATRIX &global,
478  const FullMatrix<number> &local,
479  const unsigned int block_row,
480  const unsigned int block_col,
481  const std::vector<types::global_dof_index> &dof1,
482  const std::vector<types::global_dof_index> &dof2,
483  const unsigned int level1,
484  const unsigned int level2);
485 
490  void assemble_in(
491  MATRIX &global,
492  const FullMatrix<number> &local,
493  const unsigned int block_row,
494  const unsigned int block_col,
495  const std::vector<types::global_dof_index> &dof1,
496  const std::vector<types::global_dof_index> &dof2,
497  const unsigned int level1,
498  const unsigned int level2);
499 
504  void assemble_out(
505  MATRIX &global,
506  const FullMatrix<number> &local,
507  const unsigned int block_row,
508  const unsigned int block_col,
509  const std::vector<types::global_dof_index> &dof1,
510  const std::vector<types::global_dof_index> &dof2,
511  const unsigned int level1,
512  const unsigned int level2);
513 
519  MatrixPtrVectorPtr matrices;
520 
526  MatrixPtrVectorPtr flux_down;
527 
533  MatrixPtrVectorPtr flux_up;
534 
540  MatrixPtrVectorPtr interface_out;
541 
547  MatrixPtrVectorPtr interface_in;
548 
553 
558 
559 
569  const double threshold;
570 
571  };
572 
573 //----------------------------------------------------------------------//
574 
575  template <class VECTOR>
576  inline void
579  {
580  block_info = b;
581  residuals = m;
582  }
583 
584  template <class VECTOR>
585  inline void
587  const ConstraintMatrix &c)
588  {
589  constraints = &c;
590  }
591 
592 
593  template <class VECTOR>
594  template <class DOFINFO>
595  inline void
597  DOFINFO &info, bool) const
598  {
599  info.initialize_vectors(residuals.size());
600  }
601 
602  template <class VECTOR>
603  inline void
605  VECTOR &global,
606  const BlockVector<double> &local,
607  const std::vector<types::global_dof_index> &dof)
608  {
609  if (constraints == 0)
610  {
611  for (unsigned int b=0; b<local.n_blocks(); ++b)
612  for (unsigned int j=0; j<local.block(b).size(); ++j)
613  {
614  // The coordinates of
615  // the current entry in
616  // DoFHandler
617  // numbering, which
618  // differs from the
619  // block-wise local
620  // numbering we use in
621  // our local vectors
622  const unsigned int jcell = this->block_info->local().local_to_global(b, j);
623  global(dof[jcell]) += local.block(b)(j);
624  }
625  }
626  else
627  constraints->distribute_local_to_global(local, dof, global);
628  }
629 
630 
631  template <class VECTOR>
632  template <class DOFINFO>
633  inline void
635  const DOFINFO &info)
636  {
637  for (unsigned int i=0; i<residuals.size(); ++i)
638  assemble(*residuals(i), info.vector(i), info.indices);
639  }
640 
641 
642  template <class VECTOR>
643  template <class DOFINFO>
644  inline void
646  const DOFINFO &info1,
647  const DOFINFO &info2)
648  {
649  for (unsigned int i=0; i<residuals.size(); ++i)
650  {
651  assemble(*residuals(i), info1.vector(i), info1.indices);
652  assemble(*residuals(i), info2.vector(i), info2.indices);
653  }
654  }
655 
656 
657 //----------------------------------------------------------------------//
658 
659  template <class MATRIX, typename number>
660  inline
662  double threshold)
663  :
664  threshold(threshold)
665  {}
666 
667 
668  template <class MATRIX, typename number>
669  inline void
671  const BlockInfo *b,
673  {
674  block_info = b;
675  matrices = &m;
676  }
677 
678 
679 
680  template <class MATRIX, typename number>
681  inline void
683  const ConstraintMatrix &c)
684  {
685  constraints = &c;
686  }
687 
688 
689 
690  template <class MATRIX ,typename number>
691  template <class DOFINFO>
692  inline void
694  DOFINFO &info,
695  bool face) const
696  {
697  info.initialize_matrices(*matrices, face);
698  }
699 
700 
701 
702  template <class MATRIX, typename number>
703  inline void
705  MatrixBlock<MATRIX> &global,
706  const FullMatrix<number> &local,
707  const unsigned int block_row,
708  const unsigned int block_col,
709  const std::vector<types::global_dof_index> &dof1,
710  const std::vector<types::global_dof_index> &dof2)
711  {
712  if (constraints == 0)
713  {
714  for (unsigned int j=0; j<local.n_rows(); ++j)
715  for (unsigned int k=0; k<local.n_cols(); ++k)
716  if (std::fabs(local(j,k)) >= threshold)
717  {
718  // The coordinates of
719  // the current entry in
720  // DoFHandler
721  // numbering, which
722  // differs from the
723  // block-wise local
724  // numbering we use in
725  // our local matrices
726  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
727  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
728 
729  global.add(dof1[jcell], dof2[kcell], local(j,k));
730  }
731  }
732  else
733  {
734  const BlockIndices &bi = this->block_info->local();
735  std::vector<types::global_dof_index> sliced_row_indices (bi.block_size(block_row));
736  for (unsigned int i=0; i<sliced_row_indices.size(); ++i)
737  sliced_row_indices[i] = dof1[bi.block_start(block_row)+i];
738 
739  std::vector<types::global_dof_index> sliced_col_indices (bi.block_size(block_col));
740  for (unsigned int i=0; i<sliced_col_indices.size(); ++i)
741  sliced_col_indices[i] = dof2[bi.block_start(block_col)+i];
742 
743  constraints->distribute_local_to_global(local,
744  sliced_row_indices, sliced_col_indices, global);
745  }
746  }
747 
748 
749  template <class MATRIX, typename number>
750  template <class DOFINFO>
751  inline void
753  const DOFINFO &info)
754  {
755  for (unsigned int i=0; i<matrices->size(); ++i)
756  {
757  // Row and column index of
758  // the block we are dealing with
759  const types::global_dof_index row = matrices->block(i).row;
760  const types::global_dof_index col = matrices->block(i).column;
761 
762  assemble(matrices->block(i), info.matrix(i,false).matrix, row, col, info.indices, info.indices);
763  }
764  }
765 
766 
767  template <class MATRIX, typename number>
768  template <class DOFINFO>
769  inline void
771  const DOFINFO &info1,
772  const DOFINFO &info2)
773  {
774  for (unsigned int i=0; i<matrices->size(); ++i)
775  {
776  // Row and column index of
777  // the block we are dealing with
778  const types::global_dof_index row = matrices->block(i).row;
779  const types::global_dof_index col = matrices->block(i).column;
780 
781  assemble(matrices->block(i), info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices);
782  assemble(matrices->block(i), info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices);
783  assemble(matrices->block(i), info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices);
784  assemble(matrices->block(i), info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices);
785  }
786  }
787 
788 
789 // ----------------------------------------------------------------------//
790 
791  template <class MATRIX, typename number>
792  inline
794  double threshold)
795  :
796  threshold(threshold)
797  {}
798 
799 
800  template <class MATRIX, typename number>
801  inline void
803  const BlockInfo *b,
804  MatrixPtrVector &m)
805  {
806  block_info = b;
807  AssertDimension(block_info->local().size(), block_info->global().size());
808  matrices = &m;
809  }
810 
811 
812  template <class MATRIX, typename number>
813  inline void
815  const MGConstrainedDoFs &mg_c)
816  {
817  mg_constrained_dofs = &mg_c;
818  }
819 
820 
821  template <class MATRIX ,typename number>
822  template <class DOFINFO>
823  inline void
825  DOFINFO &info,
826  bool face) const
827  {
828  info.initialize_matrices(*matrices, face);
829  }
830 
831 
832 
833  template <class MATRIX, typename number>
834  inline void
836  MatrixPtrVector &up,
837  MatrixPtrVector &down)
838  {
839  flux_up = up;
840  flux_down = down;
841  }
842 
843 
844  template <class MATRIX, typename number>
845  inline void
847  MatrixPtrVector &in,
848  MatrixPtrVector &out)
849  {
850  interface_in = in;
851  interface_out = out;
852  }
853 
854 
855  template <class MATRIX, typename number>
856  inline void
858  MATRIX &global,
859  const FullMatrix<number> &local,
860  const unsigned int block_row,
861  const unsigned int block_col,
862  const std::vector<types::global_dof_index> &dof1,
863  const std::vector<types::global_dof_index> &dof2,
864  const unsigned int level1,
865  const unsigned int level2,
866  bool transpose)
867  {
868  for (unsigned int j=0; j<local.n_rows(); ++j)
869  for (unsigned int k=0; k<local.n_cols(); ++k)
870  if (std::fabs(local(j,k)) >= threshold)
871  {
872  // The coordinates of
873  // the current entry in
874  // DoFHandler
875  // numbering, which
876  // differs from the
877  // block-wise local
878  // numbering we use in
879  // our local matrices
880  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
881  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
882 
883  // The global dof
884  // indices to assemble
885  // in. Since we may
886  // have face matrices
887  // coupling two
888  // different cells, we
889  // provide two sets of
890  // dof indices.
891  const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
892  const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
893 
894  if (mg_constrained_dofs == 0)
895  {
896  if (transpose)
897  global.add(kglobal, jglobal, local(j,k));
898  else
899  global.add(jglobal, kglobal, local(j,k));
900  }
901  else
902  {
903  if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
904  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
905  {
906  if (mg_constrained_dofs->set_boundary_values())
907  {
908  if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
909  !mg_constrained_dofs->is_boundary_index(level2, kglobal))
910  ||
911  (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
912  mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
913  jglobal == kglobal))
914  {
915  if (transpose)
916  global.add(kglobal, jglobal, local(j,k));
917  else
918  global.add(jglobal, kglobal, local(j,k));
919  }
920  }
921  else
922  {
923  if (transpose)
924  global.add(kglobal, jglobal, local(j,k));
925  else
926  global.add(jglobal, kglobal, local(j,k));
927  }
928  }
929  }
930  }
931  }
932 
933 
934  template <class MATRIX, typename number>
935  inline void
937  MATRIX &global,
938  const FullMatrix<number> &local,
939  const unsigned int block_row,
940  const unsigned int block_col,
941  const std::vector<types::global_dof_index> &dof1,
942  const std::vector<types::global_dof_index> &dof2,
943  const unsigned int level1,
944  const unsigned int level2)
945  {
946  for (unsigned int j=0; j<local.n_rows(); ++j)
947  for (unsigned int k=0; k<local.n_cols(); ++k)
948  if (std::fabs(local(j,k)) >= threshold)
949  {
950  // The coordinates of
951  // the current entry in
952  // DoFHandler
953  // numbering, which
954  // differs from the
955  // block-wise local
956  // numbering we use in
957  // our local matrices
958  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
959  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
960 
961  // The global dof
962  // indices to assemble
963  // in. Since we may
964  // have face matrices
965  // coupling two
966  // different cells, we
967  // provide two sets of
968  // dof indices.
969  const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
970  const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
971 
972  if (mg_constrained_dofs == 0)
973  global.add(jglobal, kglobal, local(j,k));
974  else
975  {
976  if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
977  !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
978  {
979  if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
980  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
981  global.add(jglobal, kglobal, local(j,k));
982  }
983  }
984  }
985  }
986 
987  template <class MATRIX, typename number>
988  inline void
990  MATRIX &global,
991  const FullMatrix<number> &local,
992  const unsigned int block_row,
993  const unsigned int block_col,
994  const std::vector<types::global_dof_index> &dof1,
995  const std::vector<types::global_dof_index> &dof2,
996  const unsigned int level1,
997  const unsigned int level2)
998  {
999  for (unsigned int j=0; j<local.n_rows(); ++j)
1000  for (unsigned int k=0; k<local.n_cols(); ++k)
1001  if (std::fabs(local(j,k)) >= threshold)
1002  {
1003  // The coordinates of
1004  // the current entry in
1005  // DoFHandler
1006  // numbering, which
1007  // differs from the
1008  // block-wise local
1009  // numbering we use in
1010  // our local matrices
1011  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1012  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1013 
1014  // The global dof
1015  // indices to assemble
1016  // in. Since we may
1017  // have face matrices
1018  // coupling two
1019  // different cells, we
1020  // provide two sets of
1021  // dof indices.
1022  const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1023  const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1024 
1025  if (mg_constrained_dofs == 0)
1026  global.add(jglobal, kglobal, local(j,k));
1027  else
1028  {
1029  if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
1030  !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
1031  {
1032  if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1033  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1034  global.add(jglobal, kglobal, local(j,k));
1035  }
1036  }
1037  }
1038  }
1039 
1040  template <class MATRIX, typename number>
1041  inline void
1043  MATRIX &global,
1044  const FullMatrix<number> &local,
1045  const unsigned int block_row,
1046  const unsigned int block_col,
1047  const std::vector<types::global_dof_index> &dof1,
1048  const std::vector<types::global_dof_index> &dof2,
1049  const unsigned int level1,
1050  const unsigned int level2)
1051  {
1052  for (unsigned int j=0; j<local.n_rows(); ++j)
1053  for (unsigned int k=0; k<local.n_cols(); ++k)
1054  if (std::fabs(local(k,j)) >= threshold)
1055  {
1056  // The coordinates of
1057  // the current entry in
1058  // DoFHandler
1059  // numbering, which
1060  // differs from the
1061  // block-wise local
1062  // numbering we use in
1063  // our local matrices
1064  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1065  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1066 
1067  // The global dof
1068  // indices to assemble
1069  // in. Since we may
1070  // have face matrices
1071  // coupling two
1072  // different cells, we
1073  // provide two sets of
1074  // dof indices.
1075  const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1076  const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1077 
1078  if (mg_constrained_dofs == 0)
1079  global.add(jglobal, kglobal, local(k,j));
1080  else
1081  {
1082  if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
1083  !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
1084  {
1085  if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1086  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1087  global.add(jglobal, kglobal, local(k,j));
1088  }
1089  }
1090  }
1091  }
1092 
1093  template <class MATRIX, typename number>
1094  inline void
1096  MATRIX &global,
1097  const FullMatrix<number> &local,
1098  const unsigned int block_row,
1099  const unsigned int block_col,
1100  const std::vector<types::global_dof_index> &dof1,
1101  const std::vector<types::global_dof_index> &dof2,
1102  const unsigned int level1,
1103  const unsigned int level2)
1104  {
1105 // AssertDimension(local.n(), dof1.size());
1106 // AssertDimension(local.m(), dof2.size());
1107 
1108  for (unsigned int j=0; j<local.n_rows(); ++j)
1109  for (unsigned int k=0; k<local.n_cols(); ++k)
1110  if (std::fabs(local(j,k)) >= threshold)
1111  {
1112  // The coordinates of
1113  // the current entry in
1114  // DoFHandler
1115  // numbering, which
1116  // differs from the
1117  // block-wise local
1118  // numbering we use in
1119  // our local matrices
1120  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1121  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1122 
1123  // The global dof
1124  // indices to assemble
1125  // in. Since we may
1126  // have face matrices
1127  // coupling two
1128  // different cells, we
1129  // provide two sets of
1130  // dof indices.
1131  const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1132  const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1133 
1134  if (mg_constrained_dofs == 0)
1135  global.add(jglobal, kglobal, local(j,k));
1136  else
1137  {
1138  if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1139  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1140  {
1141  if (mg_constrained_dofs->set_boundary_values())
1142  {
1143  if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1144  !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
1145  ||
1146  (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1147  mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
1148  jglobal == kglobal))
1149  global.add(jglobal, kglobal, local(j,k));
1150  }
1151  else
1152  global.add(jglobal, kglobal, local(j,k));
1153  }
1154  }
1155  }
1156  }
1157 
1158  template <class MATRIX, typename number>
1159  inline void
1161  MATRIX &global,
1162  const FullMatrix<number> &local,
1163  const unsigned int block_row,
1164  const unsigned int block_col,
1165  const std::vector<types::global_dof_index> &dof1,
1166  const std::vector<types::global_dof_index> &dof2,
1167  const unsigned int level1,
1168  const unsigned int level2)
1169  {
1170 // AssertDimension(local.n(), dof1.size());
1171 // AssertDimension(local.m(), dof2.size());
1172 
1173  for (unsigned int j=0; j<local.n_rows(); ++j)
1174  for (unsigned int k=0; k<local.n_cols(); ++k)
1175  if (std::fabs(local(k,j)) >= threshold)
1176  {
1177  // The coordinates of
1178  // the current entry in
1179  // DoFHandler
1180  // numbering, which
1181  // differs from the
1182  // block-wise local
1183  // numbering we use in
1184  // our local matrices
1185  const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1186  const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1187 
1188  // The global dof
1189  // indices to assemble
1190  // in. Since we may
1191  // have face matrices
1192  // coupling two
1193  // different cells, we
1194  // provide two sets of
1195  // dof indices.
1196  const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1197  const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1198 
1199  if (mg_constrained_dofs == 0)
1200  global.add(jglobal, kglobal, local(k,j));
1201  else
1202  {
1203  if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1204  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1205  {
1206  if (mg_constrained_dofs->set_boundary_values())
1207  {
1208  if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1209  !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
1210  ||
1211  (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1212  mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
1213  jglobal == kglobal))
1214  global.add(jglobal, kglobal, local(k,j));
1215  }
1216  else
1217  global.add(jglobal, kglobal, local(k,j));
1218  }
1219  }
1220  }
1221  }
1222 
1223 
1224  template <class MATRIX, typename number>
1225  template <class DOFINFO>
1226  inline void
1228  {
1229  const unsigned int level = info.cell->level();
1230 
1231  for (unsigned int i=0; i<matrices->size(); ++i)
1232  {
1233  // Row and column index of
1234  // the block we are dealing with
1235  const unsigned int row = matrices->block(i)[level].row;
1236  const unsigned int col = matrices->block(i)[level].column;
1237 
1238  assemble(matrices->block(i)[level].matrix, info.matrix(i,false).matrix, row, col,
1239  info.indices, info.indices, level, level);
1240  if (mg_constrained_dofs != 0)
1241  {
1242  if (interface_in != 0)
1243  assemble_in(interface_in->block(i)[level], info.matrix(i,false).matrix, row, col,
1244  info.indices, info.indices, level, level);
1245  if (interface_out != 0)
1246  assemble_in(interface_out->block(i)[level], info.matrix(i,false).matrix, row, col,
1247  info.indices, info.indices, level, level);
1248 
1249  assemble_in(matrices->block_in(i)[level], info.matrix(i,false).matrix, row, col,
1250  info.indices, info.indices, level, level);
1251  assemble_out(matrices->block_out(i)[level], info.matrix(i,false).matrix, row, col,
1252  info.indices, info.indices, level, level);
1253  }
1254  }
1255  }
1256 
1257 
1258  template <class MATRIX, typename number>
1259  template <class DOFINFO>
1260  inline void
1262  const DOFINFO &info1,
1263  const DOFINFO &info2)
1264  {
1265  const unsigned int level1 = info1.cell->level();
1266  const unsigned int level2 = info2.cell->level();
1267 
1268  for (unsigned int i=0; i<matrices->size(); ++i)
1269  {
1271 
1272  // Row and column index of
1273  // the block we are dealing with
1274  const unsigned int row = o[level1].row;
1275  const unsigned int col = o[level1].column;
1276 
1277  if (level1 == level2)
1278  {
1279  if (mg_constrained_dofs == 0)
1280  {
1281  assemble(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
1282  info1.indices, info1.indices, level1, level1);
1283  assemble(o[level1].matrix, info1.matrix(i,true).matrix, row, col,
1284  info1.indices, info2.indices, level1, level2);
1285  assemble(o[level1].matrix, info2.matrix(i,false).matrix, row, col,
1286  info2.indices, info2.indices, level2, level2);
1287  assemble(o[level1].matrix, info2.matrix(i,true).matrix, row, col,
1288  info2.indices, info1.indices, level2, level1);
1289  }
1290  else
1291  {
1292  assemble_fluxes(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
1293  info1.indices, info1.indices, level1, level1);
1294  assemble_fluxes(o[level1].matrix, info1.matrix(i,true).matrix, row, col,
1295  info1.indices, info2.indices, level1, level2);
1296  assemble_fluxes(o[level1].matrix, info2.matrix(i,false).matrix, row, col,
1297  info2.indices, info2.indices, level2, level2);
1298  assemble_fluxes(o[level1].matrix, info2.matrix(i,true).matrix, row, col,
1299  info2.indices, info1.indices, level2, level1);
1300  }
1301  }
1302  else
1303  {
1304  Assert(level1 > level2, ExcNotImplemented());
1305  if (flux_up->size() != 0)
1306  {
1307  // Do not add M22,
1308  // which is done by
1309  // the coarser cell
1310  assemble_fluxes(o[level1].matrix, info1.matrix(i,false).matrix, row, col,
1311  info1.indices, info1.indices, level1, level1);
1312  assemble_up(flux_up->block(i)[level1].matrix, info1.matrix(i,true).matrix, row, col,
1313  info1.indices, info2.indices, level1, level2);
1314  assemble_down(flux_down->block(i)[level1].matrix, info2.matrix(i,true).matrix, row, col,
1315  info2.indices, info1.indices, level2, level1);
1316  }
1317  }
1318  }
1319  }
1320  }
1321 }
1322 
1323 DEAL_II_NAMESPACE_CLOSE
1324 
1325 #endif
void assemble_in(MATRIX &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition: assembler.h:1095
SmartPointer< const BlockInfo, MatrixLocalBlocksToGlobalBlocks< MATRIX, number > > block_info
Definition: assembler.h:295
const value_type & block_in(size_type i) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
MeshWorker::Assember objects for systems without block structure.
void initialize_info(DOFINFO &info, bool face) const
Definition: assembler.h:596
SmartPointer< MatrixBlockVector< MATRIX >, MatrixLocalBlocksToGlobalBlocks< MATRIX, number > > matrices
Definition: assembler.h:289
unsigned int size() const
MatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
Definition: assembler.h:661
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:55
void add(const size_type i, const size_type j, const typename MATRIX::value_type value)
Definition: matrix_block.h:768
SmartPointer< const ConstraintMatrix, ResidualLocalBlocksToGlobalBlocks< VECTOR > > constraints
Definition: assembler.h:181
void initialize_interfaces(MatrixPtrVector &interface_in, MatrixPtrVector &interface_out)
Definition: assembler.h:846
SmartPointer< const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks< MATRIX, number > > block_info
Definition: assembler.h:552
size_type block_start(const unsigned int i) const
void initialize_info(DOFINFO &info, bool face) const
Definition: assembler.h:693
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
void assemble_down(MATRIX &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition: assembler.h:1042
void assemble_fluxes(MATRIX &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition: assembler.h:936
SmartPointer< const ConstraintMatrix, MatrixLocalBlocksToGlobalBlocks< MATRIX, number > > constraints
Definition: assembler.h:300
void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down)
Definition: assembler.h:835
size_type block_size(const unsigned int i) const
void initialize(const BlockInfo *block_info, MatrixPtrVector &matrices)
Definition: assembler.h:802
NamedData< SmartPointer< VECTOR, ResidualLocalBlocksToGlobalBlocks< VECTOR > > > residuals
Definition: assembler.h:170
void initialize_info(DOFINFO &info, bool face) const
Definition: assembler.h:824
void assemble_out(MATRIX &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition: assembler.h:1160
const value_type & block(size_type i) const
::ExceptionBase & ExcNotImplemented()
A small class collecting the different BlockIndices involved in global, multilevel and local computat...
Definition: block_info.h:93
const value_type & block_out(size_type i) const
void initialize(const BlockInfo *block_info, MatrixBlockVector< MATRIX > &matrices)
Definition: assembler.h:670
BlockType & block(const unsigned int i)
void initialize(const BlockInfo *block_info, NamedData< VECTOR * > &residuals)
Definition: assembler.h:577
SmartPointer< const MGConstrainedDoFs, MGMatrixLocalBlocksToGlobalBlocks< MATRIX, number > > mg_constrained_dofs
Definition: assembler.h:557
SmartPointer< const BlockInfo, ResidualLocalBlocksToGlobalBlocks< VECTOR > > block_info
Definition: assembler.h:176
void assemble_up(MATRIX &global, const FullMatrix< number > &local, const unsigned int block_row, const unsigned int block_col, const std::vector< types::global_dof_index > &dof1, const std::vector< types::global_dof_index > &dof2, const unsigned int level1, const unsigned int level2)
Definition: assembler.h:989