18 #ifndef __deal2__mesh_worker_assembler_h
19 #define __deal2__mesh_worker_assembler_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>
28 #include <deal.II/multigrid/mg_constrained_dofs.h>
109 template <
class VECTOR>
137 template <
class DOFINFO>
145 template<
class DOFINFO>
152 template<
class DOFINFO>
154 const DOFINFO &info2);
162 const std::vector<types::global_dof_index> &dof);
212 template <
class MATRIX,
typename number =
double>
251 template <
class DOFINFO>
259 template<
class DOFINFO>
266 template<
class DOFINFO>
268 const DOFINFO &info2);
278 const unsigned int block_row,
279 const unsigned int block_col,
280 const std::vector<types::global_dof_index> &dof1,
281 const std::vector<types::global_dof_index> &dof2);
339 template <
class MATRIX,
typename number =
double>
409 template <
class DOFINFO>
417 template<
class DOFINFO>
424 template<
class DOFINFO>
426 const DOFINFO &info2);
436 const unsigned int block_row,
437 const unsigned int block_col,
438 const std::vector<types::global_dof_index> &dof1,
439 const std::vector<types::global_dof_index> &dof2,
440 const unsigned int level1,
441 const unsigned int level2,
442 bool transpose =
false);
451 const unsigned int block_row,
452 const unsigned int block_col,
453 const std::vector<types::global_dof_index> &dof1,
454 const std::vector<types::global_dof_index> &dof2,
455 const unsigned int level1,
456 const unsigned int level2);
465 const unsigned int block_row,
466 const unsigned int block_col,
467 const std::vector<types::global_dof_index> &dof1,
468 const std::vector<types::global_dof_index> &dof2,
469 const unsigned int level1,
470 const unsigned int level2);
479 const unsigned int block_row,
480 const unsigned int block_col,
481 const std::vector<types::global_dof_index> &dof1,
482 const std::vector<types::global_dof_index> &dof2,
483 const unsigned int level1,
484 const unsigned int level2);
493 const unsigned int block_row,
494 const unsigned int block_col,
495 const std::vector<types::global_dof_index> &dof1,
496 const std::vector<types::global_dof_index> &dof2,
497 const unsigned int level1,
498 const unsigned int level2);
507 const unsigned int block_row,
508 const unsigned int block_col,
509 const std::vector<types::global_dof_index> &dof1,
510 const std::vector<types::global_dof_index> &dof2,
511 const unsigned int level1,
512 const unsigned int level2);
575 template <
class VECTOR>
584 template <
class VECTOR>
593 template <
class VECTOR>
594 template <
class DOFINFO>
597 DOFINFO &info,
bool)
const
599 info.initialize_vectors(residuals.size());
602 template <
class VECTOR>
607 const std::vector<types::global_dof_index> &dof)
609 if (constraints == 0)
611 for (
unsigned int b=0; b<local.
n_blocks(); ++b)
612 for (
unsigned int j=0; j<local.
block(b).size(); ++j)
622 const unsigned int jcell = this->block_info->local().local_to_global(b, j);
623 global(dof[jcell]) += local.
block(b)(j);
627 constraints->distribute_local_to_global(local, dof, global);
631 template <
class VECTOR>
632 template <
class DOFINFO>
637 for (
unsigned int i=0; i<residuals.size(); ++i)
638 assemble(*residuals(i), info.vector(i), info.indices);
642 template <
class VECTOR>
643 template <
class DOFINFO>
646 const DOFINFO &info1,
647 const DOFINFO &info2)
649 for (
unsigned int i=0; i<residuals.size(); ++i)
651 assemble(*residuals(i), info1.vector(i), info1.indices);
652 assemble(*residuals(i), info2.vector(i), info2.indices);
659 template <
class MATRIX,
typename number>
668 template <
class MATRIX,
typename number>
680 template <
class MATRIX,
typename number>
690 template <
class MATRIX ,
typename number>
691 template <
class DOFINFO>
697 info.initialize_matrices(*matrices, face);
702 template <
class MATRIX,
typename number>
707 const unsigned int block_row,
708 const unsigned int block_col,
709 const std::vector<types::global_dof_index> &dof1,
710 const std::vector<types::global_dof_index> &dof2)
712 if (constraints == 0)
714 for (
unsigned int j=0; j<local.n_rows(); ++j)
715 for (
unsigned int k=0; k<local.n_cols(); ++k)
716 if (std::fabs(local(j,k)) >= threshold)
726 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
727 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
729 global.
add(dof1[jcell], dof2[kcell], local(j,k));
735 std::vector<types::global_dof_index> sliced_row_indices (bi.
block_size(block_row));
736 for (
unsigned int i=0; i<sliced_row_indices.size(); ++i)
737 sliced_row_indices[i] = dof1[bi.
block_start(block_row)+i];
739 std::vector<types::global_dof_index> sliced_col_indices (bi.
block_size(block_col));
740 for (
unsigned int i=0; i<sliced_col_indices.size(); ++i)
741 sliced_col_indices[i] = dof2[bi.
block_start(block_col)+i];
743 constraints->distribute_local_to_global(local,
744 sliced_row_indices, sliced_col_indices, global);
749 template <
class MATRIX,
typename number>
750 template <
class DOFINFO>
755 for (
unsigned int i=0; i<matrices->size(); ++i)
762 assemble(matrices->block(i), info.matrix(i,
false).matrix, row, col, info.indices, info.indices);
767 template <
class MATRIX,
typename number>
768 template <
class DOFINFO>
771 const DOFINFO &info1,
772 const DOFINFO &info2)
774 for (
unsigned int i=0; i<matrices->size(); ++i)
781 assemble(matrices->block(i), info1.matrix(i,
false).matrix, row, col, info1.indices, info1.indices);
782 assemble(matrices->block(i), info1.matrix(i,
true).matrix, row, col, info1.indices, info2.indices);
783 assemble(matrices->block(i), info2.matrix(i,
false).matrix, row, col, info2.indices, info2.indices);
784 assemble(matrices->block(i), info2.matrix(i,
true).matrix, row, col, info2.indices, info1.indices);
791 template <
class MATRIX,
typename number>
800 template <
class MATRIX,
typename number>
807 AssertDimension(block_info->local().size(), block_info->global().size());
812 template <
class MATRIX,
typename number>
817 mg_constrained_dofs = &mg_c;
821 template <
class MATRIX ,
typename number>
822 template <
class DOFINFO>
828 info.initialize_matrices(*matrices, face);
833 template <
class MATRIX,
typename number>
844 template <
class MATRIX,
typename number>
855 template <
class MATRIX,
typename number>
860 const unsigned int block_row,
861 const unsigned int block_col,
862 const std::vector<types::global_dof_index> &dof1,
863 const std::vector<types::global_dof_index> &dof2,
864 const unsigned int level1,
865 const unsigned int level2,
868 for (
unsigned int j=0; j<local.n_rows(); ++j)
869 for (
unsigned int k=0; k<local.n_cols(); ++k)
870 if (std::fabs(local(j,k)) >= threshold)
880 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
881 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
891 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
892 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
894 if (mg_constrained_dofs == 0)
897 global.add(kglobal, jglobal, local(j,k));
899 global.add(jglobal, kglobal, local(j,k));
903 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
904 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
906 if (mg_constrained_dofs->set_boundary_values())
908 if ((!mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
909 !mg_constrained_dofs->is_boundary_index(level2, kglobal))
911 (mg_constrained_dofs->is_boundary_index(level1, jglobal) &&
912 mg_constrained_dofs->is_boundary_index(level2, kglobal) &&
916 global.add(kglobal, jglobal, local(j,k));
918 global.add(jglobal, kglobal, local(j,k));
924 global.add(kglobal, jglobal, local(j,k));
926 global.add(jglobal, kglobal, local(j,k));
934 template <
class MATRIX,
typename number>
939 const unsigned int block_row,
940 const unsigned int block_col,
941 const std::vector<types::global_dof_index> &dof1,
942 const std::vector<types::global_dof_index> &dof2,
943 const unsigned int level1,
944 const unsigned int level2)
946 for (
unsigned int j=0; j<local.n_rows(); ++j)
947 for (
unsigned int k=0; k<local.n_cols(); ++k)
948 if (std::fabs(local(j,k)) >= threshold)
958 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
959 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
969 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
970 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
972 if (mg_constrained_dofs == 0)
973 global.add(jglobal, kglobal, local(j,k));
976 if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
977 !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
979 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
980 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
981 global.add(jglobal, kglobal, local(j,k));
987 template <
class MATRIX,
typename number>
992 const unsigned int block_row,
993 const unsigned int block_col,
994 const std::vector<types::global_dof_index> &dof1,
995 const std::vector<types::global_dof_index> &dof2,
996 const unsigned int level1,
997 const unsigned int level2)
999 for (
unsigned int j=0; j<local.n_rows(); ++j)
1000 for (
unsigned int k=0; k<local.n_cols(); ++k)
1001 if (std::fabs(local(j,k)) >= threshold)
1011 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1012 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1022 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1023 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1025 if (mg_constrained_dofs == 0)
1026 global.add(jglobal, kglobal, local(j,k));
1029 if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
1030 !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
1032 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1033 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1034 global.add(jglobal, kglobal, local(j,k));
1040 template <
class MATRIX,
typename number>
1045 const unsigned int block_row,
1046 const unsigned int block_col,
1047 const std::vector<types::global_dof_index> &dof1,
1048 const std::vector<types::global_dof_index> &dof2,
1049 const unsigned int level1,
1050 const unsigned int level2)
1052 for (
unsigned int j=0; j<local.n_rows(); ++j)
1053 for (
unsigned int k=0; k<local.n_cols(); ++k)
1054 if (std::fabs(local(k,j)) >= threshold)
1064 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1065 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1075 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1076 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1078 if (mg_constrained_dofs == 0)
1079 global.add(jglobal, kglobal, local(k,j));
1082 if (!mg_constrained_dofs->non_refinement_edge_index(level1, jglobal) &&
1083 !mg_constrained_dofs->non_refinement_edge_index(level2, kglobal))
1085 if (!mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1086 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1087 global.add(jglobal, kglobal, local(k,j));
1093 template <
class MATRIX,
typename number>
1098 const unsigned int block_row,
1099 const unsigned int block_col,
1100 const std::vector<types::global_dof_index> &dof1,
1101 const std::vector<types::global_dof_index> &dof2,
1102 const unsigned int level1,
1103 const unsigned int level2)
1108 for (
unsigned int j=0; j<local.n_rows(); ++j)
1109 for (
unsigned int k=0; k<local.n_cols(); ++k)
1110 if (std::fabs(local(j,k)) >= threshold)
1120 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1121 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1131 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1132 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1134 if (mg_constrained_dofs == 0)
1135 global.add(jglobal, kglobal, local(j,k));
1138 if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1139 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1141 if (mg_constrained_dofs->set_boundary_values())
1143 if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1144 !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
1146 (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1147 mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
1148 jglobal == kglobal))
1149 global.add(jglobal, kglobal, local(j,k));
1152 global.add(jglobal, kglobal, local(j,k));
1158 template <
class MATRIX,
typename number>
1163 const unsigned int block_row,
1164 const unsigned int block_col,
1165 const std::vector<types::global_dof_index> &dof1,
1166 const std::vector<types::global_dof_index> &dof2,
1167 const unsigned int level1,
1168 const unsigned int level2)
1173 for (
unsigned int j=0; j<local.n_rows(); ++j)
1174 for (
unsigned int k=0; k<local.n_cols(); ++k)
1175 if (std::fabs(local(k,j)) >= threshold)
1185 const unsigned int jcell = this->block_info->local().local_to_global(block_row, j);
1186 const unsigned int kcell = this->block_info->local().local_to_global(block_col, k);
1196 const unsigned int jglobal = this->block_info->level(level1).global_to_local(dof1[jcell]).second;
1197 const unsigned int kglobal = this->block_info->level(level2).global_to_local(dof2[kcell]).second;
1199 if (mg_constrained_dofs == 0)
1200 global.add(jglobal, kglobal, local(k,j));
1203 if (mg_constrained_dofs->at_refinement_edge(level1, jglobal) &&
1204 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1206 if (mg_constrained_dofs->set_boundary_values())
1208 if ((!mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1209 !mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal))
1211 (mg_constrained_dofs->at_refinement_edge_boundary(level1, jglobal) &&
1212 mg_constrained_dofs->at_refinement_edge_boundary(level2, kglobal) &&
1213 jglobal == kglobal))
1214 global.add(jglobal, kglobal, local(k,j));
1217 global.add(jglobal, kglobal, local(k,j));
1224 template <
class MATRIX,
typename number>
1225 template <
class DOFINFO>
1229 const unsigned int level = info.cell->level();
1231 for (
unsigned int i=0; i<matrices->size(); ++i)
1235 const unsigned int row = matrices->block(i)[level].row;
1236 const unsigned int col = matrices->block(i)[level].column;
1238 assemble(matrices->block(i)[level].matrix, info.matrix(i,
false).matrix, row, col,
1239 info.indices, info.indices, level, level);
1240 if (mg_constrained_dofs != 0)
1242 if (interface_in != 0)
1243 assemble_in(interface_in->block(i)[level], info.matrix(i,
false).matrix, row, col,
1244 info.indices, info.indices, level, level);
1245 if (interface_out != 0)
1246 assemble_in(interface_out->block(i)[level], info.matrix(i,
false).matrix, row, col,
1247 info.indices, info.indices, level, level);
1249 assemble_in(matrices->block_in(i)[level], info.matrix(i,
false).matrix, row, col,
1250 info.indices, info.indices, level, level);
1251 assemble_out(matrices->block_out(i)[level], info.matrix(i,
false).matrix, row, col,
1252 info.indices, info.indices, level, level);
1258 template <
class MATRIX,
typename number>
1259 template <
class DOFINFO>
1262 const DOFINFO &info1,
1263 const DOFINFO &info2)
1265 const unsigned int level1 = info1.cell->level();
1266 const unsigned int level2 = info2.cell->level();
1268 for (
unsigned int i=0; i<matrices->size(); ++i)
1274 const unsigned int row = o[level1].row;
1275 const unsigned int col = o[level1].column;
1277 if (level1 == level2)
1279 if (mg_constrained_dofs == 0)
1281 assemble(o[level1].matrix, info1.matrix(i,
false).matrix, row, col,
1282 info1.indices, info1.indices, level1, level1);
1283 assemble(o[level1].matrix, info1.matrix(i,
true).matrix, row, col,
1284 info1.indices, info2.indices, level1, level2);
1285 assemble(o[level1].matrix, info2.matrix(i,
false).matrix, row, col,
1286 info2.indices, info2.indices, level2, level2);
1287 assemble(o[level1].matrix, info2.matrix(i,
true).matrix, row, col,
1288 info2.indices, info1.indices, level2, level1);
1292 assemble_fluxes(o[level1].matrix, info1.matrix(i,
false).matrix, row, col,
1293 info1.indices, info1.indices, level1, level1);
1294 assemble_fluxes(o[level1].matrix, info1.matrix(i,
true).matrix, row, col,
1295 info1.indices, info2.indices, level1, level2);
1296 assemble_fluxes(o[level1].matrix, info2.matrix(i,
false).matrix, row, col,
1297 info2.indices, info2.indices, level2, level2);
1298 assemble_fluxes(o[level1].matrix, info2.matrix(i,
true).matrix, row, col,
1299 info2.indices, info1.indices, level2, level1);
1305 if (flux_up->size() != 0)
1310 assemble_fluxes(o[level1].matrix, info1.matrix(i,
false).matrix, row, col,
1311 info1.indices, info1.indices, level1, level1);
1312 assemble_up(flux_up->block(i)[level1].matrix, info1.matrix(i,
true).matrix, row, col,
1313 info1.indices, info2.indices, level1, level2);
1314 assemble_down(flux_down->block(i)[level1].matrix, info2.matrix(i,
true).matrix, row, col,
1315 info2.indices, info1.indices, level2, level1);
1323 DEAL_II_NAMESPACE_CLOSE
void assemble_in(MATRIX &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)
SmartPointer< const BlockInfo, MatrixLocalBlocksToGlobalBlocks< MATRIX, number > > block_info
#define AssertDimension(dim1, dim2)
MeshWorker::Assember objects for systems without block structure.
MatrixPtrVectorPtr interface_out
void initialize_info(DOFINFO &info, bool face) const
SmartPointer< MatrixBlockVector< MATRIX >, MatrixLocalBlocksToGlobalBlocks< MATRIX, number > > matrices
MatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
void add(const size_type i, const size_type j, const typename MATRIX::value_type value)
void assemble(const DOFINFO &info)
SmartPointer< const ConstraintMatrix, ResidualLocalBlocksToGlobalBlocks< VECTOR > > constraints
unsigned int n_blocks() const
MatrixPtrVectorPtr interface_in
void initialize_interfaces(MatrixPtrVector &interface_in, MatrixPtrVector &interface_out)
MatrixPtrVectorPtr flux_down
SmartPointer< const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks< MATRIX, number > > block_info
size_type block_start(const unsigned int i) const
void initialize_info(DOFINFO &info, bool face) const
void assemble(const DOFINFO &info)
MatrixPtrVectorPtr matrices
unsigned int global_dof_index
#define Assert(cond, exc)
void assemble_down(MATRIX &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)
void assemble_fluxes(MATRIX &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)
SmartPointer< const ConstraintMatrix, MatrixLocalBlocksToGlobalBlocks< MATRIX, number > > constraints
void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down)
size_type block_size(const unsigned int i) const
void initialize(const BlockInfo *block_info, MatrixPtrVector &matrices)
NamedData< SmartPointer< VECTOR, ResidualLocalBlocksToGlobalBlocks< VECTOR > > > residuals
void initialize_info(DOFINFO &info, bool face) const
void assemble_out(MATRIX &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)
::ExceptionBase & ExcNotImplemented()
A small class collecting the different BlockIndices involved in global, multilevel and local computat...
void assemble(const DOFINFO &info)
MatrixPtrVectorPtr flux_up
void initialize(const BlockInfo *block_info, MatrixBlockVector< MATRIX > &matrices)
MGMatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
BlockType & block(const unsigned int i)
void initialize(const BlockInfo *block_info, NamedData< VECTOR * > &residuals)
SmartPointer< const MGConstrainedDoFs, MGMatrixLocalBlocksToGlobalBlocks< MATRIX, number > > mg_constrained_dofs
SmartPointer< const BlockInfo, ResidualLocalBlocksToGlobalBlocks< VECTOR > > block_info
void assemble_up(MATRIX &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)