17 #ifndef dealii_mesh_worker_assembler_h 18 #define dealii_mesh_worker_assembler_h 20 #include <deal.II/algorithms/any_data.h> 22 #include <deal.II/base/mg_level_object.h> 23 #include <deal.II/base/smartpointer.h> 25 #include <deal.II/lac/block_vector.h> 27 #include <deal.II/meshworker/dof_info.h> 28 #include <deal.II/meshworker/functional.h> 29 #include <deal.II/meshworker/simple.h> 31 #include <deal.II/multigrid/mg_constrained_dofs.h> 34 DEAL_II_NAMESPACE_OPEN
110 template <
typename VectorType>
134 template <
class DOFINFO>
142 template <
class DOFINFO>
149 template <
class DOFINFO>
151 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
160 const std::vector<types::global_dof_index> &dof);
210 template <
typename MatrixType,
typename number =
double>
242 template <
class DOFINFO>
250 template <
class DOFINFO>
257 template <
class DOFINFO>
259 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
268 const unsigned int block_row,
269 const unsigned int block_col,
270 const std::vector<types::global_dof_index> &dof1,
271 const std::vector<types::global_dof_index> &dof2);
326 template <
typename MatrixType,
typename number =
double>
377 template <
class DOFINFO>
385 template <
class DOFINFO>
392 template <
class DOFINFO>
394 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
403 const unsigned int block_row,
404 const unsigned int block_col,
405 const std::vector<types::global_dof_index> &dof1,
406 const std::vector<types::global_dof_index> &dof2,
407 const unsigned int level1,
408 const unsigned int level2,
417 const unsigned int block_row,
418 const unsigned int block_col,
419 const std::vector<types::global_dof_index> &dof1,
420 const std::vector<types::global_dof_index> &dof2,
421 const unsigned int level1,
422 const unsigned int level2);
430 const unsigned int block_row,
431 const unsigned int block_col,
432 const std::vector<types::global_dof_index> &dof1,
433 const std::vector<types::global_dof_index> &dof2,
434 const unsigned int level1,
435 const unsigned int level2);
443 const unsigned int block_row,
444 const unsigned int block_col,
445 const std::vector<types::global_dof_index> &dof1,
446 const std::vector<types::global_dof_index> &dof2,
447 const unsigned int level1,
448 const unsigned int level2);
456 const unsigned int block_row,
457 const unsigned int block_col,
458 const std::vector<types::global_dof_index> &dof1,
459 const std::vector<types::global_dof_index> &dof2,
460 const unsigned int level1,
461 const unsigned int level2);
469 const unsigned int block_row,
470 const unsigned int block_col,
471 const std::vector<types::global_dof_index> &dof1,
472 const std::vector<types::global_dof_index> &dof2,
473 const unsigned int level1,
474 const unsigned int level2);
530 template <
typename VectorType>
540 template <
typename VectorType>
549 template <
typename VectorType>
550 template <
class DOFINFO>
556 info.initialize_vectors(residuals.size());
559 template <
typename VectorType>
564 const std::vector<types::global_dof_index> &dof)
566 if (constraints == 0)
568 for (
unsigned int b = 0; b < local.
n_blocks(); ++b)
569 for (
unsigned int j = 0; j < local.
block(b).size(); ++j)
579 const unsigned int jcell =
580 this->block_info->local().local_to_global(b, j);
581 global(dof[jcell]) += local.
block(b)(j);
585 constraints->distribute_local_to_global(local, dof, global);
589 template <
typename VectorType>
590 template <
class DOFINFO>
594 for (
unsigned int i = 0; i < residuals.size(); ++i)
595 assemble(*(residuals.entry<VectorType>(i)),
601 template <
typename VectorType>
602 template <
class DOFINFO>
605 const DOFINFO &info1,
606 const DOFINFO &info2)
608 for (
unsigned int i = 0; i < residuals.size(); ++i)
610 assemble(*(residuals.entry<VectorType>(i)),
613 assemble(*(residuals.entry<VectorType>(i)),
622 template <
typename MatrixType,
typename number>
625 : threshold(threshold)
629 template <
typename MatrixType,
typename number>
641 template <
typename MatrixType,
typename number>
651 template <
typename MatrixType,
typename number>
652 template <
class DOFINFO>
658 info.initialize_matrices(*matrices, face);
663 template <
typename MatrixType,
typename number>
668 const unsigned int block_row,
669 const unsigned int block_col,
670 const std::vector<types::global_dof_index> &dof1,
671 const std::vector<types::global_dof_index> &dof2)
673 if (constraints ==
nullptr)
675 for (
unsigned int j = 0; j < local.n_rows(); ++j)
676 for (
unsigned int k = 0; k < local.n_cols(); ++k)
677 if (std::fabs(local(j, k)) >= threshold)
687 const unsigned int jcell =
688 this->block_info->local().local_to_global(block_row, j);
689 const unsigned int kcell =
690 this->block_info->local().local_to_global(block_col, k);
692 global.
add(dof1[jcell], dof2[kcell], local(j, k));
698 std::vector<types::global_dof_index> sliced_row_indices(
700 for (
unsigned int i = 0; i < sliced_row_indices.size(); ++i)
701 sliced_row_indices[i] = dof1[bi.
block_start(block_row) + i];
703 std::vector<types::global_dof_index> sliced_col_indices(
705 for (
unsigned int i = 0; i < sliced_col_indices.size(); ++i)
706 sliced_col_indices[i] = dof2[bi.
block_start(block_col) + i];
708 constraints->distribute_local_to_global(local,
716 template <
typename MatrixType,
typename number>
717 template <
class DOFINFO>
722 for (
unsigned int i = 0; i < matrices->size(); ++i)
729 assemble(matrices->block(i),
730 info.matrix(i,
false).matrix,
739 template <
typename MatrixType,
typename number>
740 template <
class DOFINFO>
743 const DOFINFO &info1,
744 const DOFINFO &info2)
746 for (
unsigned int i = 0; i < matrices->size(); ++i)
753 assemble(matrices->block(i),
754 info1.matrix(i,
false).matrix,
759 assemble(matrices->block(i),
760 info1.matrix(i,
true).matrix,
765 assemble(matrices->block(i),
766 info2.matrix(i,
false).matrix,
771 assemble(matrices->block(i),
772 info2.matrix(i,
true).matrix,
783 template <
typename MatrixType,
typename number>
786 : threshold(threshold)
790 template <
typename MatrixType,
typename number>
797 AssertDimension(block_info->local().size(), block_info->global().size());
802 template <
typename MatrixType,
typename number>
807 mg_constrained_dofs = &mg_c;
811 template <
typename MatrixType,
typename number>
812 template <
class DOFINFO>
818 info.initialize_matrices(*matrices, face);
823 template <
typename MatrixType,
typename number>
834 template <
typename MatrixType,
typename number>
844 template <
typename MatrixType,
typename number>
849 const unsigned int block_row,
850 const unsigned int block_col,
851 const std::vector<types::global_dof_index> &dof1,
852 const std::vector<types::global_dof_index> &dof2,
853 const unsigned int level1,
854 const unsigned int level2,
857 for (
unsigned int j = 0; j < local.n_rows(); ++j)
858 for (
unsigned int k = 0; k < local.n_cols(); ++k)
859 if (std::fabs(local(j, k)) >= threshold)
869 const unsigned int jcell =
870 this->block_info->local().local_to_global(block_row, j);
871 const unsigned int kcell =
872 this->block_info->local().local_to_global(block_col, k);
882 const unsigned int jglobal = this->block_info->level(level1)
883 .global_to_local(dof1[jcell])
885 const unsigned int kglobal = this->block_info->level(level2)
886 .global_to_local(dof2[kcell])
889 if (mg_constrained_dofs == 0)
892 global.add(kglobal, jglobal, local(j, k));
894 global.add(jglobal, kglobal, local(j, k));
898 if (!mg_constrained_dofs->at_refinement_edge(level1,
900 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
902 if (mg_constrained_dofs->set_boundary_values())
904 if ((!mg_constrained_dofs->is_boundary_index(
906 !mg_constrained_dofs->is_boundary_index(
908 (mg_constrained_dofs->is_boundary_index(
910 mg_constrained_dofs->is_boundary_index(
915 global.add(kglobal, jglobal, local(j, k));
917 global.add(jglobal, kglobal, local(j, k));
923 global.add(kglobal, jglobal, local(j, k));
925 global.add(jglobal, kglobal, local(j, k));
933 template <
typename MatrixType,
typename number>
938 const unsigned int block_row,
939 const unsigned int block_col,
940 const std::vector<types::global_dof_index> &dof1,
941 const std::vector<types::global_dof_index> &dof2,
942 const unsigned int level1,
943 const unsigned int level2)
945 for (
unsigned int j = 0; j < local.n_rows(); ++j)
946 for (
unsigned int k = 0; k < local.n_cols(); ++k)
947 if (std::fabs(local(j, k)) >= threshold)
957 const unsigned int jcell =
958 this->block_info->local().local_to_global(block_row, j);
959 const unsigned int kcell =
960 this->block_info->local().local_to_global(block_col, k);
970 const unsigned int jglobal = this->block_info->level(level1)
971 .global_to_local(dof1[jcell])
973 const unsigned int kglobal = this->block_info->level(level2)
974 .global_to_local(dof2[kcell])
977 if (mg_constrained_dofs == 0)
978 global.add(jglobal, kglobal, local(j, k));
981 if (!mg_constrained_dofs->non_refinement_edge_index(
983 !mg_constrained_dofs->non_refinement_edge_index(level2,
986 if (!mg_constrained_dofs->at_refinement_edge(level1,
988 !mg_constrained_dofs->at_refinement_edge(level2,
990 global.add(jglobal, kglobal, local(j, k));
996 template <
typename MatrixType,
typename number>
1001 const unsigned int block_row,
1002 const unsigned int block_col,
1003 const std::vector<types::global_dof_index> &dof1,
1004 const std::vector<types::global_dof_index> &dof2,
1005 const unsigned int level1,
1006 const unsigned int level2)
1008 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1009 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1010 if (std::fabs(local(j, k)) >= threshold)
1020 const unsigned int jcell =
1021 this->block_info->local().local_to_global(block_row, j);
1022 const unsigned int kcell =
1023 this->block_info->local().local_to_global(block_col, k);
1033 const unsigned int jglobal = this->block_info->level(level1)
1034 .global_to_local(dof1[jcell])
1036 const unsigned int kglobal = this->block_info->level(level2)
1037 .global_to_local(dof2[kcell])
1040 if (mg_constrained_dofs == 0)
1041 global.add(jglobal, kglobal, local(j, k));
1044 if (!mg_constrained_dofs->non_refinement_edge_index(
1046 !mg_constrained_dofs->non_refinement_edge_index(level2,
1049 if (!mg_constrained_dofs->at_refinement_edge(level1,
1051 !mg_constrained_dofs->at_refinement_edge(level2,
1053 global.add(jglobal, kglobal, local(j, k));
1059 template <
typename MatrixType,
typename number>
1062 MatrixType & global,
1064 const unsigned int block_row,
1065 const unsigned int block_col,
1066 const std::vector<types::global_dof_index> &dof1,
1067 const std::vector<types::global_dof_index> &dof2,
1068 const unsigned int level1,
1069 const unsigned int level2)
1071 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1072 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1073 if (std::fabs(local(k, j)) >= threshold)
1083 const unsigned int jcell =
1084 this->block_info->local().local_to_global(block_row, j);
1085 const unsigned int kcell =
1086 this->block_info->local().local_to_global(block_col, k);
1096 const unsigned int jglobal = this->block_info->level(level1)
1097 .global_to_local(dof1[jcell])
1099 const unsigned int kglobal = this->block_info->level(level2)
1100 .global_to_local(dof2[kcell])
1103 if (mg_constrained_dofs == 0)
1104 global.add(jglobal, kglobal, local(k, j));
1107 if (!mg_constrained_dofs->non_refinement_edge_index(
1109 !mg_constrained_dofs->non_refinement_edge_index(level2,
1112 if (!mg_constrained_dofs->at_refinement_edge(level1,
1114 !mg_constrained_dofs->at_refinement_edge(level2,
1116 global.add(jglobal, kglobal, local(k, j));
1122 template <
typename MatrixType,
typename number>
1125 MatrixType & global,
1127 const unsigned int block_row,
1128 const unsigned int block_col,
1129 const std::vector<types::global_dof_index> &dof1,
1130 const std::vector<types::global_dof_index> &dof2,
1131 const unsigned int level1,
1132 const unsigned int level2)
1137 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1138 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1139 if (std::fabs(local(j, k)) >= threshold)
1149 const unsigned int jcell =
1150 this->block_info->local().local_to_global(block_row, j);
1151 const unsigned int kcell =
1152 this->block_info->local().local_to_global(block_col, k);
1162 const unsigned int jglobal = this->block_info->level(level1)
1163 .global_to_local(dof1[jcell])
1165 const unsigned int kglobal = this->block_info->level(level2)
1166 .global_to_local(dof2[kcell])
1169 if (mg_constrained_dofs == 0)
1170 global.add(jglobal, kglobal, local(j, k));
1173 if (mg_constrained_dofs->at_refinement_edge(level1,
1175 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1177 if (mg_constrained_dofs->set_boundary_values())
1179 if ((!mg_constrained_dofs->is_boundary_index(
1181 !mg_constrained_dofs->is_boundary_index(
1182 level2, kglobal)) ||
1183 (mg_constrained_dofs->is_boundary_index(
1185 mg_constrained_dofs->is_boundary_index(
1187 jglobal == kglobal))
1188 global.add(jglobal, kglobal, local(j, k));
1191 global.add(jglobal, kglobal, local(j, k));
1197 template <
typename MatrixType,
typename number>
1200 MatrixType & global,
1202 const unsigned int block_row,
1203 const unsigned int block_col,
1204 const std::vector<types::global_dof_index> &dof1,
1205 const std::vector<types::global_dof_index> &dof2,
1206 const unsigned int level1,
1207 const unsigned int level2)
1212 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1213 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1214 if (std::fabs(local(k, j)) >= threshold)
1224 const unsigned int jcell =
1225 this->block_info->local().local_to_global(block_row, j);
1226 const unsigned int kcell =
1227 this->block_info->local().local_to_global(block_col, k);
1237 const unsigned int jglobal = this->block_info->level(level1)
1238 .global_to_local(dof1[jcell])
1240 const unsigned int kglobal = this->block_info->level(level2)
1241 .global_to_local(dof2[kcell])
1244 if (mg_constrained_dofs == 0)
1245 global.add(jglobal, kglobal, local(k, j));
1248 if (mg_constrained_dofs->at_refinement_edge(level1,
1250 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1252 if (mg_constrained_dofs->set_boundary_values())
1254 if ((!mg_constrained_dofs->is_boundary_index(
1256 !mg_constrained_dofs->is_boundary_index(
1257 level2, kglobal)) ||
1258 (mg_constrained_dofs->is_boundary_index(
1260 mg_constrained_dofs->is_boundary_index(
1262 jglobal == kglobal))
1263 global.add(jglobal, kglobal, local(k, j));
1266 global.add(jglobal, kglobal, local(k, j));
1273 template <
typename MatrixType,
typename number>
1274 template <
class DOFINFO>
1277 const DOFINFO &info)
1279 const unsigned int level = info.cell->level();
1281 for (
unsigned int i = 0; i < matrices->size(); ++i)
1285 const unsigned int row = matrices->block(i)[level].row;
1286 const unsigned int col = matrices->block(i)[level].column;
1288 assemble(matrices->block(i)[level].matrix,
1289 info.matrix(i,
false).matrix,
1296 if (mg_constrained_dofs != 0)
1298 if (interface_in != 0)
1299 assemble_in(interface_in->block(i)[level],
1300 info.matrix(i,
false).matrix,
1307 if (interface_out != 0)
1308 assemble_in(interface_out->block(i)[level],
1309 info.matrix(i,
false).matrix,
1317 assemble_in(matrices->block_in(i)[level],
1318 info.matrix(i,
false).matrix,
1325 assemble_out(matrices->block_out(i)[level],
1326 info.matrix(i,
false).matrix,
1338 template <
typename MatrixType,
typename number>
1339 template <
class DOFINFO>
1342 const DOFINFO &info1,
1343 const DOFINFO &info2)
1345 const unsigned int level1 = info1.cell->level();
1346 const unsigned int level2 = info2.cell->level();
1348 for (
unsigned int i = 0; i < matrices->size(); ++i)
1354 const unsigned int row = o[level1].row;
1355 const unsigned int col = o[level1].column;
1357 if (level1 == level2)
1359 if (mg_constrained_dofs == 0)
1361 assemble(o[level1].matrix,
1362 info1.matrix(i,
false).matrix,
1369 assemble(o[level1].matrix,
1370 info1.matrix(i,
true).matrix,
1377 assemble(o[level1].matrix,
1378 info2.matrix(i,
false).matrix,
1385 assemble(o[level1].matrix,
1386 info2.matrix(i,
true).matrix,
1396 assemble_fluxes(o[level1].matrix,
1397 info1.matrix(i,
false).matrix,
1404 assemble_fluxes(o[level1].matrix,
1405 info1.matrix(i,
true).matrix,
1412 assemble_fluxes(o[level1].matrix,
1413 info2.matrix(i,
false).matrix,
1420 assemble_fluxes(o[level1].matrix,
1421 info2.matrix(i,
true).matrix,
1433 if (flux_up->size() != 0)
1438 assemble_fluxes(o[level1].matrix,
1439 info1.matrix(i,
false).matrix,
1446 assemble_up(flux_up->block(i)[level1].matrix,
1447 info1.matrix(i,
true).matrix,
1454 assemble_down(flux_down->block(i)[level1].matrix,
1455 info2.matrix(i,
true).matrix,
1469 DEAL_II_NAMESPACE_CLOSE
SmartPointer< const BlockInfo, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
MatrixPtrVectorPtr interface_out
SmartPointer< const AffineConstraints< typename VectorType::value_type >, ResidualLocalBlocksToGlobalBlocks< VectorType > > constraints
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)
#define AssertDimension(dim1, dim2)
void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down)
void initialize(const BlockInfo *block_info, AnyData &residuals)
SmartPointer< const MGConstrainedDoFs, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > mg_constrained_dofs
void initialize_info(DOFINFO &info, bool face) const
SmartPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > constraints
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
MGMatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
size_type block_size(const unsigned int i) const
void initialize_info(DOFINFO &info, bool face) const
MatrixPtrVectorPtr flux_up
#define Assert(cond, exc)
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)
SmartPointer< const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
void assemble(const DOFINFO &info)
MatrixPtrVectorPtr flux_down
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
void initialize(const BlockInfo *block_info, MatrixBlockVector< MatrixType > &matrices)
SmartPointer< const BlockInfo, ResidualLocalBlocksToGlobalBlocks< VectorType > > block_info
unsigned int global_dof_index
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)
unsigned int n_blocks() const
void initialize_interfaces(MatrixPtrVector &interface_in, MatrixPtrVector &interface_out)
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)
void assemble(const DOFINFO &info)
size_type block_start(const unsigned int i) const
A small class collecting the different BlockIndices involved in global, multilevel and local computat...
void assemble(const DOFINFO &info)
static ::ExceptionBase & ExcNotImplemented()
MatrixPtrVectorPtr interface_in
SmartPointer< MatrixBlockVector< MatrixType >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > matrices
BlockType & block(const unsigned int i)
void initialize_info(DOFINFO &info, bool face) const
MatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
MatrixPtrVectorPtr matrices
void initialize(const BlockInfo *block_info, MatrixPtrVector &matrices)
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)