17 #ifndef dealii_mesh_worker_assembler_h
18 #define dealii_mesh_worker_assembler_h
112 template <
typename VectorType>
136 template <
class DOFINFO>
144 template <
class DOFINFO>
151 template <
class DOFINFO>
153 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
162 const std::vector<types::global_dof_index> &dof);
212 template <
typename MatrixType,
typename number =
double>
244 template <
class DOFINFO>
252 template <
class DOFINFO>
259 template <
class DOFINFO>
261 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
270 const unsigned int block_row,
271 const unsigned int block_col,
272 const std::vector<types::global_dof_index> &dof1,
273 const std::vector<types::global_dof_index> &dof2);
328 template <
typename MatrixType,
typename number =
double>
379 template <
class DOFINFO>
387 template <
class DOFINFO>
394 template <
class DOFINFO>
396 assemble(
const DOFINFO &info1,
const DOFINFO &info2);
405 const unsigned int block_row,
406 const unsigned int block_col,
407 const std::vector<types::global_dof_index> &dof1,
408 const std::vector<types::global_dof_index> &dof2,
409 const unsigned int level1,
410 const unsigned int level2,
419 const unsigned int block_row,
420 const unsigned int block_col,
421 const std::vector<types::global_dof_index> &dof1,
422 const std::vector<types::global_dof_index> &dof2,
423 const unsigned int level1,
424 const unsigned int level2);
432 const unsigned int block_row,
433 const unsigned int block_col,
434 const std::vector<types::global_dof_index> &dof1,
435 const std::vector<types::global_dof_index> &dof2,
436 const unsigned int level1,
437 const unsigned int level2);
445 const unsigned int block_row,
446 const unsigned int block_col,
447 const std::vector<types::global_dof_index> &dof1,
448 const std::vector<types::global_dof_index> &dof2,
449 const unsigned int level1,
450 const unsigned int level2);
458 const unsigned int block_row,
459 const unsigned int block_col,
460 const std::vector<types::global_dof_index> &dof1,
461 const std::vector<types::global_dof_index> &dof2,
462 const unsigned int level1,
463 const unsigned int level2);
471 const unsigned int block_row,
472 const unsigned int block_col,
473 const std::vector<types::global_dof_index> &dof1,
474 const std::vector<types::global_dof_index> &dof2,
475 const unsigned int level1,
476 const unsigned int level2);
532 template <
typename VectorType>
542 template <
typename VectorType>
551 template <
typename VectorType>
552 template <
class DOFINFO>
558 info.initialize_vectors(residuals.size());
561 template <
typename VectorType>
566 const std::vector<types::global_dof_index> &dof)
568 if (constraints == 0)
570 for (
unsigned int b = 0;
b < local.
n_blocks(); ++
b)
571 for (
unsigned int j = 0; j < local.
block(
b).size(); ++j)
581 const unsigned int jcell =
582 this->block_info->local().local_to_global(
b, j);
583 global(dof[jcell]) += local.
block(
b)(j);
587 constraints->distribute_local_to_global(local, dof, global);
591 template <
typename VectorType>
592 template <
class DOFINFO>
596 for (
unsigned int i = 0; i < residuals.size(); ++i)
603 template <
typename VectorType>
604 template <
class DOFINFO>
607 const DOFINFO &info1,
608 const DOFINFO &info2)
610 for (
unsigned int i = 0; i < residuals.size(); ++i)
624 template <
typename MatrixType,
typename number>
627 : threshold(threshold)
631 template <
typename MatrixType,
typename number>
643 template <
typename MatrixType,
typename number>
653 template <
typename MatrixType,
typename number>
654 template <
class DOFINFO>
660 info.initialize_matrices(*matrices, face);
665 template <
typename MatrixType,
typename number>
670 const unsigned int block_row,
671 const unsigned int block_col,
672 const std::vector<types::global_dof_index> &dof1,
673 const std::vector<types::global_dof_index> &dof2)
675 if (constraints ==
nullptr)
677 for (
unsigned int j = 0; j < local.n_rows(); ++j)
678 for (
unsigned int k = 0; k < local.n_cols(); ++k)
689 const unsigned int jcell =
690 this->block_info->local().local_to_global(block_row, j);
691 const unsigned int kcell =
692 this->block_info->local().local_to_global(block_col, k);
694 global.
add(dof1[jcell], dof2[kcell], local(j, k));
700 std::vector<types::global_dof_index> sliced_row_indices(
702 for (
unsigned int i = 0; i < sliced_row_indices.size(); ++i)
703 sliced_row_indices[i] = dof1[bi.
block_start(block_row) + i];
705 std::vector<types::global_dof_index> sliced_col_indices(
707 for (
unsigned int i = 0; i < sliced_col_indices.size(); ++i)
708 sliced_col_indices[i] = dof2[bi.
block_start(block_col) + i];
710 constraints->distribute_local_to_global(local,
718 template <
typename MatrixType,
typename number>
719 template <
class DOFINFO>
724 for (
unsigned int i = 0; i < matrices->size(); ++i)
732 info.matrix(i,
false).matrix,
741 template <
typename MatrixType,
typename number>
742 template <
class DOFINFO>
745 const DOFINFO &info1,
746 const DOFINFO &info2)
748 for (
unsigned int i = 0; i < matrices->size(); ++i)
756 info1.matrix(i,
false).matrix,
762 info1.matrix(i,
true).matrix,
768 info2.matrix(i,
false).matrix,
774 info2.matrix(i,
true).matrix,
785 template <
typename MatrixType,
typename number>
788 : threshold(threshold)
792 template <
typename MatrixType,
typename number>
799 AssertDimension(block_info->local().size(), block_info->global().size());
804 template <
typename MatrixType,
typename number>
809 mg_constrained_dofs = &mg_c;
813 template <
typename MatrixType,
typename number>
814 template <
class DOFINFO>
820 info.initialize_matrices(*matrices, face);
825 template <
typename MatrixType,
typename number>
836 template <
typename MatrixType,
typename number>
846 template <
typename MatrixType,
typename number>
851 const unsigned int block_row,
852 const unsigned int block_col,
853 const std::vector<types::global_dof_index> &dof1,
854 const std::vector<types::global_dof_index> &dof2,
855 const unsigned int level1,
856 const unsigned int level2,
859 for (
unsigned int j = 0; j < local.n_rows(); ++j)
860 for (
unsigned int k = 0; k < local.n_cols(); ++k)
871 const unsigned int jcell =
872 this->block_info->local().local_to_global(block_row, j);
873 const unsigned int kcell =
874 this->block_info->local().local_to_global(block_col, k);
884 const unsigned int jglobal = this->block_info->level(level1)
885 .global_to_local(dof1[jcell])
887 const unsigned int kglobal = this->block_info->level(level2)
888 .global_to_local(dof2[kcell])
891 if (mg_constrained_dofs == 0)
894 global.add(kglobal, jglobal, local(j, k));
896 global.add(jglobal, kglobal, local(j, k));
900 if (!mg_constrained_dofs->at_refinement_edge(level1,
902 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
904 if (mg_constrained_dofs->set_boundary_values())
906 if ((!mg_constrained_dofs->is_boundary_index(
908 !mg_constrained_dofs->is_boundary_index(
910 (mg_constrained_dofs->is_boundary_index(
912 mg_constrained_dofs->is_boundary_index(
917 global.add(kglobal, jglobal, local(j, k));
919 global.add(jglobal, kglobal, local(j, k));
925 global.add(kglobal, jglobal, local(j, k));
927 global.add(jglobal, kglobal, local(j, k));
935 template <
typename MatrixType,
typename number>
940 const unsigned int block_row,
941 const unsigned int block_col,
942 const std::vector<types::global_dof_index> &dof1,
943 const std::vector<types::global_dof_index> &dof2,
944 const unsigned int level1,
945 const unsigned int level2)
947 for (
unsigned int j = 0; j < local.n_rows(); ++j)
948 for (
unsigned int k = 0; k < local.n_cols(); ++k)
959 const unsigned int jcell =
960 this->block_info->local().local_to_global(block_row, j);
961 const unsigned int kcell =
962 this->block_info->local().local_to_global(block_col, k);
972 const unsigned int jglobal = this->block_info->level(level1)
973 .global_to_local(dof1[jcell])
975 const unsigned int kglobal = this->block_info->level(level2)
976 .global_to_local(dof2[kcell])
979 if (mg_constrained_dofs == 0)
980 global.add(jglobal, kglobal, local(j, k));
983 if (!mg_constrained_dofs->non_refinement_edge_index(
985 !mg_constrained_dofs->non_refinement_edge_index(level2,
988 if (!mg_constrained_dofs->at_refinement_edge(level1,
990 !mg_constrained_dofs->at_refinement_edge(level2,
992 global.add(jglobal, kglobal, local(j, k));
998 template <
typename MatrixType,
typename number>
1001 MatrixType & global,
1003 const unsigned int block_row,
1004 const unsigned int block_col,
1005 const std::vector<types::global_dof_index> &dof1,
1006 const std::vector<types::global_dof_index> &dof2,
1007 const unsigned int level1,
1008 const unsigned int level2)
1010 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1011 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1012 if (
std::fabs(local(j, k)) >= threshold)
1022 const unsigned int jcell =
1023 this->block_info->local().local_to_global(block_row, j);
1024 const unsigned int kcell =
1025 this->block_info->local().local_to_global(block_col, k);
1035 const unsigned int jglobal = this->block_info->level(level1)
1036 .global_to_local(dof1[jcell])
1038 const unsigned int kglobal = this->block_info->level(level2)
1039 .global_to_local(dof2[kcell])
1042 if (mg_constrained_dofs == 0)
1043 global.add(jglobal, kglobal, local(j, k));
1046 if (!mg_constrained_dofs->non_refinement_edge_index(
1048 !mg_constrained_dofs->non_refinement_edge_index(level2,
1051 if (!mg_constrained_dofs->at_refinement_edge(level1,
1053 !mg_constrained_dofs->at_refinement_edge(level2,
1055 global.add(jglobal, kglobal, local(j, k));
1061 template <
typename MatrixType,
typename number>
1064 MatrixType & global,
1066 const unsigned int block_row,
1067 const unsigned int block_col,
1068 const std::vector<types::global_dof_index> &dof1,
1069 const std::vector<types::global_dof_index> &dof2,
1070 const unsigned int level1,
1071 const unsigned int level2)
1073 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1074 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1075 if (
std::fabs(local(k, j)) >= threshold)
1085 const unsigned int jcell =
1086 this->block_info->local().local_to_global(block_row, j);
1087 const unsigned int kcell =
1088 this->block_info->local().local_to_global(block_col, k);
1098 const unsigned int jglobal = this->block_info->level(level1)
1099 .global_to_local(dof1[jcell])
1101 const unsigned int kglobal = this->block_info->level(level2)
1102 .global_to_local(dof2[kcell])
1105 if (mg_constrained_dofs == 0)
1106 global.add(jglobal, kglobal, local(k, j));
1109 if (!mg_constrained_dofs->non_refinement_edge_index(
1111 !mg_constrained_dofs->non_refinement_edge_index(level2,
1114 if (!mg_constrained_dofs->at_refinement_edge(level1,
1116 !mg_constrained_dofs->at_refinement_edge(level2,
1118 global.add(jglobal, kglobal, local(k, j));
1124 template <
typename MatrixType,
typename number>
1127 MatrixType & global,
1129 const unsigned int block_row,
1130 const unsigned int block_col,
1131 const std::vector<types::global_dof_index> &dof1,
1132 const std::vector<types::global_dof_index> &dof2,
1133 const unsigned int level1,
1134 const unsigned int level2)
1139 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1140 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1141 if (
std::fabs(local(j, k)) >= threshold)
1151 const unsigned int jcell =
1152 this->block_info->local().local_to_global(block_row, j);
1153 const unsigned int kcell =
1154 this->block_info->local().local_to_global(block_col, k);
1164 const unsigned int jglobal = this->block_info->level(level1)
1165 .global_to_local(dof1[jcell])
1167 const unsigned int kglobal = this->block_info->level(level2)
1168 .global_to_local(dof2[kcell])
1171 if (mg_constrained_dofs == 0)
1172 global.add(jglobal, kglobal, local(j, k));
1175 if (mg_constrained_dofs->at_refinement_edge(level1,
1177 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1179 if (mg_constrained_dofs->set_boundary_values())
1181 if ((!mg_constrained_dofs->is_boundary_index(
1183 !mg_constrained_dofs->is_boundary_index(
1184 level2, kglobal)) ||
1185 (mg_constrained_dofs->is_boundary_index(
1187 mg_constrained_dofs->is_boundary_index(
1189 jglobal == kglobal))
1190 global.add(jglobal, kglobal, local(j, k));
1193 global.add(jglobal, kglobal, local(j, k));
1199 template <
typename MatrixType,
typename number>
1202 MatrixType & global,
1204 const unsigned int block_row,
1205 const unsigned int block_col,
1206 const std::vector<types::global_dof_index> &dof1,
1207 const std::vector<types::global_dof_index> &dof2,
1208 const unsigned int level1,
1209 const unsigned int level2)
1214 for (
unsigned int j = 0; j < local.n_rows(); ++j)
1215 for (
unsigned int k = 0; k < local.n_cols(); ++k)
1216 if (
std::fabs(local(k, j)) >= threshold)
1226 const unsigned int jcell =
1227 this->block_info->local().local_to_global(block_row, j);
1228 const unsigned int kcell =
1229 this->block_info->local().local_to_global(block_col, k);
1239 const unsigned int jglobal = this->block_info->level(level1)
1240 .global_to_local(dof1[jcell])
1242 const unsigned int kglobal = this->block_info->level(level2)
1243 .global_to_local(dof2[kcell])
1246 if (mg_constrained_dofs == 0)
1247 global.add(jglobal, kglobal, local(k, j));
1250 if (mg_constrained_dofs->at_refinement_edge(level1,
1252 !mg_constrained_dofs->at_refinement_edge(level2, kglobal))
1254 if (mg_constrained_dofs->set_boundary_values())
1256 if ((!mg_constrained_dofs->is_boundary_index(
1258 !mg_constrained_dofs->is_boundary_index(
1259 level2, kglobal)) ||
1260 (mg_constrained_dofs->is_boundary_index(
1262 mg_constrained_dofs->is_boundary_index(
1264 jglobal == kglobal))
1265 global.add(jglobal, kglobal, local(k, j));
1268 global.add(jglobal, kglobal, local(k, j));
1275 template <
typename MatrixType,
typename number>
1276 template <
class DOFINFO>
1279 const DOFINFO &info)
1281 const unsigned int level = info.cell->level();
1283 for (
unsigned int i = 0; i < matrices->size(); ++i)
1287 const unsigned int row = matrices->block(i)[
level].row;
1288 const unsigned int col = matrices->block(i)[
level].column;
1291 info.matrix(i,
false).matrix,
1298 if (mg_constrained_dofs != 0)
1300 if (interface_in != 0)
1301 assemble_in(interface_in->block(i)[
level],
1302 info.matrix(i,
false).matrix,
1309 if (interface_out != 0)
1310 assemble_in(interface_out->block(i)[
level],
1311 info.matrix(i,
false).matrix,
1319 assemble_in(matrices->block_in(i)[
level],
1320 info.matrix(i,
false).matrix,
1327 assemble_out(matrices->block_out(i)[
level],
1328 info.matrix(i,
false).matrix,
1340 template <
typename MatrixType,
typename number>
1341 template <
class DOFINFO>
1344 const DOFINFO &info1,
1345 const DOFINFO &info2)
1347 const unsigned int level1 = info1.cell->level();
1348 const unsigned int level2 = info2.cell->level();
1350 for (
unsigned int i = 0; i < matrices->size(); ++i)
1356 const unsigned int row = o[level1].row;
1357 const unsigned int col = o[level1].column;
1359 if (level1 == level2)
1361 if (mg_constrained_dofs == 0)
1364 info1.matrix(i,
false).matrix,
1372 info1.matrix(i,
true).matrix,
1380 info2.matrix(i,
false).matrix,
1388 info2.matrix(i,
true).matrix,
1398 assemble_fluxes(o[level1].
matrix,
1399 info1.matrix(i,
false).matrix,
1406 assemble_fluxes(o[level1].
matrix,
1407 info1.matrix(i,
true).matrix,
1414 assemble_fluxes(o[level1].
matrix,
1415 info2.matrix(i,
false).matrix,
1422 assemble_fluxes(o[level1].
matrix,
1423 info2.matrix(i,
true).matrix,
1435 if (flux_up->size() != 0)
1440 assemble_fluxes(o[level1].
matrix,
1441 info1.matrix(i,
false).matrix,
1448 assemble_up(flux_up->block(i)[level1].matrix,
1449 info1.matrix(i,
true).matrix,
1456 assemble_down(flux_down->block(i)[level1].matrix,
1457 info2.matrix(i,
true).matrix,