18 #ifndef __deal2__mesh_worker_simple_h
19 #define __deal2__mesh_worker_simple_h
21 #include <deal.II/base/named_data.h>
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/base/mg_level_object.h>
24 #include <deal.II/lac/block_vector.h>
25 #include <deal.II/meshworker/dof_info.h>
26 #include <deal.II/meshworker/functional.h>
27 #include <deal.II/multigrid/mg_constrained_dofs.h>
40 DEAL_II_NAMESPACE_OPEN
59 template <
class VECTOR>
94 template <
class DOFINFO>
104 template <
class DOFINFO>
110 template<
class DOFINFO>
112 const DOFINFO &info2);
151 template <
class MATRIX>
198 template <
class DOFINFO>
204 template<
class DOFINFO>
211 template<
class DOFINFO>
213 const DOFINFO &info2);
220 const std::vector<types::global_dof_index> &i1,
221 const std::vector<types::global_dof_index> &i2);
254 template <
class MATRIX>
311 template <
class DOFINFO>
317 template<
class DOFINFO>
324 template<
class DOFINFO>
326 const DOFINFO &info2);
333 const std::vector<types::global_dof_index> &i1,
334 const std::vector<types::global_dof_index> &i2);
341 const std::vector<types::global_dof_index> &i1,
342 const std::vector<types::global_dof_index> &i2,
343 const unsigned int level);
351 const std::vector<types::global_dof_index> &i1,
352 const std::vector<types::global_dof_index> &i2,
360 const std::vector<types::global_dof_index> &i1,
361 const std::vector<types::global_dof_index> &i2,
370 const std::vector<types::global_dof_index> &i1,
371 const std::vector<types::global_dof_index> &i2,
380 const std::vector<types::global_dof_index> &i1,
381 const std::vector<types::global_dof_index> &i2,
437 template <
class MATRIX,
class VECTOR>
478 template <
class DOFINFO>
486 template<
class DOFINFO>
495 template<
class DOFINFO>
497 const DOFINFO &info2);
503 template <
class VECTOR>
510 template <
class VECTOR>
518 template <
class MATRIX>
524 template <
class VECTOR>
525 template <
class DOFINFO>
529 info.initialize_vectors(residuals.size());
533 template <
class VECTOR>
534 template <
class DOFINFO>
538 for (
unsigned int k=0; k<residuals.size(); ++k)
540 if (constraints == 0)
542 for (
unsigned int i=0; i<info.vector(k).block(0).size(); ++i)
543 (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
547 if (info.indices_by_block.size() == 0)
548 constraints->distribute_local_to_global(info.vector(k).block(0), info.indices, (*residuals(k)));
550 for (
unsigned int i=0; i != info.vector(k).n_blocks(); ++i)
551 constraints->distribute_local_to_global(info.vector(k).block(i), info.indices_by_block[i], (*residuals(k)));
556 template <
class VECTOR>
557 template <
class DOFINFO>
560 const DOFINFO &info2)
562 for (
unsigned int k=0; k<residuals.size(); ++k)
564 if (constraints == 0)
566 for (
unsigned int i=0; i<info1.vector(k).block(0).size(); ++i)
567 (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
568 for (
unsigned int i=0; i<info2.vector(k).block(0).size(); ++i)
569 (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
573 if (info1.indices_by_block.size() == 0 && info2.indices_by_block.size() == 0)
575 constraints->distribute_local_to_global
576 (info1.vector(k).block(0), info1.indices, (*residuals(k)));
577 constraints->distribute_local_to_global
578 (info2.vector(k).block(0), info2.indices, (*residuals(k)));
580 else if (info1.indices_by_block.size() != 0 && info2.indices_by_block.size() != 0)
582 for (
unsigned int i=0; i<info1.vector(k).n_blocks(); ++i)
584 constraints->distribute_local_to_global
585 (info1.vector(k).block(i), info1.indices_by_block[i], (*residuals(k)));
586 constraints->distribute_local_to_global
587 (info2.vector(k).block(i), info2.indices_by_block[i], (*residuals(k)));
601 template <
class MATRIX>
609 template <
class MATRIX>
617 template <
class MATRIX>
625 template <
class MATRIX>
631 template <
class MATRIX >
632 template <
class DOFINFO>
636 const unsigned int n = info.indices_by_block.size();
639 info.initialize_matrices(1, face);
642 info.initialize_matrices(n*n, face);
644 for (
unsigned int i=0; i<n; ++i)
645 for (
unsigned int j=0; j<n; ++j,++k)
647 info.matrix(k,
false).row = i;
648 info.matrix(k,
false).column = j;
651 info.matrix(k,
true).row = i;
652 info.matrix(k,
true).column = j;
660 template <
class MATRIX>
663 const std::vector<types::global_dof_index> &i1,
664 const std::vector<types::global_dof_index> &i2)
669 if (constraints == 0)
671 for (
unsigned int j=0; j<i1.size(); ++j)
672 for (
unsigned int k=0; k<i2.size(); ++k)
673 if (std::fabs(M(j,k)) >= threshold)
674 matrix->add(i1[j], i2[k], M(j,k));
677 constraints->distribute_local_to_global(M, i1, i2, *matrix);
681 template <
class MATRIX>
682 template <
class DOFINFO>
688 if (info.indices_by_block.size() == 0)
689 assemble(info.matrix(0,
false).matrix, info.indices, info.indices);
692 for (
unsigned int k=0; k<info.n_matrices(); ++k)
694 assemble(info.matrix(k,
false).matrix,
695 info.indices_by_block[info.matrix(k,
false).row],
696 info.indices_by_block[info.matrix(k,
false).column]);
702 template <
class MATRIX>
703 template <
class DOFINFO>
710 if (info1.indices_by_block.size() == 0 && info2.indices_by_block.size() == 0)
712 assemble(info1.matrix(0,
false).matrix, info1.indices, info1.indices);
713 assemble(info1.matrix(0,
true).matrix, info1.indices, info2.indices);
714 assemble(info2.matrix(0,
false).matrix, info2.indices, info2.indices);
715 assemble(info2.matrix(0,
true).matrix, info2.indices, info1.indices);
717 else if (info1.indices_by_block.size() != 0 && info2.indices_by_block.size() != 0)
718 for (
unsigned int k=0; k<info1.n_matrices(); ++k)
720 const unsigned int row = info1.matrix(k,
false).row;
721 const unsigned int column = info1.matrix(k,
false).column;
723 assemble(info1.matrix(k,
false).matrix,
724 info1.indices_by_block[row], info1.indices_by_block[column]);
725 assemble(info1.matrix(k,
true).matrix,
726 info1.indices_by_block[row], info2.indices_by_block[column]);
727 assemble(info2.matrix(k,
false).matrix,
728 info2.indices_by_block[row], info2.indices_by_block[column]);
729 assemble(info2.matrix(k,
true).matrix,
730 info2.indices_by_block[row], info1.indices_by_block[column]);
741 template <
class MATRIX>
749 template <
class MATRIX>
756 template <
class MATRIX>
760 mg_constrained_dofs = &c;
763 template <
class MATRIX>
769 template <
class MATRIX>
779 template <
class MATRIX>
785 interface_out = &out;
789 template <
class MATRIX >
790 template <
class DOFINFO>
794 const unsigned int n = info.indices_by_block.size();
797 info.initialize_matrices(1, face);
800 info.initialize_matrices(n*n, face);
802 for (
unsigned int i=0; i<n; ++i)
803 for (
unsigned int j=0; j<n; ++j,++k)
805 info.matrix(k,
false).row = i;
806 info.matrix(k,
false).column = j;
809 info.matrix(k,
true).row = i;
810 info.matrix(k,
true).column = j;
817 template <
class MATRIX>
822 const std::vector<types::global_dof_index> &i1,
823 const std::vector<types::global_dof_index> &i2)
828 if (mg_constrained_dofs == 0)
830 for (
unsigned int j=0; j<i1.size(); ++j)
831 for (
unsigned int k=0; k<i2.size(); ++k)
832 if (std::fabs(M(j,k)) >= threshold)
833 G.add(i1[j], i2[k], M(j,k));
837 for (
unsigned int j=0; j<i1.size(); ++j)
838 for (
unsigned int k=0; k<i2.size(); ++k)
839 if (std::fabs(M(j,k)) >= threshold)
841 if (!mg_constrained_dofs->continuity_across_refinement_edges())
842 G.add(i1[j], i2[k], M(j,k));
848 template <
class MATRIX>
853 const std::vector<types::global_dof_index> &i1,
854 const std::vector<types::global_dof_index> &i2,
855 const unsigned int level)
860 if (mg_constrained_dofs == 0)
862 for (
unsigned int j=0; j<i1.size(); ++j)
863 for (
unsigned int k=0; k<i2.size(); ++k)
864 if (std::fabs(M(j,k)) >= threshold)
865 G.add(i1[j], i2[k], M(j,k));
869 for (
unsigned int j=0; j<i1.size(); ++j)
870 for (
unsigned int k=0; k<i2.size(); ++k)
871 if (std::fabs(M(j,k)) >= threshold)
872 if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
873 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
875 if (mg_constrained_dofs->set_boundary_values())
879 if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
880 !mg_constrained_dofs->is_boundary_index(level, i2[k]))
882 (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
883 mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
885 G.add(i1[j], i2[k], M(j,k));
888 G.add(i1[j], i2[k], M(j,k));
894 template <
class MATRIX>
899 const std::vector<types::global_dof_index> &i1,
900 const std::vector<types::global_dof_index> &i2,
901 const unsigned int level)
906 if (mg_constrained_dofs == 0)
908 for (
unsigned int j=0; j<i1.size(); ++j)
909 for (
unsigned int k=0; k<i2.size(); ++k)
910 if (std::fabs(M(k,j)) >= threshold)
911 G.add(i1[j], i2[k], M(k,j));
915 for (
unsigned int j=0; j<i1.size(); ++j)
916 for (
unsigned int k=0; k<i2.size(); ++k)
917 if (std::fabs(M(k,j)) >= threshold)
918 if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
919 G.add(i1[j], i2[k], M(k,j));
923 template <
class MATRIX>
928 const std::vector<types::global_dof_index> &i1,
929 const std::vector<types::global_dof_index> &i2,
930 const unsigned int level)
935 if (mg_constrained_dofs == 0)
937 for (
unsigned int j=0; j<i1.size(); ++j)
938 for (
unsigned int k=0; k<i2.size(); ++k)
939 if (std::fabs(M(j,k)) >= threshold)
940 G.add(i1[j], i2[k], M(j,k));
944 for (
unsigned int j=0; j<i1.size(); ++j)
945 for (
unsigned int k=0; k<i2.size(); ++k)
946 if (std::fabs(M(j,k)) >= threshold)
947 if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
948 G.add(i1[j], i2[k], M(j,k));
952 template <
class MATRIX>
957 const std::vector<types::global_dof_index> &i1,
958 const std::vector<types::global_dof_index> &i2,
959 const unsigned int level)
965 for (
unsigned int j=0; j<i1.size(); ++j)
966 for (
unsigned int k=0; k<i2.size(); ++k)
967 if (std::fabs(M(j,k)) >= threshold)
979 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
980 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
982 if (mg_constrained_dofs->set_boundary_values())
984 if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
985 !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
987 (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
988 mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
990 G.add(i1[j], i2[k], M(j,k));
993 G.add(i1[j], i2[k], M(j,k));
998 template <
class MATRIX>
1003 const std::vector<types::global_dof_index> &i1,
1004 const std::vector<types::global_dof_index> &i2,
1005 const unsigned int level)
1011 for (
unsigned int j=0; j<i1.size(); ++j)
1012 for (
unsigned int k=0; k<i2.size(); ++k)
1013 if (std::fabs(M(k,j)) >= threshold)
1014 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1015 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1017 if (mg_constrained_dofs->set_boundary_values())
1019 if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
1020 !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
1022 (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
1023 mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
1025 G.add(i1[j], i2[k], M(k,j));
1028 G.add(i1[j], i2[k], M(k,j));
1033 template <
class MATRIX>
1034 template <
class DOFINFO>
1039 const unsigned int level = info.cell->level();
1041 if (info.indices_by_block.size() == 0)
1043 assemble((*matrix)[level], info.matrix(0,
false).matrix,
1044 info.indices, info.indices, level);
1045 if (mg_constrained_dofs != 0)
1047 assemble_in((*interface_in)[level], info.matrix(0,
false).matrix,
1048 info.indices, info.indices, level);
1049 assemble_out((*interface_out)[level],info.matrix(0,
false).matrix,
1050 info.indices, info.indices, level);
1054 for (
unsigned int k=0; k<info.n_matrices(); ++k)
1056 const unsigned int row = info.matrix(k,
false).row;
1057 const unsigned int column = info.matrix(k,
false).column;
1059 assemble((*matrix)[level], info.matrix(k,
false).matrix,
1060 info.indices_by_block[row], info.indices_by_block[column], level);
1062 if (mg_constrained_dofs != 0)
1064 assemble_in((*interface_in)[level], info.matrix(k,
false).matrix,
1065 info.indices_by_block[row], info.indices_by_block[column], level);
1066 assemble_out((*interface_out)[level],info.matrix(k,
false).matrix,
1067 info.indices_by_block[column], info.indices_by_block[row], level);
1073 template <
class MATRIX>
1074 template <
class DOFINFO>
1077 const DOFINFO &info2)
1081 const unsigned int level1 = info1.cell->level();
1082 const unsigned int level2 = info2.cell->level();
1084 if (info1.indices_by_block.size() == 0)
1086 if (level1 == level2)
1088 if (mg_constrained_dofs == 0)
1090 assemble((*matrix)[level1], info1.matrix(0,
false).matrix, info1.indices, info1.indices);
1091 assemble((*matrix)[level1], info1.matrix(0,
true).matrix, info1.indices, info2.indices);
1092 assemble((*matrix)[level1], info2.matrix(0,
false).matrix, info2.indices, info2.indices);
1093 assemble((*matrix)[level1], info2.matrix(0,
true).matrix, info2.indices, info1.indices);
1097 assemble((*matrix)[level1], info1.matrix(0,
false).matrix, info1.indices, info1.indices, level1);
1098 assemble((*matrix)[level1], info1.matrix(0,
true).matrix, info1.indices, info2.indices, level1);
1099 assemble((*matrix)[level1], info2.matrix(0,
false).matrix, info2.indices, info2.indices, level1);
1100 assemble((*matrix)[level1], info2.matrix(0,
true).matrix, info2.indices, info1.indices, level1);
1109 assemble((*matrix)[level1], info1.matrix(0,
false).matrix, info1.indices, info1.indices);
1112 assemble_up((*flux_up)[level1],info1.matrix(0,
true).matrix, info2.indices, info1.indices, level1);
1113 assemble_down((*flux_down)[level1], info2.matrix(0,
true).matrix, info2.indices, info1.indices, level1);
1118 for (
unsigned int k=0; k<info1.n_matrices(); ++k)
1120 const unsigned int row = info1.matrix(k,
false).row;
1121 const unsigned int column = info1.matrix(k,
false).column;
1123 if (level1 == level2)
1125 if (mg_constrained_dofs == 0)
1127 assemble((*matrix)[level1], info1.matrix(k,
false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
1128 assemble((*matrix)[level1], info1.matrix(k,
true).matrix, info1.indices_by_block[row], info2.indices_by_block[column]);
1129 assemble((*matrix)[level1], info2.matrix(k,
false).matrix, info2.indices_by_block[row], info2.indices_by_block[column]);
1130 assemble((*matrix)[level1], info2.matrix(k,
true).matrix, info2.indices_by_block[row], info1.indices_by_block[column]);
1134 assemble((*matrix)[level1], info1.matrix(k,
false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
1135 assemble((*matrix)[level1], info1.matrix(k,
true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
1136 assemble((*matrix)[level1], info2.matrix(k,
false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
1137 assemble((*matrix)[level1], info2.matrix(k,
true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
1146 assemble((*matrix)[level1], info1.matrix(k,
false).matrix, info1.indices_by_block[row], info1.indices_by_block[column]);
1149 assemble_up((*flux_up)[level1],info1.matrix(k,
true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
1150 assemble_down((*flux_down)[level1], info2.matrix(k,
true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
1158 template <
class MATRIX,
class VECTOR>
1165 template <
class MATRIX,
class VECTOR>
1171 data.
add(p,
"right hand side");
1177 template <
class MATRIX,
class VECTOR>
1186 template <
class MATRIX,
class VECTOR>
1187 template <
class DOFINFO>
1197 template <
class MATRIX,
class VECTOR>
1198 template <
class DOFINFO>
1207 template <
class MATRIX,
class VECTOR>
1208 template <
class DOFINFO>
1211 const DOFINFO &info2)
1219 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
void assemble(const DOFINFO &info)
#define AssertDimension(dim1, dim2)
void initialize_info(DOFINFO &info, bool face) const
void initialize(MATRIX &m)
void assemble(const DOFINFO &info)
void add(DATA &v, const std::string &name)
void initialize_info(DOFINFO &info, bool face) const
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > flux_up
::ExceptionBase & ExcMessage(std::string arg1)
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
MatrixSimple(double threshold=1.e-12)
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > interface_in
void initialize_local_blocks(const BlockIndices &)
void initialize_local_blocks(const BlockIndices &)
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > matrix
void initialize_info(DOFINFO &info, bool face) const
void initialize_fluxes(MGLevelObject< MATRIX > &flux_up, MGLevelObject< MATRIX > &flux_down)
void assemble_down(MATRIX &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
MGMatrixSimple(double threshold=1.e-12)
SmartPointer< const MGConstrainedDoFs, MGMatrixSimple< MATRIX > > mg_constrained_dofs
#define Assert(cond, exc)
SmartPointer< const ConstraintMatrix, MatrixSimple< MATRIX > > constraints
void assemble_out(MATRIX &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
SystemSimple(double threshold=1.e-12)
void initialize_local_blocks(const BlockIndices &)
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > flux_down
SmartPointer< MGLevelObject< MATRIX >, MGMatrixSimple< MATRIX > > interface_out
NamedData< SmartPointer< VECTOR, ResidualSimple< VECTOR > > > residuals
void initialize(MATRIX &m, VECTOR &rhs)
void assemble(const DOFINFO &info)
void initialize_info(DOFINFO &info, bool face) const
void assemble_up(MATRIX &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
void assemble_in(MATRIX &G, const FullMatrix< double > &M, const std::vector< types::global_dof_index > &i1, const std::vector< types::global_dof_index > &i2, const unsigned int level=numbers::invalid_unsigned_int)
SmartPointer< const ConstraintMatrix, ResidualSimple< VECTOR > > constraints
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcInternalError()
SmartPointer< MATRIX, MatrixSimple< MATRIX > > matrix
void assemble(const DOFINFO &info)
void initialize(MGLevelObject< MATRIX > &m)
void initialize_interfaces(MGLevelObject< MATRIX > &interface_in, MGLevelObject< MATRIX > &interface_out)