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> 59 template <
class VECTOR>
94 template <
class DOFINFO>
104 template <
class DOFINFO>
110 template<
class DOFINFO>
112 const DOFINFO &info2);
151 template <
class MATRIX>
164 void initialize(
MATRIX &m);
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);
349 void assemble_up(
MATRIX &G,
351 const std::vector<types::global_dof_index> &i1,
352 const std::vector<types::global_dof_index> &i2,
358 void assemble_down(
MATRIX &G,
360 const std::vector<types::global_dof_index> &i1,
361 const std::vector<types::global_dof_index> &i2,
368 void assemble_in(
MATRIX &G,
370 const std::vector<types::global_dof_index> &i1,
371 const std::vector<types::global_dof_index> &i2,
378 void assemble_out(
MATRIX &G,
380 const std::vector<types::global_dof_index> &i1,
381 const std::vector<types::global_dof_index> &i2,
437 template <
class MATRIX,
class VECTOR>
454 void initialize(
MATRIX &m, VECTOR &rhs);
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)
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)
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)
576 (info1.vector(k).block(0), info1.indices, (*
residuals(k)));
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)
585 (info1.vector(k).block(i), info1.indices_by_block[i], (*
residuals(k)));
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)
671 for (
unsigned int j=0; j<i1.size(); ++j)
672 for (
unsigned int k=0; k<i2.size(); ++k)
674 matrix->add(i1[j], i2[k], M(j,k));
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>
763 template <
class MATRIX>
769 template <
class MATRIX>
779 template <
class MATRIX>
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)
830 for (
unsigned int j=0; j<i1.size(); ++j)
831 for (
unsigned int k=0; k<i2.size(); ++k)
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)
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)
862 for (
unsigned int j=0; j<i1.size(); ++j)
863 for (
unsigned int k=0; k<i2.size(); ++k)
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)
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)
908 for (
unsigned int j=0; j<i1.size(); ++j)
909 for (
unsigned int k=0; k<i2.size(); ++k)
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)
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)
937 for (
unsigned int j=0; j<i1.size(); ++j)
938 for (
unsigned int k=0; k<i2.size(); ++k)
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)
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)
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)
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)
1044 info.indices, info.indices, level);
1048 info.indices, info.indices, level);
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;
1060 info.indices_by_block[row], info.indices_by_block[column], level);
1065 info.indices_by_block[row], info.indices_by_block[column], level);
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)
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);
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)
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)