18 #ifndef __deal2__matrix_free_h
19 #define __deal2__matrix_free_h
22 #include <deal.II/base/parallel.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/base/vectorization.h>
25 #include <deal.II/base/template_constraints.h>
26 #include <deal.II/fe/fe.h>
27 #include <deal.II/fe/mapping.h>
28 #include <deal.II/fe/mapping_q1.h>
29 #include <deal.II/lac/vector.h>
30 #include <deal.II/lac/parallel_vector.h>
31 #include <deal.II/lac/block_vector_base.h>
32 #include <deal.II/lac/constraint_matrix.h>
33 #include <deal.II/dofs/dof_handler.h>
34 #include <deal.II/multigrid/mg_dof_handler.h>
35 #include <deal.II/hp/dof_handler.h>
36 #include <deal.II/hp/q_collection.h>
37 #include <deal.II/matrix_free/helper_functions.h>
38 #include <deal.II/matrix_free/shape_info.h>
39 #include <deal.II/matrix_free/dof_info.h>
40 #include <deal.II/matrix_free/mapping_info.h>
42 #ifdef DEAL_II_WITH_THREADS
44 #include <tbb/task_scheduler_init.h>
45 #include <tbb/parallel_for.h>
46 #include <tbb/blocked_range.h>
54 DEAL_II_NAMESPACE_OPEN
106 template <
int dim,
typename Number=
double>
299 template <
typename DH,
typename Quadrature>
301 const DH &dof_handler,
311 template <
typename DH,
typename Quadrature>
313 const DH &dof_handler,
322 template <
typename DH,
typename Quadrature>
323 void reinit (
const DH &dof_handler,
354 template <
typename DH,
typename Quadrature>
356 const std::vector<const DH *> &dof_handler,
357 const std::vector<const ConstraintMatrix *> &constraint,
358 const std::vector<IndexSet> &locally_owned_set,
359 const std::vector<Quadrature> &quad,
367 template <
typename DH,
typename Quadrature>
369 const std::vector<const DH *> &dof_handler,
370 const std::vector<const ConstraintMatrix *> &constraint,
371 const std::vector<Quadrature> &quad,
378 template <
typename DH,
typename Quadrature>
379 void reinit (
const std::vector<const DH *> &dof_handler,
380 const std::vector<const ConstraintMatrix *> &constraint,
381 const std::vector<Quadrature> &quad,
391 template <
typename DH,
typename Quadrature>
393 const std::vector<const DH *> &dof_handler,
394 const std::vector<const ConstraintMatrix *> &constraint,
402 template <
typename DH,
typename Quadrature>
403 void reinit (
const std::vector<const DH *> &dof_handler,
404 const std::vector<const ConstraintMatrix *> &constraint,
440 template <
typename OutVector,
typename InVector>
444 const std::pair<
unsigned int,
445 unsigned int> &)> &cell_operation,
447 const InVector &src)
const;
458 template <
typename CLASS,
typename OutVector,
typename InVector>
462 const std::pair<
unsigned int,
463 unsigned int> &)
const,
464 const CLASS *owning_class,
466 const InVector &src)
const;
471 template <
typename CLASS,
typename OutVector,
typename InVector>
475 const std::pair<
unsigned int,
479 const InVector &src)
const;
488 std::pair<unsigned int,unsigned int>
490 const unsigned int fe_degree,
491 const unsigned int vector_component = 0)
const;
499 std::pair<unsigned int,unsigned int>
501 const unsigned int fe_index,
502 const unsigned int vector_component = 0)
const;
518 template <
typename VectorType>
520 const unsigned int vector_component=0)
const;
530 template <
typename Number2>
532 const unsigned int vector_component=0)
const;
544 const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &
568 const std::vector<unsigned int> &
575 void renumber_dofs (std::vector<types::global_dof_index> &renumbering,
576 const unsigned int vector_component = 0);
630 const unsigned int vector_number,
631 const unsigned int fe_component = 0)
const;
644 typename hp::DoFHandler<dim>::active_cell_iterator
646 const unsigned int vector_number,
647 const unsigned int fe_component = 0)
const;
679 const unsigned int hp_active_fe_index = 0)
const;
686 const unsigned int hp_active_fe_index = 0)
const;
694 const unsigned int hp_active_fe_index = 0)
const;
702 const unsigned int hp_active_fe_index = 0)
const;
709 const unsigned int hp_active_fe_index = 0)
const;
716 const unsigned int hp_active_fe_index = 0)
const;
740 template <
typename STREAM>
747 void print (std::ostream &out)
const;
771 get_mapping_info ()
const;
776 const internal::MatrixFreeFunctions::DoFInfo &
777 get_dof_info (
const unsigned int fe_component = 0)
const;
804 const unsigned int quad_index = 0,
805 const unsigned int hp_active_fe_index = 0,
806 const unsigned int hp_active_quad_index = 0)
const;
818 const std::vector<const ConstraintMatrix *> &constraint,
819 const std::vector<IndexSet> &locally_owned_set,
828 const std::vector<const ConstraintMatrix *> &constraint,
829 const std::vector<IndexSet> &locally_owned_set,
841 const std::vector<IndexSet> &locally_owned_set);
847 const unsigned int level);
853 const unsigned int level);
863 std::vector<SmartPointer<const DoFHandler<dim> > > dof_handler;
864 std::vector<SmartPointer<const hp::DoFHandler<dim> > > hp_dof_handler;
865 enum ActiveDoFHandler { usual, hp } active_dof_handler;
866 unsigned int n_dof_handlers;
879 std::vector<internal::MatrixFreeFunctions::DoFInfo>
dof_info;
943 template <
int dim,
typename Number>
944 template <
typename VectorType>
948 const unsigned int comp)
const
951 vec.reinit(dof_info[comp].vector_partitioner->size());
956 template <
int dim,
typename Number>
957 template <
typename Number2>
961 const unsigned int comp)
const
964 vec.
reinit(dof_info[comp].vector_partitioner);
969 template <
int dim,
typename Number>
971 const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &
975 return dof_info[comp].vector_partitioner;
980 template <
int dim,
typename Number>
982 const std::vector<unsigned int> &
986 return dof_info[comp].constrained_dofs;
991 template <
int dim,
typename Number>
997 return dof_handlers.n_dof_handlers;
1002 template <
int dim,
typename Number>
1012 template <
int dim,
typename Number>
1022 template <
int dim,
typename Number>
1027 return size_info.n_macro_cells;
1032 template <
int dim,
typename Number>
1037 return size_info.n_active_cells;
1042 template <
int dim,
typename Number>
1047 return mapping_info;
1052 template <
int dim,
typename Number>
1054 const internal::MatrixFreeFunctions::DoFInfo &
1058 return dof_info[dof_index];
1063 template <
int dim,
typename Number>
1068 return constraint_pool_row_index.size()-1;
1073 template <
int dim,
typename Number>
1079 return constraint_pool_data.empty() ? 0 :
1080 &constraint_pool_data[0] + constraint_pool_row_index[row];
1085 template <
int dim,
typename Number>
1091 return constraint_pool_data.empty() ? 0 :
1092 &constraint_pool_data[0] + constraint_pool_row_index[row+1];
1097 template <
int dim,
typename Number>
1099 std::pair<unsigned int,unsigned int>
1101 (
const std::pair<unsigned int,unsigned int> &range,
1102 const unsigned int degree,
1103 const unsigned int vector_component)
const
1106 if (dof_info[vector_component].cell_active_fe_index.empty())
1108 AssertDimension (dof_info[vector_component].fe_index_conversion.size(),1);
1109 if (dof_info[vector_component].fe_index_conversion[0].first == degree)
1112 return std::pair<unsigned int,unsigned int> (range.second,range.second);
1115 const unsigned int fe_index =
1116 dof_info[vector_component].fe_index_from_degree(degree);
1117 if (fe_index >= dof_info[vector_component].max_fe_index)
1118 return std::pair<unsigned int,unsigned int>(range.second, range.second);
1120 return create_cell_subrange_hp_by_index (range, fe_index, vector_component);
1125 template <
int dim,
typename Number>
1127 std::pair<unsigned int,unsigned int>
1129 (
const std::pair<unsigned int,unsigned int> &range,
1130 const unsigned int fe_index,
1131 const unsigned int vector_component)
const
1134 const std::vector<unsigned int> &fe_indices =
1135 dof_info[vector_component].cell_active_fe_index;
1136 if (fe_indices.size() == 0)
1143 for (
unsigned int i=range.first+1; i<range.second; ++i)
1144 Assert (fe_indices[i] >= fe_indices[i-1],
1145 ExcMessage (
"Cell range must be over sorted range of fe indices in hp case!"));
1149 std::pair<unsigned int,unsigned int> return_range;
1150 return_range.first =
1151 std::lower_bound(&fe_indices[0] + range.first,
1152 &fe_indices[0] + range.second, fe_index)
1154 return_range.second =
1155 std::lower_bound(&fe_indices[0] + return_range.first,
1156 &fe_indices[0] + range.second,
1157 fe_index + 1)-&fe_indices[0];
1158 Assert(return_range.first >= range.first &&
1160 return return_range;
1166 template <
int dim,
typename Number>
1170 const unsigned int vector_component)
1173 dof_info[vector_component].renumber_dofs (renumbering);
1178 template <
int dim,
typename Number>
1184 if (dof_handlers.active_dof_handler == DoFHandlers::usual)
1187 dof_handlers.n_dof_handlers);
1188 return *dof_handlers.dof_handler[dof_index];
1201 template <
int dim,
typename Number>
1205 const unsigned int vector_number,
1206 const unsigned int dof_index)
const
1213 const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
1214 if (irreg_filled > 0)
1219 if (dof_handlers.active_dof_handler == DoFHandlers::usual)
1222 dof_handlers.n_dof_handlers);
1223 dofh = dof_handlers.dof_handler[dof_index];
1228 "for underlying DoFHandler!"));
1231 std::pair<unsigned int,unsigned int> index =
1232 cell_level_index[macro_cell_number*vectorization_length+vector_number];
1234 (&dofh->
get_tria(), index.first, index.second, dofh);
1239 template <
int dim,
typename Number>
1241 typename hp::DoFHandler<dim>::active_cell_iterator
1243 const unsigned int vector_number,
1244 const unsigned int dof_index)
const
1251 const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
1252 if (irreg_filled > 0)
1256 Assert (dof_handlers.active_dof_handler == DoFHandlers::hp,
1259 std::pair<unsigned int,unsigned int> index =
1260 cell_level_index[macro_cell_number*vectorization_length+vector_number];
1261 return typename hp::DoFHandler<dim>::cell_iterator
1262 (&dofh->
get_tria(), index.first, index.second, dofh);
1267 template <
int dim,
typename Number>
1273 return dof_info[0].row_starts[macro_cell][2] > 0;
1278 template <
int dim,
typename Number>
1284 const unsigned int n_filled = dof_info[0].row_starts[macro_cell][2];
1293 template <
int dim,
typename Number>
1297 const unsigned int active_fe_index)
const
1300 return dof_info[dof_index].dofs_per_cell[active_fe_index];
1305 template <
int dim,
typename Number>
1309 const unsigned int active_fe_index)
const
1312 mapping_info.mapping_data_gen.size());
1313 return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index];
1318 template <
int dim,
typename Number>
1322 const unsigned int active_fe_index)
const
1325 return dof_info[dof_index].dofs_per_face[active_fe_index];
1330 template <
int dim,
typename Number>
1334 const unsigned int active_fe_index)
const
1337 mapping_info.mapping_data_gen.size());
1338 return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index];
1343 template <
int dim,
typename Number>
1349 return dof_info[dof_index].vector_partitioner->locally_owned_range();
1354 template <
int dim,
typename Number>
1360 return dof_info[dof_index].vector_partitioner->ghost_indices();
1365 template <
int dim,
typename Number>
1369 const unsigned int index_quad,
1370 const unsigned int active_fe_index,
1371 const unsigned int active_quad_index)
const
1377 return shape_info(index_fe, index_quad,
1378 active_fe_index, active_quad_index);
1383 template <
int dim,
typename Number>
1387 const unsigned int active_fe_index)
const
1390 return mapping_info.mapping_data_gen[quad_index].
1391 quadrature[active_fe_index];
1396 template <
int dim,
typename Number>
1400 const unsigned int active_fe_index)
const
1403 return mapping_info.mapping_data_gen[quad_index].
1404 face_quadrature[active_fe_index];
1409 template <
int dim,
typename Number>
1414 return indices_are_initialized;
1419 template <
int dim,
typename Number>
1424 return mapping_is_initialized;
1435 template <
typename DH>
1437 std::vector<IndexSet>
1438 extract_locally_owned_index_sets (
const std::vector<const DH *> &dofh,
1439 const unsigned int level)
1441 std::vector<IndexSet> locally_owned_set;
1442 locally_owned_set.reserve (dofh.size());
1443 for (
unsigned int j=0; j<dofh.size(); j++)
1445 locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
1449 IndexSet new_set (dofh[j]->n_dofs(level));
1450 new_set.add_range (0, dofh[j]->n_dofs(level));
1451 locally_owned_set.push_back(new_set);
1453 return locally_owned_set;
1460 template <
int dim,
typename Number>
1461 template <
typename DH,
typename Quad>
1463 reinit(
const DH &dof_handler,
1469 std::vector<const DH *> dof_handlers;
1470 std::vector<const ConstraintMatrix *> constraints;
1471 std::vector<Quad> quads;
1473 dof_handlers.push_back(&dof_handler);
1474 constraints.push_back (&constraints_in);
1475 quads.push_back (quad);
1477 std::vector<IndexSet> locally_owned_sets =
1478 internal::MatrixFree::extract_locally_owned_index_sets
1480 reinit(mapping, dof_handlers,constraints, locally_owned_sets, quads,
1486 template <
int dim,
typename Number>
1487 template <
typename DH,
typename Quad>
1490 const DH &dof_handler,
1495 std::vector<const DH *> dof_handlers;
1496 std::vector<const ConstraintMatrix *> constraints;
1497 std::vector<Quad> quads;
1499 dof_handlers.push_back(&dof_handler);
1500 constraints.push_back (&constraints_in);
1501 quads.push_back (quad);
1503 std::vector<IndexSet> locally_owned_sets =
1504 internal::MatrixFree::extract_locally_owned_index_sets
1506 reinit(mapping, dof_handlers,constraints,locally_owned_sets, quads,
1512 template <
int dim,
typename Number>
1513 template <
typename DH,
typename Quad>
1515 reinit(
const std::vector<const DH *> &dof_handler,
1516 const std::vector<const ConstraintMatrix *> &constraint,
1517 const std::vector<Quad> &quad,
1521 std::vector<IndexSet> locally_owned_set =
1522 internal::MatrixFree::extract_locally_owned_index_sets
1524 reinit(mapping, dof_handler,constraint,locally_owned_set,
1525 static_cast<const std::vector<Quadrature<1>
>&>(quad),
1531 template <
int dim,
typename Number>
1532 template <
typename DH,
typename Quad>
1534 reinit(
const std::vector<const DH *> &dof_handler,
1535 const std::vector<const ConstraintMatrix *> &constraint,
1540 std::vector<Quad> quads;
1541 quads.push_back(quad);
1542 std::vector<IndexSet> locally_owned_set =
1543 internal::MatrixFree::extract_locally_owned_index_sets
1545 reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
1551 template <
int dim,
typename Number>
1552 template <
typename DH,
typename Quad>
1555 const std::vector<const DH *> &dof_handler,
1556 const std::vector<const ConstraintMatrix *> &constraint,
1560 std::vector<Quad> quads;
1561 quads.push_back(quad);
1562 std::vector<IndexSet> locally_owned_set =
1563 internal::MatrixFree::extract_locally_owned_index_sets
1565 reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
1571 template <
int dim,
typename Number>
1572 template <
typename DH,
typename Quad>
1575 const std::vector<const DH *> &dof_handler,
1576 const std::vector<const ConstraintMatrix *> &constraint,
1577 const std::vector<Quad> &quad,
1580 std::vector<IndexSet> locally_owned_set =
1581 internal::MatrixFree::extract_locally_owned_index_sets
1583 reinit(mapping, dof_handler,constraint,locally_owned_set,
1584 quad, additional_data);
1598 template <
typename DH>
1600 std::vector<const ::DoFHandler<DH::dimension> *>
1601 resolve_dof_handler (
const std::vector<const DH *> &dof_handler)
1603 std::vector<const ::DoFHandler<DH::dimension> *> conversion(dof_handler.size());
1604 for (
unsigned int i=0; i<dof_handler.size(); ++i)
1605 conversion[i] =
static_cast<const ::DoFHandler<DH::dimension> *
>(dof_handler[i]);
1611 std::vector<const ::hp::DoFHandler<dim> *>
1612 resolve_dof_handler (
const std::vector<const ::hp::DoFHandler<dim> *> &dof_handler)
1621 template <
int dim,
typename Number>
1622 template <
typename DH,
typename Quad>
1625 const std::vector<const DH *> &dof_handler,
1626 const std::vector<const ConstraintMatrix *> &constraint,
1627 const std::vector<IndexSet> &locally_owned_set,
1628 const std::vector<Quad> &quad,
1632 std::vector<hp::QCollection<1> > quad_hp;
1633 for (
unsigned int q=0; q<quad.size(); ++q)
1635 internal_reinit (mapping,
1636 internal::MatrixFree::resolve_dof_handler(dof_handler),
1637 constraint, locally_owned_set, quad_hp, additional_data);
1653 template<
typename VectorStruct>
1654 bool update_ghost_values_start_block (
const VectorStruct &vec,
1655 const unsigned int channel,
1657 template<
typename VectorStruct>
1658 void reset_ghost_values_block (
const VectorStruct &vec,
1659 const bool zero_out_ghosts,
1661 template<
typename VectorStruct>
1662 void update_ghost_values_finish_block (
const VectorStruct &vec,
1664 template<
typename VectorStruct>
1665 void compress_start_block (
const VectorStruct &vec,
1666 const unsigned int channel,
1668 template<
typename VectorStruct>
1669 void compress_finish_block (
const VectorStruct &vec,
1672 template<
typename VectorStruct>
1673 bool update_ghost_values_start_block (
const VectorStruct &,
1679 template<
typename VectorStruct>
1680 void reset_ghost_values_block (
const VectorStruct &,
1684 template<
typename VectorStruct>
1685 void update_ghost_values_finish_block (
const VectorStruct &,
1688 template<
typename VectorStruct>
1689 void compress_start_block (
const VectorStruct &,
1693 template<
typename VectorStruct>
1694 void compress_finish_block (
const VectorStruct &,
1702 template<
typename VectorStruct>
1704 bool update_ghost_values_start (
const VectorStruct &vec,
1705 const unsigned int channel = 0)
1708 update_ghost_values_start_block(vec, channel,
1714 template<
typename Number>
1717 const unsigned int channel = 0)
1721 return return_value;
1726 template <
typename VectorStruct>
1728 bool update_ghost_values_start (
const std::vector<VectorStruct> &vec)
1730 bool return_value =
false;
1731 for (
unsigned int comp=0; comp<vec.size(); comp++)
1732 return_value = update_ghost_values_start(vec[comp], comp);
1733 return return_value;
1738 template <
typename VectorStruct>
1740 bool update_ghost_values_start (
const std::vector<VectorStruct *> &vec)
1742 bool return_value =
false;
1743 for (
unsigned int comp=0; comp<vec.size(); comp++)
1744 return_value = update_ghost_values_start(*vec[comp], comp);
1745 return return_value;
1750 template<
typename VectorStruct>
1752 bool update_ghost_values_start_block (
const VectorStruct &vec,
1753 const unsigned int channel,
1756 bool return_value =
false;
1757 for (
unsigned int i=0; i<vec.n_blocks(); ++i)
1758 return_value = update_ghost_values_start(vec.block(i), channel+509*i);
1759 return return_value;
1767 template<
typename VectorStruct>
1769 void reset_ghost_values (
const VectorStruct &vec,
1770 const bool zero_out_ghosts)
1772 reset_ghost_values_block(vec, zero_out_ghosts,
1778 template<
typename Number>
1781 const bool zero_out_ghosts)
1783 if (zero_out_ghosts)
1789 template <
typename VectorStruct>
1791 void reset_ghost_values (
const std::vector<VectorStruct> &vec,
1792 const bool zero_out_ghosts)
1794 for (
unsigned int comp=0; comp<vec.size(); comp++)
1795 reset_ghost_values(vec[comp], zero_out_ghosts);
1800 template <
typename VectorStruct>
1802 void reset_ghost_values (
const std::vector<VectorStruct *> &vec,
1803 const bool zero_out_ghosts)
1805 for (
unsigned int comp=0; comp<vec.size(); comp++)
1806 reset_ghost_values(*vec[comp], zero_out_ghosts);
1811 template<
typename VectorStruct>
1813 void reset_ghost_values_block (
const VectorStruct &vec,
1814 const bool zero_out_ghosts,
1817 for (
unsigned int i=0; i<vec.n_blocks(); ++i)
1818 reset_ghost_values(vec.block(i), zero_out_ghosts);
1823 template <
typename VectorStruct>
1825 void update_ghost_values_finish (
const VectorStruct &vec)
1827 update_ghost_values_finish_block(vec,
1833 template <
typename Number>
1842 template <
typename VectorStruct>
1844 void update_ghost_values_finish (
const std::vector<VectorStruct> &vec)
1846 for (
unsigned int comp=0; comp<vec.size(); comp++)
1847 update_ghost_values_finish(vec[comp]);
1852 template <
typename VectorStruct>
1854 void update_ghost_values_finish (
const std::vector<VectorStruct *> &vec)
1856 for (
unsigned int comp=0; comp<vec.size(); comp++)
1857 update_ghost_values_finish(*vec[comp]);
1862 template <
typename VectorStruct>
1864 void update_ghost_values_finish_block (
const VectorStruct &vec,
1867 for (
unsigned int i=0; i<vec.n_blocks(); ++i)
1868 update_ghost_values_finish(vec.block(i));
1873 template <
typename VectorStruct>
1875 void compress_start (VectorStruct &vec,
1876 const unsigned int channel = 0)
1878 compress_start_block (vec, channel,
1884 template <
typename Number>
1887 const unsigned int channel = 0)
1894 template <
typename VectorStruct>
1896 void compress_start (std::vector<VectorStruct> &vec)
1898 for (
unsigned int comp=0; comp<vec.size(); comp++)
1899 compress_start (vec[comp], comp);
1904 template <
typename VectorStruct>
1906 void compress_start (std::vector<VectorStruct *> &vec)
1908 for (
unsigned int comp=0; comp<vec.size(); comp++)
1909 compress_start (*vec[comp], comp);
1914 template <
typename VectorStruct>
1916 void compress_start_block (VectorStruct &vec,
1917 const unsigned int channel,
1920 for (
unsigned int i=0; i<vec.n_blocks(); ++i)
1921 compress_start(vec.block(i), channel + 500*i);
1926 template <
typename VectorStruct>
1928 void compress_finish (VectorStruct &vec)
1930 compress_finish_block(vec,
1936 template <
typename Number>
1945 template <
typename VectorStruct>
1947 void compress_finish (std::vector<VectorStruct> &vec)
1949 for (
unsigned int comp=0; comp<vec.size(); comp++)
1950 compress_finish(vec[comp]);
1955 template <
typename VectorStruct>
1957 void compress_finish (std::vector<VectorStruct *> &vec)
1959 for (
unsigned int comp=0; comp<vec.size(); comp++)
1960 compress_finish(*vec[comp]);
1965 template <
typename VectorStruct>
1967 void compress_finish_block (VectorStruct &vec,
1970 for (
unsigned int i=0; i<vec.n_blocks(); ++i)
1971 compress_finish(vec.block(i));
1976 #ifdef DEAL_II_WITH_THREADS
1983 template<
typename Worker,
bool blocked=false>
1984 class CellWork :
public tbb::task
1987 CellWork (
const Worker &worker_in,
1988 const unsigned int partition_in,
1993 task_info (task_info_in)
1995 tbb::task *execute ()
1997 std::pair<unsigned int, unsigned int> cell_range
1998 (task_info.partition_color_blocks_data[partition],
1999 task_info.partition_color_blocks_data[partition+1]);
2002 dummy->spawn (*dummy);
2006 tbb::empty_task *dummy;
2009 const Worker &worker;
2016 template<
typename Worker,
bool blocked=false>
2017 class PartitionWork :
public tbb::task
2020 PartitionWork (
const Worker &function_in,
2021 const unsigned int partition_in,
2024 function (function_in),
2026 task_info (task_info_in)
2028 tbb::task *execute ()
2032 std::pair<unsigned int, unsigned int> cell_range
2033 (task_info.partition_color_blocks_data
2034 [task_info.partition_color_blocks_row_index[partition]],
2035 task_info.partition_color_blocks_data
2036 [task_info.partition_color_blocks_row_index[partition+1]]);
2037 function(cell_range);
2041 tbb::empty_task *root =
new( tbb::task::allocate_root() )
2043 unsigned int evens = task_info.partition_evens[
partition];
2044 unsigned int odds = task_info.partition_odds[
partition];
2045 unsigned int n_blocked_workers =
2046 task_info.partition_n_blocked_workers[
partition];
2047 unsigned int n_workers = task_info.partition_n_workers[
partition];
2048 std::vector<CellWork<Worker,false>*> worker(n_workers);
2049 std::vector<CellWork<Worker,true>*> blocked_worker(n_blocked_workers);
2051 root->set_ref_count(evens+1);
2052 for (
unsigned int j=0; j<evens; j++)
2054 worker[j] =
new(root->allocate_child())
2055 CellWork<Worker,false>(
function,task_info.
2056 partition_color_blocks_row_index
2057 [partition] + 2*j, task_info);
2060 worker[j]->set_ref_count(2);
2061 blocked_worker[j-1]->dummy =
new(worker[j]->allocate_child())
2063 worker[j-1]->spawn(*blocked_worker[j-1]);
2066 worker[j]->set_ref_count(1);
2069 blocked_worker[j] =
new(worker[j]->allocate_child())
2070 CellWork<Worker,true>(
function,task_info.
2071 partition_color_blocks_row_index
2072 [partition] + 2*j+1, task_info);
2078 worker[evens] =
new(worker[j]->allocate_child())
2079 CellWork<Worker,false>(
function,
2081 partition_color_blocks_row_index
2082 [partition]+2*j+1,task_info);
2083 worker[j]->spawn(*worker[evens]);
2087 tbb::empty_task *child =
new(worker[j]->allocate_child())
2089 worker[j]->spawn(*child);
2094 root->wait_for_all();
2095 root->destroy(*root);
2098 dummy->spawn (*dummy);
2102 tbb::empty_task *dummy;
2105 const Worker &
function;
2116 template <
typename Worker>
2120 CellWork (
const Worker &worker_in,
2124 task_info (task_info_in)
2126 void operator()(
const tbb::blocked_range<unsigned int> &r)
const
2128 for (
unsigned int block=r.begin(); block<r.end(); block++)
2130 std::pair<unsigned int,unsigned int> cell_range;
2131 if (task_info.position_short_block<block)
2133 cell_range.first = (block-1)*task_info.block_size+
2134 task_info.block_size_last;
2135 cell_range.second = cell_range.first + task_info.block_size;
2139 cell_range.first = block*task_info.block_size;
2140 cell_range.second = cell_range.first +
2141 ((block == task_info.position_short_block)?
2142 (task_info.block_size_last):(task_info.block_size));
2144 worker (cell_range);
2148 const Worker &worker;
2153 template<
typename Worker,
bool blocked=false>
2154 class PartitionWork :
public tbb::task
2157 PartitionWork (
const Worker &worker_in,
2158 const unsigned int partition_in,
2163 task_info (task_info_in)
2165 tbb::task *execute ()
2167 unsigned int lower = task_info.partition_color_blocks_data[
partition],
2168 upper = task_info.partition_color_blocks_data[
partition+1];
2169 parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
2170 CellWork<Worker> (worker,task_info));
2172 dummy->spawn (*dummy);
2176 tbb::empty_task *dummy;
2179 const Worker &worker;
2187 template<
typename VectorStruct>
2188 class MPIComDistribute :
public tbb::task
2191 MPIComDistribute (
const VectorStruct &src_in)
2196 tbb::task *execute ()
2198 internal::update_ghost_values_finish(src);
2203 const VectorStruct &src;
2208 template<
typename VectorStruct>
2209 class MPIComCompress :
public tbb::task
2212 MPIComCompress (VectorStruct &dst_in)
2217 tbb::task *execute ()
2219 internal::compress_start(dst);
2227 #endif // DEAL_II_WITH_THREADS
2233 template <
int dim,
typename Number>
2234 template <
typename OutVector,
typename InVector>
2241 const std::pair<
unsigned int,
2242 unsigned int> &)> &cell_operation,
2244 const InVector &src)
const
2247 bool ghosts_were_not_set = internal::update_ghost_values_start (src);
2249 #ifdef DEAL_II_WITH_THREADS
2253 if (task_info.use_multithreading ==
true && task_info.n_blocks > 3)
2258 std_cxx1x::function<void (const std::pair<unsigned int,unsigned int> &range)>
2261 const Worker func = std_cxx1x::bind (std_cxx1x::ref(cell_operation),
2262 std_cxx1x::cref(*
this),
2263 std_cxx1x::ref(dst),
2264 std_cxx1x::cref(src),
2267 if (task_info.use_partition_partition ==
true)
2269 tbb::empty_task *root =
new( tbb::task::allocate_root() )
2271 unsigned int evens = task_info.evens;
2272 unsigned int odds = task_info.odds;
2273 root->set_ref_count(evens+1);
2274 unsigned int n_blocked_workers = task_info.n_blocked_workers;
2275 unsigned int n_workers = task_info.n_workers;
2276 std::vector<internal::partition::PartitionWork<Worker,false>*>
2278 std::vector<internal::partition::PartitionWork<Worker,true>*>
2279 blocked_worker(n_blocked_workers);
2280 internal::MPIComCompress<OutVector> *worker_compr =
2281 new(root->allocate_child())
2282 internal::MPIComCompress<OutVector>(dst);
2283 worker_compr->set_ref_count(1);
2284 for (
unsigned int j=0; j<evens; j++)
2288 worker[j] =
new(root->allocate_child())
2289 internal::partition::PartitionWork<Worker,false>
2290 (func,2*j,task_info);
2291 worker[j]->set_ref_count(2);
2292 blocked_worker[j-1]->dummy =
new(worker[j]->allocate_child())
2295 worker[j-1]->spawn(*blocked_worker[j-1]);
2297 worker_compr->spawn(*blocked_worker[j-1]);
2301 worker[j] =
new(worker_compr->allocate_child())
2302 internal::partition::PartitionWork<Worker,false>
2303 (func,2*j,task_info);
2304 worker[j]->set_ref_count(2);
2305 internal::MPIComDistribute<InVector> *worker_dist =
2306 new (worker[j]->allocate_child())
2307 internal::MPIComDistribute<InVector>(src);
2308 worker_dist->spawn(*worker_dist);
2312 blocked_worker[j] =
new(worker[j]->allocate_child())
2313 internal::partition::PartitionWork<Worker,true>
2314 (func,2*j+1,task_info);
2320 worker[evens] =
new(worker[j]->allocate_child())
2321 internal::partition::PartitionWork<Worker,false>
2322 (func,2*j+1,task_info);
2323 worker[j]->spawn(*worker[evens]);
2327 tbb::empty_task *child =
new(worker[j]->allocate_child())
2329 worker[j]->spawn(*child);
2334 root->wait_for_all();
2335 root->destroy(*root);
2339 unsigned int evens = task_info.evens;
2340 unsigned int odds = task_info.odds;
2346 tbb::empty_task *root =
new( tbb::task::allocate_root() ) tbb::empty_task;
2347 root->set_ref_count(evens+1);
2348 unsigned int n_blocked_workers = odds-(odds+evens+1)%2;
2349 unsigned int n_workers = task_info.partition_color_blocks_data.size()-1-
2351 std::vector<internal::color::PartitionWork<Worker,false>*> worker(n_workers);
2352 std::vector<internal::color::PartitionWork<Worker,true>*> blocked_worker(n_blocked_workers);
2353 unsigned int worker_index = 0, slice_index = 0;
2354 unsigned int spawn_index = 0, spawn_index_new = 0;
2355 int spawn_index_child = -2;
2356 internal::MPIComCompress<OutVector> *worker_compr =
new(root->allocate_child())
2357 internal::MPIComCompress<OutVector>(dst);
2358 worker_compr->set_ref_count(1);
2359 for (
unsigned int part=0;
2360 part<task_info.partition_color_blocks_row_index.size()-1; part++)
2362 spawn_index_new = worker_index;
2364 worker[worker_index] =
new(worker_compr->allocate_child())
2365 internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2367 worker[worker_index] =
new(root->allocate_child())
2368 internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2370 for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
2373 worker[worker_index]->set_ref_count(1);
2375 worker[worker_index] =
new (worker[worker_index-1]->allocate_child())
2376 internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2378 worker[worker_index]->set_ref_count(2);
2381 blocked_worker[(part-1)/2]->dummy =
2382 new (worker[worker_index]->allocate_child()) tbb::empty_task;
2384 if (spawn_index_child == -1)
2385 worker[spawn_index]->spawn(*blocked_worker[(part-1)/2]);
2387 worker[spawn_index]->spawn(*worker[spawn_index_child]);
2388 spawn_index = spawn_index_new;
2389 spawn_index_child = -2;
2393 internal::MPIComDistribute<InVector> *worker_dist =
2394 new (worker[worker_index]->allocate_child())
2395 internal::MPIComDistribute<InVector>(src);
2396 worker_dist->spawn(*worker_dist);
2400 if (part<task_info.partition_color_blocks_row_index.size()-1)
2402 if (part<task_info.partition_color_blocks_row_index.size()-2)
2404 blocked_worker[part/2] =
new(worker[worker_index-1]->allocate_child())
2405 internal::color::PartitionWork<Worker,true>(func,slice_index,task_info);
2408 task_info.partition_color_blocks_row_index[part+1])
2410 blocked_worker[part/2]->set_ref_count(1);
2411 worker[worker_index] =
new(blocked_worker[part/2]->allocate_child())
2412 internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2417 spawn_index_child = -1;
2421 for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
2425 task_info.partition_color_blocks_row_index[part])
2427 worker[worker_index]->set_ref_count(1);
2430 worker[worker_index] =
new (worker[worker_index-1]->allocate_child())
2431 internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2433 spawn_index_child = worker_index;
2438 tbb::empty_task *
final =
new (worker[worker_index-1]->allocate_child())
2440 worker[spawn_index]->spawn(*
final);
2441 spawn_index_child = worker_index-1;
2445 worker[spawn_index]->spawn(*worker[spawn_index_child]);
2446 root->wait_for_all();
2447 root->destroy(*root);
2454 internal::update_ghost_values_finish(src);
2456 for (
unsigned int color=0;
2457 color < task_info.partition_color_blocks_row_index[1];
2460 unsigned int lower = task_info.partition_color_blocks_data[color],
2461 upper = task_info.partition_color_blocks_data[color+1];
2462 parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
2463 internal::color::CellWork<Worker>
2467 internal::compress_start(dst);
2475 std::pair<unsigned int,unsigned int> cell_range;
2479 cell_range.first = 0;
2480 cell_range.second = size_info.boundary_cells_start;
2481 cell_operation (*
this, dst, src, cell_range);
2486 internal::update_ghost_values_finish(src);
2489 if (size_info.boundary_cells_end > size_info.boundary_cells_start)
2491 cell_range.first = size_info.boundary_cells_start;
2492 cell_range.second = size_info.boundary_cells_end;
2493 cell_operation (*
this, dst, src, cell_range);
2496 internal::compress_start(dst);
2499 if (size_info.n_macro_cells > size_info.boundary_cells_end)
2501 cell_range.first = size_info.boundary_cells_end;
2502 cell_range.second = size_info.n_macro_cells;
2503 cell_operation (*
this, dst, src, cell_range);
2508 internal::compress_finish(dst);
2509 internal::reset_ghost_values(src, ghosts_were_not_set);
2514 template <
int dim,
typename Number>
2515 template <
typename CLASS,
typename OutVector,
typename InVector>
2522 const std::pair<
unsigned int,
2523 unsigned int> &)
const,
2524 const CLASS *owning_class,
2526 const InVector &src)
const
2530 std_cxx1x::function<void (const MatrixFree<dim,Number> &,
2533 const std::pair<
unsigned int,
2535 function = std_cxx1x::bind<void>(function_pointer,
2536 std_cxx1x::cref(*owning_class),
2541 cell_loop (
function, dst, src);
2546 template <
int dim,
typename Number>
2547 template <
typename CLASS,
typename OutVector,
typename InVector>
2554 const std::pair<
unsigned int,
2556 CLASS *owning_class,
2558 const InVector &src)
const
2562 std_cxx1x::function<void (const MatrixFree<dim,Number> &,
2565 const std::pair<
unsigned int,
2567 function = std_cxx1x::bind<void>(function_pointer,
2568 std_cxx1x::ref(*owning_class),
2573 cell_loop (
function, dst, src);
2577 #endif // ifndef DOXYGEN
2581 DEAL_II_NAMESPACE_CLOSE
Transformed quadrature weights.
void internal_reinit(const Mapping< dim > &mapping, const std::vector< const DoFHandler< dim > * > &dof_handler, const std::vector< const ConstraintMatrix * > &constraint, const std::vector< IndexSet > &locally_owned_set, const std::vector< hp::QCollection< 1 > > &quad, const AdditionalData additional_data)
unsigned int get_n_q_points_face(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
internal::MatrixFreeFunctions::TaskInfo task_info
const IndexSet & get_ghost_set(const unsigned int fe_component=0) const
bool indices_are_initialized
static const unsigned int invalid_unsigned_int
bool at_irregular_cell(const unsigned int macro_cell_number) const
std::vector< unsigned int > constraint_pool_row_index
#define AssertDimension(dim1, dim2)
AdditionalData(const MPI_Comm mpi_communicator=MPI_COMM_SELF, const TasksParallelScheme tasks_parallel_scheme=partition_partition, const unsigned int tasks_block_size=0, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const unsigned int level_mg_handler=numbers::invalid_unsigned_int, const bool store_plain_indices=true, const bool initialize_indices=true, const bool initialize_mapping=true)
internal::MatrixFreeFunctions::MappingInfo< dim, Number > mapping_info
const DoFHandler< dim > & get_dof_handler(const unsigned int fe_component=0) const
void compress_start(const unsigned int communication_channel=0,::VectorOperation::values operation=VectorOperation::add)
unsigned int get_n_q_points(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int fe_component=0) const
bool mapping_initialized() const
void initialize_dof_handlers(const std::vector< const DoFHandler< dim > * > &dof_handlers, const unsigned int level)
void compress_finish(::VectorOperation::values operation)
::ExceptionBase & ExcMessage(std::string arg1)
void reinit(const size_type size, const bool fast=false)
std::vector< Number > constraint_pool_data
#define AssertIndexRange(index, range)
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
void copy_from(const MatrixFree< dim, Number > &matrix_free_base)
unsigned int n_constraint_pool_entries() const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
const std_cxx1x::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int vector_component=0) const
const IndexSet & get_locally_owned_set(const unsigned int fe_component=0) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int fe_component=0, const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
bool mapping_is_initialized
unsigned int n_macro_cells() const
internal::MatrixFreeFunctions::SizeInfo size_info
void initialize_dof_vector(VectorType &vec, const unsigned int vector_component=0) const
void print_memory_consumption(STREAM &out) const
#define Assert(cond, exc)
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int fe_component=0) const
unsigned int n_components() const
bool has_ghost_elements() const
bool indices_initialized() const
void update_ghost_values_finish() const
const internal::MatrixFreeFunctions::SizeInfo & get_size_info() const
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
TasksParallelScheme tasks_parallel_scheme
const Quadrature< dim-1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
MPI_Comm mpi_communicator
void update_ghost_values_start(const unsigned int communication_channel=0) const
const Number * constraint_pool_end(const unsigned int pool_index) const
std::size_t memory_consumption() const
unsigned int tasks_block_size
unsigned int get_dofs_per_cell(const unsigned int fe_component=0, const unsigned int hp_active_fe_index=0) const
const Triangulation< dim, spacedim > & get_tria() const
unsigned int n_physical_cells() const
unsigned int get_dofs_per_face(const unsigned int fe_component=0, const unsigned int hp_active_fe_index=0) const
void reinit(const Mapping< dim > &mapping, const DH &dof_handler, const ConstraintMatrix &constraint, const IndexSet &locally_owned_dofs, const Quadrature &quad, const AdditionalData additional_data=AdditionalData())
unsigned int level_mg_handler
std::pair< unsigned int, unsigned int > create_cell_subrange_hp(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_degree, const unsigned int vector_component=0) const
Shape function gradients.
hp::DoFHandler< dim >::active_cell_iterator get_hp_cell_iterator(const unsigned int macro_cell_number, const unsigned int vector_number, const unsigned int fe_component=0) const
UpdateFlags mapping_update_flags
::ExceptionBase & ExcNotImplemented()
const Number * constraint_pool_begin(const unsigned int pool_index) const
const Triangulation< dim, spacedim > & get_tria() const
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
::ExceptionBase & ExcInternalError()
const std::vector< unsigned int > & get_constrained_dofs(const unsigned int fe_component=0) const
void renumber_dofs(std::vector< types::global_dof_index > &renumbering, const unsigned int vector_component=0)
void initialize_indices(const std::vector< const ConstraintMatrix * > &constraint, const std::vector< IndexSet > &locally_owned_set)
std::pair< unsigned int, unsigned int > create_cell_subrange_hp_by_index(const std::pair< unsigned int, unsigned int > &range, const unsigned int fe_index, const unsigned int vector_component=0) const
void print(std::ostream &out) const
unsigned int n_components_filled(const unsigned int macro_cell_number) const
void cell_loop(const std_cxx1x::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src) const