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