Reference documentation for deal.II version Git e159ac89de 2020-06-06 19:38:41 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
assembler.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_mesh_worker_assembler_h
18 #define dealii_mesh_worker_assembler_h
19 
20 #include <deal.II/base/config.h>
21 
23 
26 
28 
32 
34 
35 
37 
38 namespace MeshWorker
39 {
88  namespace Assembler
89  {
110  template <typename VectorType>
112  {
113  public:
117  void
119 
123  void
124  initialize(
126 
134  template <class DOFINFO>
135  void
136  initialize_info(DOFINFO &info, bool face) const;
137 
138 
142  template <class DOFINFO>
143  void
144  assemble(const DOFINFO &info);
145 
149  template <class DOFINFO>
150  void
151  assemble(const DOFINFO &info1, const DOFINFO &info2);
152 
153  private:
157  void
158  assemble(VectorType & global,
159  const BlockVector<double> & local,
160  const std::vector<types::global_dof_index> &dof);
161 
166 
170  SmartPointer<const BlockInfo,
173 
180  };
181 
182 
209  template <typename MatrixType, typename number = double>
211  {
212  public:
217  MatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
218 
223  void
226 
230  void
231  initialize(
233 
241  template <class DOFINFO>
242  void
243  initialize_info(DOFINFO &info, bool face) const;
244 
245 
249  template <class DOFINFO>
250  void
251  assemble(const DOFINFO &info);
252 
256  template <class DOFINFO>
257  void
258  assemble(const DOFINFO &info1, const DOFINFO &info2);
259 
260  private:
264  void
266  const FullMatrix<number> & local,
267  const unsigned int block_row,
268  const unsigned int block_col,
269  const std::vector<types::global_dof_index> &dof1,
270  const std::vector<types::global_dof_index> &dof2);
271 
278 
282  SmartPointer<const BlockInfo,
292 
298  const double threshold;
299  };
300 
324  template <typename MatrixType, typename number = double>
326  {
327  public:
329  using MatrixPtrVectorPtr =
332 
337  MGMatrixLocalBlocksToGlobalBlocks(double threshold = 1.e-12);
338 
343  void
344  initialize(const BlockInfo *block_info, MatrixPtrVector &matrices);
345 
349  void
350  initialize(const MGConstrainedDoFs &mg_constrained_dofs);
351 
357  void
358  initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down);
359 
365  void
366  initialize_interfaces(MatrixPtrVector &interface_in,
367  MatrixPtrVector &interface_out);
375  template <class DOFINFO>
376  void
377  initialize_info(DOFINFO &info, bool face) const;
378 
379 
383  template <class DOFINFO>
384  void
385  assemble(const DOFINFO &info);
386 
390  template <class DOFINFO>
391  void
392  assemble(const DOFINFO &info1, const DOFINFO &info2);
393 
394  private:
398  void
399  assemble(MatrixType & global,
400  const FullMatrix<number> & local,
401  const unsigned int block_row,
402  const unsigned int block_col,
403  const std::vector<types::global_dof_index> &dof1,
404  const std::vector<types::global_dof_index> &dof2,
405  const unsigned int level1,
406  const unsigned int level2,
407  bool transpose = false);
408 
412  void
413  assemble_fluxes(MatrixType & global,
414  const FullMatrix<number> & local,
415  const unsigned int block_row,
416  const unsigned int block_col,
417  const std::vector<types::global_dof_index> &dof1,
418  const std::vector<types::global_dof_index> &dof2,
419  const unsigned int level1,
420  const unsigned int level2);
421 
425  void
426  assemble_up(MatrixType & global,
427  const FullMatrix<number> & local,
428  const unsigned int block_row,
429  const unsigned int block_col,
430  const std::vector<types::global_dof_index> &dof1,
431  const std::vector<types::global_dof_index> &dof2,
432  const unsigned int level1,
433  const unsigned int level2);
434 
438  void
439  assemble_down(MatrixType & global,
440  const FullMatrix<number> & local,
441  const unsigned int block_row,
442  const unsigned int block_col,
443  const std::vector<types::global_dof_index> &dof1,
444  const std::vector<types::global_dof_index> &dof2,
445  const unsigned int level1,
446  const unsigned int level2);
447 
451  void
452  assemble_in(MatrixType & global,
453  const FullMatrix<number> & local,
454  const unsigned int block_row,
455  const unsigned int block_col,
456  const std::vector<types::global_dof_index> &dof1,
457  const std::vector<types::global_dof_index> &dof2,
458  const unsigned int level1,
459  const unsigned int level2);
460 
464  void
465  assemble_out(MatrixType & global,
466  const FullMatrix<number> & local,
467  const unsigned int block_row,
468  const unsigned int block_col,
469  const std::vector<types::global_dof_index> &dof1,
470  const std::vector<types::global_dof_index> &dof2,
471  const unsigned int level1,
472  const unsigned int level2);
473 
478 
484 
490 
496 
502 
506  SmartPointer<const BlockInfo,
509 
516 
517 
523  const double threshold;
524  };
525 
526  //----------------------------------------------------------------------//
527 
528  template <typename VectorType>
529  inline void
531  const BlockInfo *b,
532  AnyData & m)
533  {
534  block_info = b;
535  residuals = m;
536  }
537 
538  template <typename VectorType>
539  inline void
542  {
543  constraints = &c;
544  }
545 
546 
547  template <typename VectorType>
548  template <class DOFINFO>
549  inline void
551  DOFINFO &info,
552  bool) const
553  {
554  info.initialize_vectors(residuals.size());
555  }
556 
557  template <typename VectorType>
558  inline void
560  VectorType & global,
561  const BlockVector<double> & local,
562  const std::vector<types::global_dof_index> &dof)
563  {
564  if (constraints == 0)
565  {
566  for (unsigned int b = 0; b < local.n_blocks(); ++b)
567  for (unsigned int j = 0; j < local.block(b).size(); ++j)
568  {
569  // The coordinates of
570  // the current entry in
571  // DoFHandler
572  // numbering, which
573  // differs from the
574  // block-wise local
575  // numbering we use in
576  // our local vectors
577  const unsigned int jcell =
578  this->block_info->local().local_to_global(b, j);
579  global(dof[jcell]) += local.block(b)(j);
580  }
581  }
582  else
583  constraints->distribute_local_to_global(local, dof, global);
584  }
585 
586 
587  template <typename VectorType>
588  template <class DOFINFO>
589  inline void
591  {
592  for (unsigned int i = 0; i < residuals.size(); ++i)
594  info.vector(i),
595  info.indices);
596  }
597 
598 
599  template <typename VectorType>
600  template <class DOFINFO>
601  inline void
603  const DOFINFO &info1,
604  const DOFINFO &info2)
605  {
606  for (unsigned int i = 0; i < residuals.size(); ++i)
607  {
609  info1.vector(i),
610  info1.indices);
612  info2.vector(i),
613  info2.indices);
614  }
615  }
616 
617 
618  //----------------------------------------------------------------------//
619 
620  template <typename MatrixType, typename number>
623  : threshold(threshold)
624  {}
625 
626 
627  template <typename MatrixType, typename number>
628  inline void
630  const BlockInfo * b,
632  {
633  block_info = b;
634  matrices = &m;
635  }
636 
637 
638 
639  template <typename MatrixType, typename number>
640  inline void
643  {
644  constraints = &c;
645  }
646 
647 
648 
649  template <typename MatrixType, typename number>
650  template <class DOFINFO>
651  inline void
653  DOFINFO &info,
654  bool face) const
655  {
656  info.initialize_matrices(*matrices, face);
657  }
658 
659 
660 
661  template <typename MatrixType, typename number>
662  inline void
664  MatrixBlock<MatrixType> & global,
665  const FullMatrix<number> & local,
666  const unsigned int block_row,
667  const unsigned int block_col,
668  const std::vector<types::global_dof_index> &dof1,
669  const std::vector<types::global_dof_index> &dof2)
670  {
671  if (constraints == nullptr)
672  {
673  for (unsigned int j = 0; j < local.n_rows(); ++j)
674  for (unsigned int k = 0; k < local.n_cols(); ++k)
675  if (std::fabs(local(j, k)) >= threshold)
676  {
677  // The coordinates of
678  // the current entry in
679  // DoFHandler
680  // numbering, which
681  // differs from the
682  // block-wise local
683  // numbering we use in
684  // our local matrices
685  const unsigned int jcell =
686  this->block_info->local().local_to_global(block_row, j);
687  const unsigned int kcell =
688  this->block_info->local().local_to_global(block_col, k);
689 
690  global.add(dof1[jcell], dof2[kcell], local(j, k));
691  }
692  }
693  else
694  {
695  const BlockIndices & bi = this->block_info->local();
696  std::vector<types::global_dof_index> sliced_row_indices(
697  bi.block_size(block_row));
698  for (unsigned int i = 0; i < sliced_row_indices.size(); ++i)
699  sliced_row_indices[i] = dof1[bi.block_start(block_row) + i];
700 
701  std::vector<types::global_dof_index> sliced_col_indices(
702  bi.block_size(block_col));
703  for (unsigned int i = 0; i < sliced_col_indices.size(); ++i)
704  sliced_col_indices[i] = dof2[bi.block_start(block_col) + i];
705 
706  constraints->distribute_local_to_global(local,
707  sliced_row_indices,
708  sliced_col_indices,
709  global);
710  }
711  }
712 
713 
714  template <typename MatrixType, typename number>
715  template <class DOFINFO>
716  inline void
718  const DOFINFO &info)
719  {
720  for (unsigned int i = 0; i < matrices->size(); ++i)
721  {
722  // Row and column index of
723  // the block we are dealing with
724  const types::global_dof_index row = matrices->block(i).row;
725  const types::global_dof_index col = matrices->block(i).column;
726 
727  assemble(matrices->block(i),
728  info.matrix(i, false).matrix,
729  row,
730  col,
731  info.indices,
732  info.indices);
733  }
734  }
735 
736 
737  template <typename MatrixType, typename number>
738  template <class DOFINFO>
739  inline void
741  const DOFINFO &info1,
742  const DOFINFO &info2)
743  {
744  for (unsigned int i = 0; i < matrices->size(); ++i)
745  {
746  // Row and column index of
747  // the block we are dealing with
748  const types::global_dof_index row = matrices->block(i).row;
749  const types::global_dof_index col = matrices->block(i).column;
750 
751  assemble(matrices->block(i),
752  info1.matrix(i, false).matrix,
753  row,
754  col,
755  info1.indices,
756  info1.indices);
757  assemble(matrices->block(i),
758  info1.matrix(i, true).matrix,
759  row,
760  col,
761  info1.indices,
762  info2.indices);
763  assemble(matrices->block(i),
764  info2.matrix(i, false).matrix,
765  row,
766  col,
767  info2.indices,
768  info2.indices);
769  assemble(matrices->block(i),
770  info2.matrix(i, true).matrix,
771  row,
772  col,
773  info2.indices,
774  info1.indices);
775  }
776  }
777 
778 
779  // ----------------------------------------------------------------------//
780 
781  template <typename MatrixType, typename number>
784  : threshold(threshold)
785  {}
786 
787 
788  template <typename MatrixType, typename number>
789  inline void
791  const BlockInfo *b,
792  MatrixPtrVector &m)
793  {
794  block_info = b;
795  AssertDimension(block_info->local().size(), block_info->global().size());
796  matrices = &m;
797  }
798 
799 
800  template <typename MatrixType, typename number>
801  inline void
803  const MGConstrainedDoFs &mg_c)
804  {
805  mg_constrained_dofs = &mg_c;
806  }
807 
808 
809  template <typename MatrixType, typename number>
810  template <class DOFINFO>
811  inline void
813  DOFINFO &info,
814  bool face) const
815  {
816  info.initialize_matrices(*matrices, face);
817  }
818 
819 
820 
821  template <typename MatrixType, typename number>
822  inline void
824  MatrixPtrVector &up,
825  MatrixPtrVector &down)
826  {
827  flux_up = up;
828  flux_down = down;
829  }
830 
831 
832  template <typename MatrixType, typename number>
833  inline void
836  {
837  interface_in = in;
838  interface_out = out;
839  }
840 
841 
842  template <typename MatrixType, typename number>
843  inline void
845  MatrixType & global,
846  const FullMatrix<number> & local,
847  const unsigned int block_row,
848  const unsigned int block_col,
849  const std::vector<types::global_dof_index> &dof1,
850  const std::vector<types::global_dof_index> &dof2,
851  const unsigned int level1,
852  const unsigned int level2,
853  bool transpose)
854  {
855  for (unsigned int j = 0; j < local.n_rows(); ++j)
856  for (unsigned int k = 0; k < local.n_cols(); ++k)
857  if (std::fabs(local(j, k)) >= threshold)
858  {
859  // The coordinates of
860  // the current entry in
861  // DoFHandler
862  // numbering, which
863  // differs from the
864  // block-wise local
865  // numbering we use in
866  // our local matrices
867  const unsigned int jcell =
868  this->block_info->local().local_to_global(block_row, j);
869  const unsigned int kcell =
870  this->block_info->local().local_to_global(block_col, k);
871 
872  // The global dof
873  // indices to assemble
874  // in. Since we may
875  // have face matrices
876  // coupling two
877  // different cells, we
878  // provide two sets of
879  // dof indices.
880  const unsigned int jglobal = this->block_info->level(level1)
881  .global_to_local(dof1[jcell])
882  .second;
883  const unsigned int kglobal = this->block_info->level(level2)
884  .global_to_local(dof2[kcell])
885  .second;
886 
887  if (mg_constrained_dofs == 0)
888  {
889  if (transpose)
890  global.add(kglobal, jglobal, local(j, k));
891  else
892  global.add(jglobal, kglobal, local(j, k));
893  }
894  else
895  {
896  if (!mg_constrained_dofs->at_refinement_edge(level1,
897  jglobal) &&
898  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
899  {
900  if (mg_constrained_dofs->set_boundary_values())
901  {
902  if ((!mg_constrained_dofs->is_boundary_index(
903  level1, jglobal) &&
904  !mg_constrained_dofs->is_boundary_index(
905  level2, kglobal)) ||
906  (mg_constrained_dofs->is_boundary_index(
907  level1, jglobal) &&
908  mg_constrained_dofs->is_boundary_index(
909  level2, kglobal) &&
910  jglobal == kglobal))
911  {
912  if (transpose)
913  global.add(kglobal, jglobal, local(j, k));
914  else
915  global.add(jglobal, kglobal, local(j, k));
916  }
917  }
918  else
919  {
920  if (transpose)
921  global.add(kglobal, jglobal, local(j, k));
922  else
923  global.add(jglobal, kglobal, local(j, k));
924  }
925  }
926  }
927  }
928  }
929 
930 
931  template <typename MatrixType, typename number>
932  inline void
934  MatrixType & global,
935  const FullMatrix<number> & local,
936  const unsigned int block_row,
937  const unsigned int block_col,
938  const std::vector<types::global_dof_index> &dof1,
939  const std::vector<types::global_dof_index> &dof2,
940  const unsigned int level1,
941  const unsigned int level2)
942  {
943  for (unsigned int j = 0; j < local.n_rows(); ++j)
944  for (unsigned int k = 0; k < local.n_cols(); ++k)
945  if (std::fabs(local(j, k)) >= threshold)
946  {
947  // The coordinates of
948  // the current entry in
949  // DoFHandler
950  // numbering, which
951  // differs from the
952  // block-wise local
953  // numbering we use in
954  // our local matrices
955  const unsigned int jcell =
956  this->block_info->local().local_to_global(block_row, j);
957  const unsigned int kcell =
958  this->block_info->local().local_to_global(block_col, k);
959 
960  // The global dof
961  // indices to assemble
962  // in. Since we may
963  // have face matrices
964  // coupling two
965  // different cells, we
966  // provide two sets of
967  // dof indices.
968  const unsigned int jglobal = this->block_info->level(level1)
969  .global_to_local(dof1[jcell])
970  .second;
971  const unsigned int kglobal = this->block_info->level(level2)
972  .global_to_local(dof2[kcell])
973  .second;
974 
975  if (mg_constrained_dofs == 0)
976  global.add(jglobal, kglobal, local(j, k));
977  else
978  {
979  if (!mg_constrained_dofs->non_refinement_edge_index(
980  level1, jglobal) &&
981  !mg_constrained_dofs->non_refinement_edge_index(level2,
982  kglobal))
983  {
984  if (!mg_constrained_dofs->at_refinement_edge(level1,
985  jglobal) &&
986  !mg_constrained_dofs->at_refinement_edge(level2,
987  kglobal))
988  global.add(jglobal, kglobal, local(j, k));
989  }
990  }
991  }
992  }
993 
994  template <typename MatrixType, typename number>
995  inline void
997  MatrixType & global,
998  const FullMatrix<number> & local,
999  const unsigned int block_row,
1000  const unsigned int block_col,
1001  const std::vector<types::global_dof_index> &dof1,
1002  const std::vector<types::global_dof_index> &dof2,
1003  const unsigned int level1,
1004  const unsigned int level2)
1005  {
1006  for (unsigned int j = 0; j < local.n_rows(); ++j)
1007  for (unsigned int k = 0; k < local.n_cols(); ++k)
1008  if (std::fabs(local(j, k)) >= threshold)
1009  {
1010  // The coordinates of
1011  // the current entry in
1012  // DoFHandler
1013  // numbering, which
1014  // differs from the
1015  // block-wise local
1016  // numbering we use in
1017  // our local matrices
1018  const unsigned int jcell =
1019  this->block_info->local().local_to_global(block_row, j);
1020  const unsigned int kcell =
1021  this->block_info->local().local_to_global(block_col, k);
1022 
1023  // The global dof
1024  // indices to assemble
1025  // in. Since we may
1026  // have face matrices
1027  // coupling two
1028  // different cells, we
1029  // provide two sets of
1030  // dof indices.
1031  const unsigned int jglobal = this->block_info->level(level1)
1032  .global_to_local(dof1[jcell])
1033  .second;
1034  const unsigned int kglobal = this->block_info->level(level2)
1035  .global_to_local(dof2[kcell])
1036  .second;
1037 
1038  if (mg_constrained_dofs == 0)
1039  global.add(jglobal, kglobal, local(j, k));
1040  else
1041  {
1042  if (!mg_constrained_dofs->non_refinement_edge_index(
1043  level1, jglobal) &&
1044  !mg_constrained_dofs->non_refinement_edge_index(level2,
1045  kglobal))
1046  {
1047  if (!mg_constrained_dofs->at_refinement_edge(level1,
1048  jglobal) &&
1049  !mg_constrained_dofs->at_refinement_edge(level2,
1050  kglobal))
1051  global.add(jglobal, kglobal, local(j, k));
1052  }
1053  }
1054  }
1055  }
1056 
1057  template <typename MatrixType, typename number>
1058  inline void
1060  MatrixType & global,
1061  const FullMatrix<number> & local,
1062  const unsigned int block_row,
1063  const unsigned int block_col,
1064  const std::vector<types::global_dof_index> &dof1,
1065  const std::vector<types::global_dof_index> &dof2,
1066  const unsigned int level1,
1067  const unsigned int level2)
1068  {
1069  for (unsigned int j = 0; j < local.n_rows(); ++j)
1070  for (unsigned int k = 0; k < local.n_cols(); ++k)
1071  if (std::fabs(local(k, j)) >= threshold)
1072  {
1073  // The coordinates of
1074  // the current entry in
1075  // DoFHandler
1076  // numbering, which
1077  // differs from the
1078  // block-wise local
1079  // numbering we use in
1080  // our local matrices
1081  const unsigned int jcell =
1082  this->block_info->local().local_to_global(block_row, j);
1083  const unsigned int kcell =
1084  this->block_info->local().local_to_global(block_col, k);
1085 
1086  // The global dof
1087  // indices to assemble
1088  // in. Since we may
1089  // have face matrices
1090  // coupling two
1091  // different cells, we
1092  // provide two sets of
1093  // dof indices.
1094  const unsigned int jglobal = this->block_info->level(level1)
1095  .global_to_local(dof1[jcell])
1096  .second;
1097  const unsigned int kglobal = this->block_info->level(level2)
1098  .global_to_local(dof2[kcell])
1099  .second;
1100 
1101  if (mg_constrained_dofs == 0)
1102  global.add(jglobal, kglobal, local(k, j));
1103  else
1104  {
1105  if (!mg_constrained_dofs->non_refinement_edge_index(
1106  level1, jglobal) &&
1107  !mg_constrained_dofs->non_refinement_edge_index(level2,
1108  kglobal))
1109  {
1110  if (!mg_constrained_dofs->at_refinement_edge(level1,
1111  jglobal) &&
1112  !mg_constrained_dofs->at_refinement_edge(level2,
1113  kglobal))
1114  global.add(jglobal, kglobal, local(k, j));
1115  }
1116  }
1117  }
1118  }
1119 
1120  template <typename MatrixType, typename number>
1121  inline void
1123  MatrixType & global,
1124  const FullMatrix<number> & local,
1125  const unsigned int block_row,
1126  const unsigned int block_col,
1127  const std::vector<types::global_dof_index> &dof1,
1128  const std::vector<types::global_dof_index> &dof2,
1129  const unsigned int level1,
1130  const unsigned int level2)
1131  {
1132  // AssertDimension(local.n(), dof1.size());
1133  // AssertDimension(local.m(), dof2.size());
1134 
1135  for (unsigned int j = 0; j < local.n_rows(); ++j)
1136  for (unsigned int k = 0; k < local.n_cols(); ++k)
1137  if (std::fabs(local(j, k)) >= threshold)
1138  {
1139  // The coordinates of
1140  // the current entry in
1141  // DoFHandler
1142  // numbering, which
1143  // differs from the
1144  // block-wise local
1145  // numbering we use in
1146  // our local matrices
1147  const unsigned int jcell =
1148  this->block_info->local().local_to_global(block_row, j);
1149  const unsigned int kcell =
1150  this->block_info->local().local_to_global(block_col, k);
1151 
1152  // The global dof
1153  // indices to assemble
1154  // in. Since we may
1155  // have face matrices
1156  // coupling two
1157  // different cells, we
1158  // provide two sets of
1159  // dof indices.
1160  const unsigned int jglobal = this->block_info->level(level1)
1161  .global_to_local(dof1[jcell])
1162  .second;
1163  const unsigned int kglobal = this->block_info->level(level2)
1164  .global_to_local(dof2[kcell])
1165  .second;
1166 
1167  if (mg_constrained_dofs == 0)
1168  global.add(jglobal, kglobal, local(j, k));
1169  else
1170  {
1171  if (mg_constrained_dofs->at_refinement_edge(level1,
1172  jglobal) &&
1173  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1174  {
1175  if (mg_constrained_dofs->set_boundary_values())
1176  {
1177  if ((!mg_constrained_dofs->is_boundary_index(
1178  level1, jglobal) &&
1179  !mg_constrained_dofs->is_boundary_index(
1180  level2, kglobal)) ||
1181  (mg_constrained_dofs->is_boundary_index(
1182  level1, jglobal) &&
1183  mg_constrained_dofs->is_boundary_index(
1184  level2, kglobal) &&
1185  jglobal == kglobal))
1186  global.add(jglobal, kglobal, local(j, k));
1187  }
1188  else
1189  global.add(jglobal, kglobal, local(j, k));
1190  }
1191  }
1192  }
1193  }
1194 
1195  template <typename MatrixType, typename number>
1196  inline void
1198  MatrixType & global,
1199  const FullMatrix<number> & local,
1200  const unsigned int block_row,
1201  const unsigned int block_col,
1202  const std::vector<types::global_dof_index> &dof1,
1203  const std::vector<types::global_dof_index> &dof2,
1204  const unsigned int level1,
1205  const unsigned int level2)
1206  {
1207  // AssertDimension(local.n(), dof1.size());
1208  // AssertDimension(local.m(), dof2.size());
1209 
1210  for (unsigned int j = 0; j < local.n_rows(); ++j)
1211  for (unsigned int k = 0; k < local.n_cols(); ++k)
1212  if (std::fabs(local(k, j)) >= threshold)
1213  {
1214  // The coordinates of
1215  // the current entry in
1216  // DoFHandler
1217  // numbering, which
1218  // differs from the
1219  // block-wise local
1220  // numbering we use in
1221  // our local matrices
1222  const unsigned int jcell =
1223  this->block_info->local().local_to_global(block_row, j);
1224  const unsigned int kcell =
1225  this->block_info->local().local_to_global(block_col, k);
1226 
1227  // The global dof
1228  // indices to assemble
1229  // in. Since we may
1230  // have face matrices
1231  // coupling two
1232  // different cells, we
1233  // provide two sets of
1234  // dof indices.
1235  const unsigned int jglobal = this->block_info->level(level1)
1236  .global_to_local(dof1[jcell])
1237  .second;
1238  const unsigned int kglobal = this->block_info->level(level2)
1239  .global_to_local(dof2[kcell])
1240  .second;
1241 
1242  if (mg_constrained_dofs == 0)
1243  global.add(jglobal, kglobal, local(k, j));
1244  else
1245  {
1246  if (mg_constrained_dofs->at_refinement_edge(level1,
1247  jglobal) &&
1248  !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1249  {
1250  if (mg_constrained_dofs->set_boundary_values())
1251  {
1252  if ((!mg_constrained_dofs->is_boundary_index(
1253  level1, jglobal) &&
1254  !mg_constrained_dofs->is_boundary_index(
1255  level2, kglobal)) ||
1256  (mg_constrained_dofs->is_boundary_index(
1257  level1, jglobal) &&
1258  mg_constrained_dofs->is_boundary_index(
1259  level2, kglobal) &&
1260  jglobal == kglobal))
1261  global.add(jglobal, kglobal, local(k, j));
1262  }
1263  else
1264  global.add(jglobal, kglobal, local(k, j));
1265  }
1266  }
1267  }
1268  }
1269 
1270 
1271  template <typename MatrixType, typename number>
1272  template <class DOFINFO>
1273  inline void
1275  const DOFINFO &info)
1276  {
1277  const unsigned int level = info.cell->level();
1278 
1279  for (unsigned int i = 0; i < matrices->size(); ++i)
1280  {
1281  // Row and column index of
1282  // the block we are dealing with
1283  const unsigned int row = matrices->block(i)[level].row;
1284  const unsigned int col = matrices->block(i)[level].column;
1285 
1286  assemble(matrices->block(i)[level].matrix,
1287  info.matrix(i, false).matrix,
1288  row,
1289  col,
1290  info.indices,
1291  info.indices,
1292  level,
1293  level);
1294  if (mg_constrained_dofs != 0)
1295  {
1296  if (interface_in != 0)
1298  info.matrix(i, false).matrix,
1299  row,
1300  col,
1301  info.indices,
1302  info.indices,
1303  level,
1304  level);
1305  if (interface_out != 0)
1307  info.matrix(i, false).matrix,
1308  row,
1309  col,
1310  info.indices,
1311  info.indices,
1312  level,
1313  level);
1314 
1316  info.matrix(i, false).matrix,
1317  row,
1318  col,
1319  info.indices,
1320  info.indices,
1321  level,
1322  level);
1324  info.matrix(i, false).matrix,
1325  row,
1326  col,
1327  info.indices,
1328  info.indices,
1329  level,
1330  level);
1331  }
1332  }
1333  }
1334 
1335 
1336  template <typename MatrixType, typename number>
1337  template <class DOFINFO>
1338  inline void
1340  const DOFINFO &info1,
1341  const DOFINFO &info2)
1342  {
1343  const unsigned int level1 = info1.cell->level();
1344  const unsigned int level2 = info2.cell->level();
1345 
1346  for (unsigned int i = 0; i < matrices->size(); ++i)
1347  {
1349 
1350  // Row and column index of
1351  // the block we are dealing with
1352  const unsigned int row = o[level1].row;
1353  const unsigned int col = o[level1].column;
1354 
1355  if (level1 == level2)
1356  {
1357  if (mg_constrained_dofs == 0)
1358  {
1359  assemble(o[level1].matrix,
1360  info1.matrix(i, false).matrix,
1361  row,
1362  col,
1363  info1.indices,
1364  info1.indices,
1365  level1,
1366  level1);
1367  assemble(o[level1].matrix,
1368  info1.matrix(i, true).matrix,
1369  row,
1370  col,
1371  info1.indices,
1372  info2.indices,
1373  level1,
1374  level2);
1375  assemble(o[level1].matrix,
1376  info2.matrix(i, false).matrix,
1377  row,
1378  col,
1379  info2.indices,
1380  info2.indices,
1381  level2,
1382  level2);
1383  assemble(o[level1].matrix,
1384  info2.matrix(i, true).matrix,
1385  row,
1386  col,
1387  info2.indices,
1388  info1.indices,
1389  level2,
1390  level1);
1391  }
1392  else
1393  {
1394  assemble_fluxes(o[level1].matrix,
1395  info1.matrix(i, false).matrix,
1396  row,
1397  col,
1398  info1.indices,
1399  info1.indices,
1400  level1,
1401  level1);
1402  assemble_fluxes(o[level1].matrix,
1403  info1.matrix(i, true).matrix,
1404  row,
1405  col,
1406  info1.indices,
1407  info2.indices,
1408  level1,
1409  level2);
1410  assemble_fluxes(o[level1].matrix,
1411  info2.matrix(i, false).matrix,
1412  row,
1413  col,
1414  info2.indices,
1415  info2.indices,
1416  level2,
1417  level2);
1418  assemble_fluxes(o[level1].matrix,
1419  info2.matrix(i, true).matrix,
1420  row,
1421  col,
1422  info2.indices,
1423  info1.indices,
1424  level2,
1425  level1);
1426  }
1427  }
1428  else
1429  {
1430  Assert(level1 > level2, ExcNotImplemented());
1431  if (flux_up->size() != 0)
1432  {
1433  // Do not add M22,
1434  // which is done by
1435  // the coarser cell
1436  assemble_fluxes(o[level1].matrix,
1437  info1.matrix(i, false).matrix,
1438  row,
1439  col,
1440  info1.indices,
1441  info1.indices,
1442  level1,
1443  level1);
1444  assemble_up(flux_up->block(i)[level1].matrix,
1445  info1.matrix(i, true).matrix,
1446  row,
1447  col,
1448  info1.indices,
1449  info2.indices,
1450  level1,
1451  level2);
1452  assemble_down(flux_down->block(i)[level1].matrix,
1453  info2.matrix(i, true).matrix,
1454  row,
1455  col,
1456  info2.indices,
1457  info1.indices,
1458  level2,
1459  level1);
1460  }
1461  }
1462  }
1463  }
1464  } // namespace Assembler
1465 } // namespace MeshWorker
1466 
1468 
1469 #endif
SmartPointer< const BlockInfo, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
Definition: assembler.h:284
SmartPointer< const AffineConstraints< typename VectorType::value_type >, ResidualLocalBlocksToGlobalBlocks< VectorType > > constraints
Definition: assembler.h:179
void assemble_down(MatrixType &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:1059
type entry(const std::string &name)
Access to stored data object by name.
Definition: any_data.h:351
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down)
Definition: assembler.h:823
Contents is actually a matrix.
void initialize(const BlockInfo *block_info, AnyData &residuals)
Definition: assembler.h:530
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SmartPointer< const MGConstrainedDoFs, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > mg_constrained_dofs
Definition: assembler.h:515
void initialize_info(DOFINFO &info, bool face) const
Definition: assembler.h:652
const value_type & block(size_type i) const
Definition: matrix_block.h:952
SmartPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > constraints
Definition: assembler.h:291
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
Definition: matrix_block.h:675
size_type block_size(const unsigned int i) const
void initialize_info(DOFINFO &info, bool face) const
Definition: assembler.h:550
const value_type & block_out(size_type i) const
Definition: matrix_block.h:984
#define Assert(cond, exc)
Definition: exceptions.h:1403
void assemble_fluxes(MatrixType &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:933
SmartPointer< const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
Definition: assembler.h:508
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
unsigned int level
Definition: grid_out.cc:4355
Expression fabs(const Expression &x)
const value_type & block_in(size_type i) const
Definition: matrix_block.h:968
void initialize(const BlockInfo *block_info, MatrixBlockVector< MatrixType > &matrices)
Definition: assembler.h:629
SmartPointer< const BlockInfo, ResidualLocalBlocksToGlobalBlocks< VectorType > > block_info
Definition: assembler.h:172
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
unsigned int size() const
Definition: matrix_block.h:920
void assemble_out(MatrixType &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:1197
void initialize_interfaces(MatrixPtrVector &interface_in, MatrixPtrVector &interface_out)
Definition: assembler.h:835
unsigned int size() const
Number of stored data objects.
Definition: any_data.h:223
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
void assemble_up(MatrixType &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:996
size_type block_start(const unsigned int i) const
A small class collecting the different BlockIndices involved in global, multilevel and local computat...
Definition: block_info.h:98
static ::ExceptionBase & ExcNotImplemented()
SmartPointer< MatrixBlockVector< MatrixType >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > matrices
Definition: assembler.h:277
BlockType & block(const unsigned int i)
void initialize_info(DOFINFO &info, bool face) const
Definition: assembler.h:812
MatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
Definition: assembler.h:622
void initialize(const BlockInfo *block_info, MatrixPtrVector &matrices)
Definition: assembler.h:790
void assemble_in(MatrixType &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:1122