16#ifndef dealii_mesh_worker_assembler_h
17#define dealii_mesh_worker_assembler_h
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);
209 template <
typename MatrixType,
typename number =
double>
241 template <
class DOFINFO>
249 template <
class DOFINFO>
256 template <
class DOFINFO>
258 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
267 const unsigned int block_row,
268 const unsigned int block_col,
269 const std::vector<types::global_dof_index> &dof1,
270 const std::vector<types::global_dof_index> &dof2);
324 template <
typename MatrixType,
typename number =
double>
375 template <
class DOFINFO>
383 template <
class DOFINFO>
390 template <
class DOFINFO>
392 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
401 const unsigned int block_row,
402 const unsigned int block_col,
403 const std::vector<types::global_dof_index> &dof1,
404 const std::vector<types::global_dof_index> &dof2,
405 const unsigned int level1,
406 const unsigned int level2,
415 const unsigned int block_row,
416 const unsigned int block_col,
417 const std::vector<types::global_dof_index> &dof1,
418 const std::vector<types::global_dof_index> &dof2,
419 const unsigned int level1,
420 const unsigned int level2);
428 const unsigned int block_row,
429 const unsigned int block_col,
430 const std::vector<types::global_dof_index> &dof1,
431 const std::vector<types::global_dof_index> &dof2,
432 const unsigned int level1,
433 const unsigned int level2);
441 const unsigned int block_row,
442 const unsigned int block_col,
443 const std::vector<types::global_dof_index> &dof1,
444 const std::vector<types::global_dof_index> &dof2,
445 const unsigned int level1,
446 const unsigned int level2);
454 const unsigned int block_row,
455 const unsigned int block_col,
456 const std::vector<types::global_dof_index> &dof1,
457 const std::vector<types::global_dof_index> &dof2,
458 const unsigned int level1,
459 const unsigned int level2);
467 const unsigned int block_row,
468 const unsigned int block_col,
469 const std::vector<types::global_dof_index> &dof1,
470 const std::vector<types::global_dof_index> &dof2,
471 const unsigned int level1,
472 const unsigned int level2);
528 template <
typename VectorType>
538 template <
typename VectorType>
547 template <
typename VectorType>
548 template <
class DOFINFO>
554 info.initialize_vectors(residuals.size());
557 template <
typename VectorType>
562 const std::vector<types::global_dof_index> &dof)
564 if (constraints == 0)
566 for (
unsigned int b = 0; b < local.
n_blocks(); ++b)
567 for (
unsigned int j = 0; j < local.
block(b).size(); ++j)
577 const unsigned int jcell =
578 this->block_info->local().local_to_global(b, j);
579 global(dof[jcell]) += local.
block(b)(j);
583 constraints->distribute_local_to_global(local, dof, global);
587 template <
typename VectorType>
588 template <
class DOFINFO>
592 for (
unsigned int i = 0; i < residuals.size(); ++i)
593 assemble(*(residuals.entry<VectorType>(i)),
599 template <
typename VectorType>
600 template <
class DOFINFO>
603 const DOFINFO &info1,
604 const DOFINFO &info2)
606 for (
unsigned int i = 0; i < residuals.size(); ++i)
608 assemble(*(residuals.entry<VectorType>(i)),
611 assemble(*(residuals.entry<VectorType>(i)),
620 template <
typename MatrixType,
typename number>
623 : threshold(threshold)
627 template <
typename MatrixType,
typename number>
639 template <
typename MatrixType,
typename number>
649 template <
typename MatrixType,
typename number>
650 template <
class DOFINFO>
656 info.initialize_matrices(*matrices, face);
661 template <
typename MatrixType,
typename number>
666 const unsigned int block_row,
667 const unsigned int block_col,
668 const std::vector<types::global_dof_index> &dof1,
669 const std::vector<types::global_dof_index> &dof2)
671 if (constraints ==
nullptr)
673 for (
unsigned int j = 0; j < local.n_rows(); ++j)
674 for (
unsigned int k = 0; k < local.n_cols(); ++k)
675 if (std::fabs(local(j, k)) >= threshold)
685 const unsigned int jcell =
686 this->block_info->local().local_to_global(block_row, j);
687 const unsigned int kcell =
688 this->block_info->local().local_to_global(block_col, k);
690 global.
add(dof1[jcell], dof2[kcell], local(j, k));
696 std::vector<types::global_dof_index> sliced_row_indices(
698 for (
unsigned int i = 0; i < sliced_row_indices.size(); ++i)
699 sliced_row_indices[i] = dof1[bi.
block_start(block_row) + i];
701 std::vector<types::global_dof_index> sliced_col_indices(
703 for (
unsigned int i = 0; i < sliced_col_indices.size(); ++i)
704 sliced_col_indices[i] = dof2[bi.
block_start(block_col) + i];
706 constraints->distribute_local_to_global(local,
714 template <
typename MatrixType,
typename number>
715 template <
class DOFINFO>
720 for (
unsigned int i = 0; i < matrices->size(); ++i)
727 assemble(matrices->block(i),
728 info.matrix(i,
false).matrix,
737 template <
typename MatrixType,
typename number>
738 template <
class DOFINFO>
741 const DOFINFO &info1,
742 const DOFINFO &info2)
744 for (
unsigned int i = 0; i < matrices->size(); ++i)
751 assemble(matrices->block(i),
752 info1.matrix(i,
false).matrix,
757 assemble(matrices->block(i),
758 info1.matrix(i,
true).matrix,
763 assemble(matrices->block(i),
764 info2.matrix(i,
false).matrix,
769 assemble(matrices->block(i),
770 info2.matrix(i,
true).matrix,
781 template <
typename MatrixType,
typename number>
784 : threshold(threshold)
788 template <
typename MatrixType,
typename number>
795 AssertDimension(block_info->local().size(), block_info->global().size());
800 template <
typename MatrixType,
typename number>
805 mg_constrained_dofs = &mg_c;
809 template <
typename MatrixType,
typename number>
810 template <
class DOFINFO>
816 info.initialize_matrices(*matrices, face);
821 template <
typename MatrixType,
typename number>
832 template <
typename MatrixType,
typename number>
842 template <
typename MatrixType,
typename number>
847 const unsigned int block_row,
848 const unsigned int block_col,
849 const std::vector<types::global_dof_index> &dof1,
850 const std::vector<types::global_dof_index> &dof2,
851 const unsigned int level1,
852 const unsigned int level2,
855 for (
unsigned int j = 0; j < local.n_rows(); ++j)
856 for (
unsigned int k = 0; k < local.n_cols(); ++k)
857 if (std::fabs(local(j, k)) >= threshold)
867 const unsigned int jcell =
868 this->block_info->local().local_to_global(block_row, j);
869 const unsigned int kcell =
870 this->block_info->local().local_to_global(block_col, k);
880 const unsigned int jglobal = this->block_info->level(level1)
881 .global_to_local(dof1[jcell])
883 const unsigned int kglobal = this->block_info->level(level2)
884 .global_to_local(dof2[kcell])
887 if (mg_constrained_dofs == 0)
890 global.
add(kglobal, jglobal, local(j, k));
892 global.add(jglobal, kglobal, local(j, k));
896 if (!mg_constrained_dofs->at_refinement_edge(level1,
898 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
900 if (mg_constrained_dofs->set_boundary_values())
902 if ((!mg_constrained_dofs->is_boundary_index(
904 !mg_constrained_dofs->is_boundary_index(
906 (mg_constrained_dofs->is_boundary_index(
908 mg_constrained_dofs->is_boundary_index(
913 global.add(kglobal, jglobal, local(j, k));
915 global.add(jglobal, kglobal, local(j, k));
921 global.add(kglobal, jglobal, local(j, k));
923 global.add(jglobal, kglobal, local(j, k));
931 template <
typename MatrixType,
typename number>
936 const unsigned int block_row,
937 const unsigned int block_col,
938 const std::vector<types::global_dof_index> &dof1,
939 const std::vector<types::global_dof_index> &dof2,
940 const unsigned int level1,
941 const unsigned int level2)
943 for (
unsigned int j = 0; j < local.n_rows(); ++j)
944 for (
unsigned int k = 0; k < local.n_cols(); ++k)
945 if (std::fabs(local(j, k)) >= threshold)
955 const unsigned int jcell =
956 this->block_info->local().local_to_global(block_row, j);
957 const unsigned int kcell =
958 this->block_info->local().local_to_global(block_col, k);
968 const unsigned int jglobal = this->block_info->level(level1)
969 .global_to_local(dof1[jcell])
971 const unsigned int kglobal = this->block_info->level(level2)
972 .global_to_local(dof2[kcell])
975 if (mg_constrained_dofs == 0)
976 global.add(jglobal, kglobal, local(j, k));
979 if (!mg_constrained_dofs->non_refinement_edge_index(
981 !mg_constrained_dofs->non_refinement_edge_index(level2,
984 if (!mg_constrained_dofs->at_refinement_edge(level1,
986 !mg_constrained_dofs->at_refinement_edge(level2,
988 global.add(jglobal, kglobal, local(j, k));
994 template <
typename MatrixType,
typename number>
999 const unsigned int block_row,
1000 const unsigned int block_col,
1001 const std::vector<types::global_dof_index> &dof1,
1002 const std::vector<types::global_dof_index> &dof2,
1003 const unsigned int level1,
1004 const unsigned int level2)
1006 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1007 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1008 if (std::fabs(local(j, k)) >= threshold)
1018 const unsigned int jcell =
1019 this->block_info->local().local_to_global(block_row, j);
1020 const unsigned int kcell =
1021 this->block_info->local().local_to_global(block_col, k);
1031 const unsigned int jglobal = this->block_info->level(level1)
1032 .global_to_local(dof1[jcell])
1034 const unsigned int kglobal = this->block_info->level(level2)
1035 .global_to_local(dof2[kcell])
1038 if (mg_constrained_dofs == 0)
1039 global.add(jglobal, kglobal, local(j, k));
1042 if (!mg_constrained_dofs->non_refinement_edge_index(
1044 !mg_constrained_dofs->non_refinement_edge_index(level2,
1047 if (!mg_constrained_dofs->at_refinement_edge(level1,
1049 !mg_constrained_dofs->at_refinement_edge(level2,
1051 global.add(jglobal, kglobal, local(j, k));
1057 template <
typename MatrixType,
typename number>
1062 const unsigned int block_row,
1063 const unsigned int block_col,
1064 const std::vector<types::global_dof_index> &dof1,
1065 const std::vector<types::global_dof_index> &dof2,
1066 const unsigned int level1,
1067 const unsigned int level2)
1069 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1070 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1071 if (std::fabs(local(k, j)) >= threshold)
1081 const unsigned int jcell =
1082 this->block_info->local().local_to_global(block_row, j);
1083 const unsigned int kcell =
1084 this->block_info->local().local_to_global(block_col, k);
1094 const unsigned int jglobal = this->block_info->level(level1)
1095 .global_to_local(dof1[jcell])
1097 const unsigned int kglobal = this->block_info->level(level2)
1098 .global_to_local(dof2[kcell])
1101 if (mg_constrained_dofs == 0)
1102 global.add(jglobal, kglobal, local(k, j));
1105 if (!mg_constrained_dofs->non_refinement_edge_index(
1107 !mg_constrained_dofs->non_refinement_edge_index(level2,
1110 if (!mg_constrained_dofs->at_refinement_edge(level1,
1112 !mg_constrained_dofs->at_refinement_edge(level2,
1114 global.add(jglobal, kglobal, local(k, j));
1120 template <
typename MatrixType,
typename number>
1125 const unsigned int block_row,
1126 const unsigned int block_col,
1127 const std::vector<types::global_dof_index> &dof1,
1128 const std::vector<types::global_dof_index> &dof2,
1129 const unsigned int level1,
1130 const unsigned int level2)
1135 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1136 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1137 if (std::fabs(local(j, k)) >= threshold)
1147 const unsigned int jcell =
1148 this->block_info->local().local_to_global(block_row, j);
1149 const unsigned int kcell =
1150 this->block_info->local().local_to_global(block_col, k);
1160 const unsigned int jglobal = this->block_info->level(level1)
1161 .global_to_local(dof1[jcell])
1163 const unsigned int kglobal = this->block_info->level(level2)
1164 .global_to_local(dof2[kcell])
1167 if (mg_constrained_dofs == 0)
1168 global.add(jglobal, kglobal, local(j, k));
1171 if (mg_constrained_dofs->at_refinement_edge(level1,
1173 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1175 if (mg_constrained_dofs->set_boundary_values())
1177 if ((!mg_constrained_dofs->is_boundary_index(
1179 !mg_constrained_dofs->is_boundary_index(
1180 level2, kglobal)) ||
1181 (mg_constrained_dofs->is_boundary_index(
1183 mg_constrained_dofs->is_boundary_index(
1185 jglobal == kglobal))
1186 global.add(jglobal, kglobal, local(j, k));
1189 global.add(jglobal, kglobal, local(j, k));
1195 template <
typename MatrixType,
typename number>
1200 const unsigned int block_row,
1201 const unsigned int block_col,
1202 const std::vector<types::global_dof_index> &dof1,
1203 const std::vector<types::global_dof_index> &dof2,
1204 const unsigned int level1,
1205 const unsigned int level2)
1210 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1211 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1212 if (std::fabs(local(k, j)) >= threshold)
1222 const unsigned int jcell =
1223 this->block_info->local().local_to_global(block_row, j);
1224 const unsigned int kcell =
1225 this->block_info->local().local_to_global(block_col, k);
1235 const unsigned int jglobal = this->block_info->level(level1)
1236 .global_to_local(dof1[jcell])
1238 const unsigned int kglobal = this->block_info->level(level2)
1239 .global_to_local(dof2[kcell])
1242 if (mg_constrained_dofs == 0)
1243 global.add(jglobal, kglobal, local(k, j));
1246 if (mg_constrained_dofs->at_refinement_edge(level1,
1248 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1250 if (mg_constrained_dofs->set_boundary_values())
1252 if ((!mg_constrained_dofs->is_boundary_index(
1254 !mg_constrained_dofs->is_boundary_index(
1255 level2, kglobal)) ||
1256 (mg_constrained_dofs->is_boundary_index(
1258 mg_constrained_dofs->is_boundary_index(
1260 jglobal == kglobal))
1261 global.add(jglobal, kglobal, local(k, j));
1264 global.add(jglobal, kglobal, local(k, j));
1271 template <
typename MatrixType,
typename number>
1272 template <
class DOFINFO>
1275 const DOFINFO &info)
1277 const unsigned int level = info.cell->level();
1279 for (
unsigned int i = 0; i < matrices->size(); ++i)
1283 const unsigned int row = matrices->block(i)[
level].row;
1284 const unsigned int col = matrices->block(i)[
level].column;
1286 assemble(matrices->block(i)[
level].matrix,
1287 info.matrix(i,
false).matrix,
1294 if (mg_constrained_dofs != 0)
1296 if (interface_in != 0)
1297 assemble_in(interface_in->block(i)[
level],
1298 info.matrix(i,
false).matrix,
1305 if (interface_out != 0)
1306 assemble_in(interface_out->block(i)[
level],
1307 info.matrix(i,
false).matrix,
1315 assemble_in(matrices->block_in(i)[
level],
1316 info.matrix(i,
false).matrix,
1323 assemble_out(matrices->block_out(i)[
level],
1324 info.matrix(i,
false).matrix,
1336 template <
typename MatrixType,
typename number>
1337 template <
class DOFINFO>
1340 const DOFINFO &info1,
1341 const DOFINFO &info2)
1343 const unsigned int level1 = info1.cell->level();
1344 const unsigned int level2 = info2.cell->level();
1346 for (
unsigned int i = 0; i < matrices->size(); ++i)
1352 const unsigned int row = o[level1].row;
1353 const unsigned int col = o[level1].column;
1355 if (level1 == level2)
1357 if (mg_constrained_dofs == 0)
1359 assemble(o[level1].matrix,
1360 info1.matrix(i,
false).matrix,
1367 assemble(o[level1].matrix,
1368 info1.matrix(i,
true).matrix,
1375 assemble(o[level1].matrix,
1376 info2.matrix(i,
false).matrix,
1383 assemble(o[level1].matrix,
1384 info2.matrix(i,
true).matrix,
1394 assemble_fluxes(o[level1].matrix,
1395 info1.matrix(i,
false).matrix,
1402 assemble_fluxes(o[level1].matrix,
1403 info1.matrix(i,
true).matrix,
1410 assemble_fluxes(o[level1].matrix,
1411 info2.matrix(i,
false).matrix,
1418 assemble_fluxes(o[level1].matrix,
1419 info2.matrix(i,
true).matrix,
1431 if (flux_up->size() != 0)
1436 assemble_fluxes(o[level1].matrix,
1437 info1.matrix(i,
false).matrix,
1444 assemble_up(flux_up->block(i)[level1].matrix,
1445 info1.matrix(i,
true).matrix,
1452 assemble_down(flux_down->block(i)[level1].matrix,
1453 info2.matrix(i,
true).matrix,
size_type block_size(const unsigned int i) const
size_type block_start(const unsigned int i) const
A small class collecting the different BlockIndices involved in global, multilevel and local computat...
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
void add(size_type row, size_type column, const std::string &name)
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
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)
void initialize_info(DOFINFO &info, bool face) const
MatrixPtrVectorPtr 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)
MGMatrixBlockVector< MatrixType > MatrixPtrVector
void assemble(const DOFINFO &info)
ObserverPointer< const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
MatrixPtrVectorPtr interface_out
void initialize(const BlockInfo *block_info, MatrixPtrVector &matrices)
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_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)
ObserverPointer< const MGConstrainedDoFs, MGMatrixLocalBlocksToGlobalBlocks< MatrixType, number > > mg_constrained_dofs
MatrixPtrVectorPtr interface_in
void initialize_interfaces(MatrixPtrVector &interface_in, MatrixPtrVector &interface_out)
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)
MatrixPtrVectorPtr flux_down
void initialize_edge_flux(MatrixPtrVector &up, MatrixPtrVector &down)
MatrixPtrVectorPtr flux_up
MGMatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
void initialize_info(DOFINFO &info, bool face) const
ObserverPointer< const BlockInfo, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > block_info
ObserverPointer< const AffineConstraints< typename MatrixType::value_type >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > constraints
void initialize(const BlockInfo *block_info, MatrixBlockVector< MatrixType > &matrices)
ObserverPointer< MatrixBlockVector< MatrixType >, MatrixLocalBlocksToGlobalBlocks< MatrixType, number > > matrices
MatrixLocalBlocksToGlobalBlocks(double threshold=1.e-12)
void assemble(const DOFINFO &info)
void initialize_info(DOFINFO &info, bool face) const
void initialize(const BlockInfo *block_info, AnyData &residuals)
void assemble(const DOFINFO &info)
ObserverPointer< const BlockInfo, ResidualLocalBlocksToGlobalBlocks< VectorType > > block_info
ObserverPointer< const AffineConstraints< typename VectorType::value_type >, ResidualLocalBlocksToGlobalBlocks< VectorType > > constraints
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)