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