Reference documentation for deal.II version 8.4.1
matrix_free.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2016 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii__matrix_free_h
18 #define dealii__matrix_free_h
19 
20 #include <deal.II/base/exceptions.h>
21 #include <deal.II/base/parallel.h>
22 #include <deal.II/base/quadrature.h>
23 #include <deal.II/base/vectorization.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/fe/fe.h>
26 #include <deal.II/fe/mapping.h>
27 #include <deal.II/fe/mapping_q1.h>
28 #include <deal.II/lac/vector.h>
29 #include <deal.II/lac/parallel_vector.h>
30 #include <deal.II/lac/block_vector_base.h>
31 #include <deal.II/lac/constraint_matrix.h>
32 #include <deal.II/dofs/dof_handler.h>
33 #include <deal.II/hp/dof_handler.h>
34 #include <deal.II/hp/q_collection.h>
35 #include <deal.II/matrix_free/helper_functions.h>
36 #include <deal.II/matrix_free/shape_info.h>
37 #include <deal.II/matrix_free/dof_info.h>
38 #include <deal.II/matrix_free/mapping_info.h>
39 
40 #ifdef DEAL_II_WITH_THREADS
41 #include <tbb/task.h>
42 #include <tbb/task_scheduler_init.h>
43 #include <tbb/parallel_for.h>
44 #include <tbb/blocked_range.h>
45 #endif
46 
47 #include <stdlib.h>
48 #include <memory>
49 #include <limits>
50 
51 
52 DEAL_II_NAMESPACE_OPEN
53 
54 
55 
103 template <int dim, typename Number=double>
105 {
106 public:
107 
137  {
141  enum TasksParallelScheme {none, partition_partition, partition_color, color};
142 
146  AdditionalData (const MPI_Comm mpi_communicator = MPI_COMM_SELF,
147  const TasksParallelScheme tasks_parallel_scheme = partition_partition,
148  const unsigned int tasks_block_size = 0,
150  const unsigned int level_mg_handler = numbers::invalid_unsigned_int,
151  const bool store_plain_indices = true,
152  const bool initialize_indices = true,
153  const bool initialize_mapping = true)
154  :
163  {};
164 
171  MPI_Comm mpi_communicator;
172 
202 
212  unsigned int tasks_block_size;
213 
226 
235  unsigned int level_mg_handler;
236 
244 
254 
262  };
263 
271  MatrixFree ();
272 
276  ~MatrixFree();
277 
295  template <typename DoFHandlerType, typename QuadratureType>
296  void reinit (const Mapping<dim> &mapping,
297  const DoFHandlerType &dof_handler,
298  const ConstraintMatrix &constraint,
299  const IndexSet &locally_owned_dofs,
300  const QuadratureType &quad,
301  const AdditionalData additional_data = AdditionalData());
302 
307  template <typename DoFHandlerType, typename QuadratureType>
308  void reinit (const Mapping<dim> &mapping,
309  const DoFHandlerType &dof_handler,
310  const ConstraintMatrix &constraint,
311  const QuadratureType &quad,
312  const AdditionalData additional_data = AdditionalData());
313 
318  template <typename DoFHandlerType, typename QuadratureType>
319  void reinit (const DoFHandlerType &dof_handler,
320  const ConstraintMatrix &constraint,
321  const QuadratureType &quad,
322  const AdditionalData additional_data = AdditionalData());
323 
350  template <typename DoFHandlerType, typename QuadratureType>
351  void reinit (const Mapping<dim> &mapping,
352  const std::vector<const DoFHandlerType *> &dof_handler,
353  const std::vector<const ConstraintMatrix *> &constraint,
354  const std::vector<IndexSet> &locally_owned_set,
355  const std::vector<QuadratureType> &quad,
356  const AdditionalData additional_data = AdditionalData());
357 
363  template <typename DoFHandlerType, typename QuadratureType>
364  void reinit (const Mapping<dim> &mapping,
365  const std::vector<const DoFHandlerType *> &dof_handler,
366  const std::vector<const ConstraintMatrix *> &constraint,
367  const std::vector<QuadratureType> &quad,
368  const AdditionalData additional_data = AdditionalData());
369 
374  template <typename DoFHandlerType, typename QuadratureType>
375  void reinit (const std::vector<const DoFHandlerType *> &dof_handler,
376  const std::vector<const ConstraintMatrix *> &constraint,
377  const std::vector<QuadratureType> &quad,
378  const AdditionalData additional_data = AdditionalData());
379 
387  template <typename DoFHandlerType, typename QuadratureType>
388  void reinit (const Mapping<dim> &mapping,
389  const std::vector<const DoFHandlerType *> &dof_handler,
390  const std::vector<const ConstraintMatrix *> &constraint,
391  const QuadratureType &quad,
392  const AdditionalData additional_data = AdditionalData());
393 
398  template <typename DoFHandlerType, typename QuadratureType>
399  void reinit (const std::vector<const DoFHandlerType *> &dof_handler,
400  const std::vector<const ConstraintMatrix *> &constraint,
401  const QuadratureType &quad,
402  const AdditionalData additional_data = AdditionalData());
403 
409  void copy_from (const MatrixFree<dim,Number> &matrix_free_base);
410 
415  void clear();
416 
418 
436  template <typename OutVector, typename InVector>
437  void cell_loop (const std_cxx11::function<void (const MatrixFree<dim,Number> &,
438  OutVector &,
439  const InVector &,
440  const std::pair<unsigned int,
441  unsigned int> &)> &cell_operation,
442  OutVector &dst,
443  const InVector &src) const;
444 
454  template <typename CLASS, typename OutVector, typename InVector>
455  void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
456  OutVector &,
457  const InVector &,
458  const std::pair<unsigned int,
459  unsigned int> &)const,
460  const CLASS *owning_class,
461  OutVector &dst,
462  const InVector &src) const;
463 
467  template <typename CLASS, typename OutVector, typename InVector>
468  void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
469  OutVector &,
470  const InVector &,
471  const std::pair<unsigned int,
472  unsigned int> &),
473  CLASS *owning_class,
474  OutVector &dst,
475  const InVector &src) const;
476 
484  std::pair<unsigned int,unsigned int>
485  create_cell_subrange_hp (const std::pair<unsigned int,unsigned int> &range,
486  const unsigned int fe_degree,
487  const unsigned int vector_component = 0) const;
488 
495  std::pair<unsigned int,unsigned int>
496  create_cell_subrange_hp_by_index (const std::pair<unsigned int,unsigned int> &range,
497  const unsigned int fe_index,
498  const unsigned int vector_component = 0) const;
499 
501 
514  template <typename VectorType>
515  void initialize_dof_vector(VectorType &vec,
516  const unsigned int vector_component=0) const;
517 
526  template <typename Number2>
528  const unsigned int vector_component=0) const;
529 
540  const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> &
541  get_vector_partitioner (const unsigned int vector_component=0) const;
542 
546  const IndexSet &
547  get_locally_owned_set (const unsigned int fe_component = 0) const;
548 
552  const IndexSet &
553  get_ghost_set (const unsigned int fe_component = 0) const;
554 
564  const std::vector<unsigned int> &
565  get_constrained_dofs (const unsigned int fe_component = 0) const;
566 
571  void renumber_dofs (std::vector<types::global_dof_index> &renumbering,
572  const unsigned int vector_component = 0);
573 
575 
583  unsigned int n_components () const;
584 
593  unsigned int n_physical_cells () const;
594 
605  unsigned int n_macro_cells () const;
606 
611  const DoFHandler<dim> &
612  get_dof_handler (const unsigned int fe_component = 0) const;
613 
625  get_cell_iterator (const unsigned int macro_cell_number,
626  const unsigned int vector_number,
627  const unsigned int fe_component = 0) const;
628 
641  get_hp_cell_iterator (const unsigned int macro_cell_number,
642  const unsigned int vector_number,
643  const unsigned int fe_component = 0) const;
644 
656  bool
657  at_irregular_cell (const unsigned int macro_cell_number) const;
658 
667  unsigned int
668  n_components_filled (const unsigned int macro_cell_number) const;
669 
673  unsigned int
674  get_dofs_per_cell (const unsigned int fe_component = 0,
675  const unsigned int hp_active_fe_index = 0) const;
676 
680  unsigned int
681  get_n_q_points (const unsigned int quad_index = 0,
682  const unsigned int hp_active_fe_index = 0) const;
683 
688  unsigned int
689  get_dofs_per_face (const unsigned int fe_component = 0,
690  const unsigned int hp_active_fe_index = 0) const;
691 
696  unsigned int
697  get_n_q_points_face (const unsigned int quad_index = 0,
698  const unsigned int hp_active_fe_index = 0) const;
699 
703  const Quadrature<dim> &
704  get_quadrature (const unsigned int quad_index = 0,
705  const unsigned int hp_active_fe_index = 0) const;
706 
710  const Quadrature<dim-1> &
711  get_face_quadrature (const unsigned int quad_index = 0,
712  const unsigned int hp_active_fe_index = 0) const;
713 
717  bool indices_initialized () const;
718 
724  bool mapping_initialized () const;
725 
730  std::size_t memory_consumption() const;
731 
736  template <typename StreamType>
737  void print_memory_consumption(StreamType &out) const;
738 
743  void print (std::ostream &out) const;
744 
746 
755  get_task_info () const;
756 
761  get_size_info () const;
762 
763  /*
764  * Returns geometry-dependent information on the cells.
765  */
767  get_mapping_info () const;
768 
772  const internal::MatrixFreeFunctions::DoFInfo &
773  get_dof_info (const unsigned int fe_component = 0) const;
774 
778  unsigned int n_constraint_pool_entries() const;
779 
784  const Number *
785  constraint_pool_begin (const unsigned int pool_index) const;
786 
792  const Number *
793  constraint_pool_end (const unsigned int pool_index) const;
794 
799  get_shape_info (const unsigned int fe_component = 0,
800  const unsigned int quad_index = 0,
801  const unsigned int hp_active_fe_index = 0,
802  const unsigned int hp_active_quad_index = 0) const;
803 
805 
806 private:
807 
812  void internal_reinit (const Mapping<dim> &mapping,
813  const std::vector<const DoFHandler<dim> *> &dof_handler,
814  const std::vector<const ConstraintMatrix *> &constraint,
815  const std::vector<IndexSet> &locally_owned_set,
816  const std::vector<hp::QCollection<1> > &quad,
817  const AdditionalData additional_data);
818 
822  void internal_reinit (const Mapping<dim> &mapping,
823  const std::vector<const hp::DoFHandler<dim>*> &dof_handler,
824  const std::vector<const ConstraintMatrix *> &constraint,
825  const std::vector<IndexSet> &locally_owned_set,
826  const std::vector<hp::QCollection<1> > &quad,
827  const AdditionalData additional_data);
828 
835  void
836  initialize_indices (const std::vector<const ConstraintMatrix *> &constraint,
837  const std::vector<IndexSet> &locally_owned_set);
838 
842  void initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handlers,
843  const unsigned int level);
844 
848  void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
849  const unsigned int level);
850 
856  struct DoFHandlers
857  {
858  DoFHandlers () : n_dof_handlers (0), level (numbers::invalid_unsigned_int) {};
859  std::vector<SmartPointer<const DoFHandler<dim> > > dof_handler;
860  std::vector<SmartPointer<const hp::DoFHandler<dim> > > hp_dof_handler;
861  enum ActiveDoFHandler { usual, hp } active_dof_handler;
862  unsigned int n_dof_handlers;
863  unsigned int level;
864  };
865 
870 
875  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
876 
883  std::vector<Number> constraint_pool_data;
884 
889  std::vector<unsigned int> constraint_pool_row_index;
890 
896 
901 
908  std::vector<std::pair<unsigned int,unsigned int> > cell_level_index;
909 
915 
920 
925 
930 };
931 
932 
933 
934 /*----------------------- Inline functions ----------------------------------*/
935 
936 #ifndef DOXYGEN
937 
938 
939 template <int dim, typename Number>
940 template <typename VectorType>
941 inline
942 void
944  const unsigned int comp) const
945 {
946  AssertIndexRange (comp, n_components());
947  vec.reinit(dof_info[comp].vector_partitioner->size());
948 }
949 
950 
951 
952 template <int dim, typename Number>
953 template <typename Number2>
954 inline
955 void
957  const unsigned int comp) const
958 {
959  AssertIndexRange (comp, n_components());
960  vec.reinit(dof_info[comp].vector_partitioner);
961 }
962 
963 
964 
965 template <int dim, typename Number>
966 inline
967 const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> &
968 MatrixFree<dim,Number>::get_vector_partitioner (const unsigned int comp) const
969 {
970  AssertIndexRange (comp, n_components());
971  return dof_info[comp].vector_partitioner;
972 }
973 
974 
975 
976 template <int dim, typename Number>
977 inline
978 const std::vector<unsigned int> &
979 MatrixFree<dim,Number>::get_constrained_dofs (const unsigned int comp) const
980 {
981  AssertIndexRange (comp, n_components());
982  return dof_info[comp].constrained_dofs;
983 }
984 
985 
986 
987 template <int dim, typename Number>
988 inline
989 unsigned int
991 {
992  AssertDimension (dof_handlers.n_dof_handlers, dof_info.size());
993  return dof_handlers.n_dof_handlers;
994 }
995 
996 
997 
998 template <int dim, typename Number>
999 inline
1002 {
1003  return task_info;
1004 }
1005 
1006 
1007 
1008 template <int dim, typename Number>
1009 inline
1012 {
1013  return size_info;
1014 }
1015 
1016 
1017 
1018 template <int dim, typename Number>
1019 inline
1020 unsigned int
1022 {
1023  return size_info.n_macro_cells;
1024 }
1025 
1026 
1027 
1028 template <int dim, typename Number>
1029 inline
1030 unsigned int
1032 {
1033  return size_info.n_active_cells;
1034 }
1035 
1036 
1037 
1038 template <int dim, typename Number>
1039 inline
1042 {
1043  return mapping_info;
1044 }
1045 
1046 
1047 
1048 template <int dim, typename Number>
1049 inline
1050 const internal::MatrixFreeFunctions::DoFInfo &
1051 MatrixFree<dim,Number>::get_dof_info (unsigned int dof_index) const
1052 {
1053  AssertIndexRange (dof_index, n_components());
1054  return dof_info[dof_index];
1055 }
1056 
1057 
1058 
1059 template <int dim, typename Number>
1060 inline
1061 unsigned int
1063 {
1064  return constraint_pool_row_index.size()-1;
1065 }
1066 
1067 
1068 
1069 template <int dim, typename Number>
1070 inline
1071 const Number *
1072 MatrixFree<dim,Number>::constraint_pool_begin (const unsigned int row) const
1073 {
1074  AssertIndexRange (row, constraint_pool_row_index.size()-1);
1075  return constraint_pool_data.empty() ? 0 :
1076  &constraint_pool_data[0] + constraint_pool_row_index[row];
1077 }
1078 
1079 
1080 
1081 template <int dim, typename Number>
1082 inline
1083 const Number *
1084 MatrixFree<dim,Number>::constraint_pool_end (const unsigned int row) const
1085 {
1086  AssertIndexRange (row, constraint_pool_row_index.size()-1);
1087  return constraint_pool_data.empty() ? 0 :
1088  &constraint_pool_data[0] + constraint_pool_row_index[row+1];
1089 }
1090 
1091 
1092 
1093 template <int dim, typename Number>
1094 inline
1095 std::pair<unsigned int,unsigned int>
1097 (const std::pair<unsigned int,unsigned int> &range,
1098  const unsigned int degree,
1099  const unsigned int vector_component) const
1100 {
1101  AssertIndexRange (vector_component, dof_info.size());
1102  if (dof_info[vector_component].cell_active_fe_index.empty())
1103  {
1104  AssertDimension (dof_info[vector_component].fe_index_conversion.size(),1);
1105  if (dof_info[vector_component].fe_index_conversion[0].first == degree)
1106  return range;
1107  else
1108  return std::pair<unsigned int,unsigned int> (range.second,range.second);
1109  }
1110 
1111  const unsigned int fe_index =
1112  dof_info[vector_component].fe_index_from_degree(degree);
1113  if (fe_index >= dof_info[vector_component].max_fe_index)
1114  return std::pair<unsigned int,unsigned int>(range.second, range.second);
1115  else
1116  return create_cell_subrange_hp_by_index (range, fe_index, vector_component);
1117 }
1118 
1119 
1120 
1121 template <int dim, typename Number>
1122 inline
1123 std::pair<unsigned int,unsigned int>
1125 (const std::pair<unsigned int,unsigned int> &range,
1126  const unsigned int fe_index,
1127  const unsigned int vector_component) const
1128 {
1129  AssertIndexRange (fe_index, dof_info[vector_component].max_fe_index);
1130  const std::vector<unsigned int> &fe_indices =
1131  dof_info[vector_component].cell_active_fe_index;
1132  if (fe_indices.size() == 0)
1133  return range;
1134  else
1135  {
1136  // the range over which we are searching must be ordered, otherwise we
1137  // got a range that spans over too many cells
1138 #ifdef DEBUG
1139  for (unsigned int i=range.first+1; i<range.second; ++i)
1140  Assert (fe_indices[i] >= fe_indices[i-1],
1141  ExcMessage ("Cell range must be over sorted range of fe indices in hp case!"));
1142  AssertIndexRange(range.first,fe_indices.size()+1);
1143  AssertIndexRange(range.second,fe_indices.size()+1);
1144 #endif
1145  std::pair<unsigned int,unsigned int> return_range;
1146  return_range.first =
1147  std::lower_bound(&fe_indices[0] + range.first,
1148  &fe_indices[0] + range.second, fe_index)
1149  -&fe_indices[0] ;
1150  return_range.second =
1151  std::lower_bound(&fe_indices[0] + return_range.first,
1152  &fe_indices[0] + range.second,
1153  fe_index + 1)-&fe_indices[0];
1154  Assert(return_range.first >= range.first &&
1155  return_range.second <= range.second, ExcInternalError());
1156  return return_range;
1157  }
1158 }
1159 
1160 
1161 
1162 template <int dim, typename Number>
1163 inline
1164 void
1165 MatrixFree<dim,Number>::renumber_dofs (std::vector<types::global_dof_index> &renumbering,
1166  const unsigned int vector_component)
1167 {
1168  AssertIndexRange(vector_component, dof_info.size());
1169  dof_info[vector_component].renumber_dofs (renumbering);
1170 }
1171 
1172 
1173 
1174 template <int dim, typename Number>
1175 inline
1176 const DoFHandler<dim> &
1177 MatrixFree<dim,Number>::get_dof_handler (const unsigned int dof_index) const
1178 {
1179  AssertIndexRange (dof_index, n_components());
1180  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
1181  {
1182  AssertDimension (dof_handlers.dof_handler.size(),
1183  dof_handlers.n_dof_handlers);
1184  return *dof_handlers.dof_handler[dof_index];
1185  }
1186  else
1187  {
1188  Assert (false, ExcNotImplemented());
1189  // put pseudo return argument to avoid compiler error, but trigger a
1190  // segfault in case this is only run in optimized mode
1191  return *dof_handlers.dof_handler[numbers::invalid_unsigned_int];
1192  }
1193 }
1194 
1195 
1196 
1197 template <int dim, typename Number>
1198 inline
1200 MatrixFree<dim,Number>::get_cell_iterator(const unsigned int macro_cell_number,
1201  const unsigned int vector_number,
1202  const unsigned int dof_index) const
1203 {
1204  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
1205 #ifdef DEBUG
1206  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
1207  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
1208  AssertIndexRange (vector_number, vectorization_length);
1209  const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
1210  if (irreg_filled > 0)
1211  AssertIndexRange (vector_number, irreg_filled);
1212 #endif
1213 
1214  const DoFHandler<dim> *dofh = 0;
1215  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
1216  {
1217  AssertDimension (dof_handlers.dof_handler.size(),
1218  dof_handlers.n_dof_handlers);
1219  dofh = dof_handlers.dof_handler[dof_index];
1220  }
1221  else
1222  {
1223  Assert (false, ExcMessage ("Cannot return DoFHandler<dim>::cell_iterator "
1224  "for underlying DoFHandler!"));
1225  }
1226 
1227  std::pair<unsigned int,unsigned int> index =
1228  cell_level_index[macro_cell_number*vectorization_length+vector_number];
1229  return typename DoFHandler<dim>::cell_iterator
1230  (&dofh->get_triangulation(), index.first, index.second, dofh);
1231 }
1232 
1233 
1234 
1235 template <int dim, typename Number>
1236 inline
1238 MatrixFree<dim,Number>::get_hp_cell_iterator(const unsigned int macro_cell_number,
1239  const unsigned int vector_number,
1240  const unsigned int dof_index) const
1241 {
1242  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
1243 #ifdef DEBUG
1244  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
1245  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
1246  AssertIndexRange (vector_number, vectorization_length);
1247  const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
1248  if (irreg_filled > 0)
1249  AssertIndexRange (vector_number, irreg_filled);
1250 #endif
1251 
1252  Assert (dof_handlers.active_dof_handler == DoFHandlers::hp,
1253  ExcNotImplemented());
1254  const hp::DoFHandler<dim> *dofh = dof_handlers.hp_dof_handler[dof_index];
1255  std::pair<unsigned int,unsigned int> index =
1256  cell_level_index[macro_cell_number*vectorization_length+vector_number];
1257  return typename hp::DoFHandler<dim>::cell_iterator
1258  (&dofh->get_triangulation(), index.first, index.second, dofh);
1259 }
1260 
1261 
1262 
1263 template <int dim, typename Number>
1264 inline
1265 bool
1266 MatrixFree<dim,Number>::at_irregular_cell (const unsigned int macro_cell) const
1267 {
1268  AssertIndexRange (macro_cell, size_info.n_macro_cells);
1269  return dof_info[0].row_starts[macro_cell][2] > 0;
1270 }
1271 
1272 
1273 
1274 template <int dim, typename Number>
1275 inline
1276 unsigned int
1277 MatrixFree<dim,Number>::n_components_filled (const unsigned int macro_cell) const
1278 {
1279  AssertIndexRange (macro_cell, size_info.n_macro_cells);
1280  const unsigned int n_filled = dof_info[0].row_starts[macro_cell][2];
1281  if (n_filled == 0)
1283  else
1284  return n_filled;
1285 }
1286 
1287 
1288 
1289 template <int dim, typename Number>
1290 inline
1291 unsigned int
1292 MatrixFree<dim,Number>::get_dofs_per_cell(const unsigned int dof_index,
1293  const unsigned int active_fe_index) const
1294 {
1295  AssertIndexRange (dof_index, dof_info.size());
1296  return dof_info[dof_index].dofs_per_cell[active_fe_index];
1297 }
1298 
1299 
1300 
1301 template <int dim, typename Number>
1302 inline
1303 unsigned int
1304 MatrixFree<dim,Number>::get_n_q_points(const unsigned int quad_index,
1305  const unsigned int active_fe_index) const
1306 {
1307  AssertIndexRange (quad_index,
1308  mapping_info.mapping_data_gen.size());
1309  return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index];
1310 }
1311 
1312 
1313 
1314 template <int dim, typename Number>
1315 inline
1316 unsigned int
1317 MatrixFree<dim,Number>::get_dofs_per_face(const unsigned int dof_index,
1318  const unsigned int active_fe_index) const
1319 {
1320  AssertIndexRange (dof_index, dof_info.size());
1321  return dof_info[dof_index].dofs_per_face[active_fe_index];
1322 }
1323 
1324 
1325 
1326 template <int dim, typename Number>
1327 inline
1328 unsigned int
1329 MatrixFree<dim,Number>::get_n_q_points_face(const unsigned int quad_index,
1330  const unsigned int active_fe_index) const
1331 {
1332  AssertIndexRange (quad_index,
1333  mapping_info.mapping_data_gen.size());
1334  return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index];
1335 }
1336 
1337 
1338 
1339 template <int dim, typename Number>
1340 inline
1341 const IndexSet &
1342 MatrixFree<dim,Number>::get_locally_owned_set(const unsigned int dof_index) const
1343 {
1344  AssertIndexRange (dof_index, dof_info.size());
1345  return dof_info[dof_index].vector_partitioner->locally_owned_range();
1346 }
1347 
1348 
1349 
1350 template <int dim, typename Number>
1351 inline
1352 const IndexSet &
1353 MatrixFree<dim,Number>::get_ghost_set(const unsigned int dof_index) const
1354 {
1355  AssertIndexRange (dof_index, dof_info.size());
1356  return dof_info[dof_index].vector_partitioner->ghost_indices();
1357 }
1358 
1359 
1360 
1361 template <int dim, typename Number>
1362 inline
1364 MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
1365  const unsigned int index_quad,
1366  const unsigned int active_fe_index,
1367  const unsigned int active_quad_index) const
1368 {
1369  AssertIndexRange (index_fe, shape_info.size(0));
1370  AssertIndexRange (index_quad, shape_info.size(1));
1371  AssertIndexRange (active_fe_index, shape_info.size(2));
1372  AssertIndexRange (active_quad_index, shape_info.size(3));
1373  return shape_info(index_fe, index_quad,
1374  active_fe_index, active_quad_index);
1375 }
1376 
1377 
1378 
1379 template <int dim, typename Number>
1380 inline
1381 const Quadrature<dim> &
1382 MatrixFree<dim,Number>::get_quadrature (const unsigned int quad_index,
1383  const unsigned int active_fe_index) const
1384 {
1385  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
1386  return mapping_info.mapping_data_gen[quad_index].
1387  quadrature[active_fe_index];
1388 }
1389 
1390 
1391 
1392 template <int dim, typename Number>
1393 inline
1394 const Quadrature<dim-1> &
1395 MatrixFree<dim,Number>::get_face_quadrature (const unsigned int quad_index,
1396  const unsigned int active_fe_index) const
1397 {
1398  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
1399  return mapping_info.mapping_data_gen[quad_index].
1400  face_quadrature[active_fe_index];
1401 }
1402 
1403 
1404 
1405 template <int dim, typename Number>
1406 inline
1407 bool
1409 {
1410  return indices_are_initialized;
1411 }
1412 
1413 
1414 
1415 template <int dim, typename Number>
1416 inline
1417 bool
1419 {
1420  return mapping_is_initialized;
1421 }
1422 
1423 
1424 
1425 // ------------------------------ reinit functions ---------------------------
1426 
1427 namespace internal
1428 {
1429  namespace MatrixFree
1430  {
1431  template <typename DoFHandlerType>
1432  inline
1433  std::vector<IndexSet>
1434  extract_locally_owned_index_sets (const std::vector<const DoFHandlerType *> &dofh,
1435  const unsigned int level)
1436  {
1437  std::vector<IndexSet> locally_owned_set;
1438  locally_owned_set.reserve (dofh.size());
1439  for (unsigned int j=0; j<dofh.size(); j++)
1440  if (level == numbers::invalid_unsigned_int)
1441  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
1442  else
1443  AssertThrow(false, ExcNotImplemented());
1444  return locally_owned_set;
1445  }
1446 
1447  template <int dim, int spacedim>
1448  inline
1449  std::vector<IndexSet>
1450  extract_locally_owned_index_sets (const std::vector<const ::DoFHandler<dim,spacedim> *> &dofh,
1451  const unsigned int level)
1452  {
1453  std::vector<IndexSet> locally_owned_set;
1454  locally_owned_set.reserve (dofh.size());
1455  for (unsigned int j=0; j<dofh.size(); j++)
1456  if (level == numbers::invalid_unsigned_int)
1457  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
1458  else
1459  locally_owned_set.push_back(dofh[j]->locally_owned_mg_dofs(level));
1460  return locally_owned_set;
1461  }
1462  }
1463 }
1464 
1465 
1466 
1467 template <int dim, typename Number>
1468 template <typename DoFHandlerType, typename QuadratureType>
1470 reinit(const DoFHandlerType &dof_handler,
1471  const ConstraintMatrix &constraints_in,
1472  const QuadratureType &quad,
1473  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1474 {
1475  std::vector<const DoFHandlerType *> dof_handlers;
1476  std::vector<const ConstraintMatrix *> constraints;
1477  std::vector<QuadratureType> quads;
1478 
1479  dof_handlers.push_back(&dof_handler);
1480  constraints.push_back (&constraints_in);
1481  quads.push_back (quad);
1482 
1483  std::vector<IndexSet> locally_owned_sets =
1484  internal::MatrixFree::extract_locally_owned_index_sets
1485  (dof_handlers, additional_data.level_mg_handler);
1486  reinit(StaticMappingQ1<dim>::mapping, dof_handlers,constraints, locally_owned_sets, quads,
1487  additional_data);
1488 }
1489 
1490 
1491 
1492 template <int dim, typename Number>
1493 template <typename DoFHandlerType, typename QuadratureType>
1495 reinit(const Mapping<dim> &mapping,
1496  const DoFHandlerType &dof_handler,
1497  const ConstraintMatrix &constraints_in,
1498  const QuadratureType &quad,
1499  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1500 {
1501  std::vector<const DoFHandlerType *> dof_handlers;
1502  std::vector<const ConstraintMatrix *> constraints;
1503  std::vector<QuadratureType> quads;
1504 
1505  dof_handlers.push_back(&dof_handler);
1506  constraints.push_back (&constraints_in);
1507  quads.push_back (quad);
1508 
1509  std::vector<IndexSet> locally_owned_sets =
1510  internal::MatrixFree::extract_locally_owned_index_sets
1511  (dof_handlers, additional_data.level_mg_handler);
1512  reinit(mapping, dof_handlers,constraints,locally_owned_sets, quads,
1513  additional_data);
1514 }
1515 
1516 
1517 
1518 template <int dim, typename Number>
1519 template <typename DoFHandlerType, typename QuadratureType>
1521 reinit(const std::vector<const DoFHandlerType *> &dof_handler,
1522  const std::vector<const ConstraintMatrix *> &constraint,
1523  const std::vector<QuadratureType> &quad,
1524  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1525 {
1526  std::vector<IndexSet> locally_owned_set =
1527  internal::MatrixFree::extract_locally_owned_index_sets
1528  (dof_handler, additional_data.level_mg_handler);
1529  reinit(StaticMappingQ1<dim>::mapping, dof_handler,constraint,locally_owned_set,
1530  static_cast<const std::vector<Quadrature<1> >&>(quad),
1531  additional_data);
1532 }
1533 
1534 
1535 
1536 template <int dim, typename Number>
1537 template <typename DoFHandlerType, typename QuadratureType>
1539 reinit(const std::vector<const DoFHandlerType *> &dof_handler,
1540  const std::vector<const ConstraintMatrix *> &constraint,
1541  const QuadratureType &quad,
1542  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1543 {
1544  std::vector<QuadratureType> quads;
1545  quads.push_back(quad);
1546  std::vector<IndexSet> locally_owned_set =
1547  internal::MatrixFree::extract_locally_owned_index_sets
1548  (dof_handler, additional_data.level_mg_handler);
1549  reinit(StaticMappingQ1<dim>::mapping, dof_handler,constraint,locally_owned_set, quads,
1550  additional_data);
1551 }
1552 
1553 
1554 
1555 template <int dim, typename Number>
1556 template <typename DoFHandlerType, typename QuadratureType>
1558 reinit(const Mapping<dim> &mapping,
1559  const std::vector<const DoFHandlerType *> &dof_handler,
1560  const std::vector<const ConstraintMatrix *> &constraint,
1561  const QuadratureType &quad,
1562  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1563 {
1564  std::vector<QuadratureType> quads;
1565  quads.push_back(quad);
1566  std::vector<IndexSet> locally_owned_set =
1567  internal::MatrixFree::extract_locally_owned_index_sets
1568  (dof_handler, additional_data.level_mg_handler);
1569  reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
1570  additional_data);
1571 }
1572 
1573 
1574 
1575 template <int dim, typename Number>
1576 template <typename DoFHandlerType, typename QuadratureType>
1578 reinit(const Mapping<dim> &mapping,
1579  const std::vector<const DoFHandlerType *> &dof_handler,
1580  const std::vector<const ConstraintMatrix *> &constraint,
1581  const std::vector<QuadratureType> &quad,
1582  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1583 {
1584  std::vector<IndexSet> locally_owned_set =
1585  internal::MatrixFree::extract_locally_owned_index_sets
1586  (dof_handler, additional_data.level_mg_handler);
1587  reinit(mapping, dof_handler,constraint,locally_owned_set,
1588  quad, additional_data);
1589 }
1590 
1591 
1592 
1593 template <int dim, typename Number>
1594 template <typename DoFHandlerType, typename QuadratureType>
1596 reinit(const Mapping<dim> &mapping,
1597  const std::vector<const DoFHandlerType *> &dof_handler,
1598  const std::vector<const ConstraintMatrix *> &constraint,
1599  const std::vector<IndexSet> &locally_owned_set,
1600  const std::vector<QuadratureType> &quad,
1601  const typename MatrixFree<dim,Number>::AdditionalData additional_data)
1602 {
1603  // find out whether we use a hp Quadrature or a standard quadrature
1604  std::vector<hp::QCollection<1> > quad_hp;
1605  for (unsigned int q=0; q<quad.size(); ++q)
1606  quad_hp.push_back (hp::QCollection<1>(quad[q]));
1607  internal_reinit (mapping,
1608  dof_handler,
1609  constraint, locally_owned_set, quad_hp, additional_data);
1610 }
1611 
1612 
1613 
1614 // ------------------------------ implementation of cell_loop ---------------
1615 
1616 // internal helper functions that define how to call MPI data exchange
1617 // functions: for generic vectors, do nothing at all. For distributed vectors,
1618 // call update_ghost_values_start function and so on. If we have collections
1619 // of vectors, just do the individual functions of the components. In order to
1620 // keep ghost values consistent (whether we are in read or write mode). the whole situation is a bit complicated by the fact
1621 // that we need to treat block vectors differently, which use some additional
1622 // helper functions to select the blocks and template magic.
1623 namespace internal
1624 {
1625  template<typename VectorStruct>
1626  bool update_ghost_values_start_block (const VectorStruct &vec,
1627  const unsigned int channel,
1629  template<typename VectorStruct>
1630  void reset_ghost_values_block (const VectorStruct &vec,
1631  const bool zero_out_ghosts,
1633  template<typename VectorStruct>
1634  void update_ghost_values_finish_block (const VectorStruct &vec,
1636  template<typename VectorStruct>
1637  void compress_start_block (const VectorStruct &vec,
1638  const unsigned int channel,
1640  template<typename VectorStruct>
1641  void compress_finish_block (const VectorStruct &vec,
1643 
1644  template<typename VectorStruct>
1645  bool update_ghost_values_start_block (const VectorStruct &,
1646  const unsigned int,
1648  {
1649  return false;
1650  }
1651  template<typename VectorStruct>
1652  void reset_ghost_values_block (const VectorStruct &,
1653  const bool,
1655  {}
1656  template<typename VectorStruct>
1657  void update_ghost_values_finish_block (const VectorStruct &,
1659  {}
1660  template<typename VectorStruct>
1661  void compress_start_block (const VectorStruct &,
1662  const unsigned int,
1664  {}
1665  template<typename VectorStruct>
1666  void compress_finish_block (const VectorStruct &,
1668  {}
1669 
1670 
1671 
1672  // returns true if the vector was in a state without ghost values before,
1673  // i.e., we need to zero out ghosts in the very end
1674  template<typename VectorStruct>
1675  inline
1676  bool update_ghost_values_start (const VectorStruct &vec,
1677  const unsigned int channel = 0)
1678  {
1679  return
1680  update_ghost_values_start_block(vec, channel,
1682  }
1683 
1684 
1685 
1686  template<typename Number>
1687  inline
1688  bool update_ghost_values_start (const parallel::distributed::Vector<Number> &vec,
1689  const unsigned int channel = 0)
1690  {
1691  bool return_value = !vec.has_ghost_elements();
1692  vec.update_ghost_values_start(channel);
1693  return return_value;
1694  }
1695 
1696 
1697 
1698  template <typename VectorStruct>
1699  inline
1700  bool update_ghost_values_start (const std::vector<VectorStruct> &vec)
1701  {
1702  bool return_value = false;
1703  for (unsigned int comp=0; comp<vec.size(); comp++)
1704  return_value = update_ghost_values_start(vec[comp], comp);
1705  return return_value;
1706  }
1707 
1708 
1709 
1710  template <typename VectorStruct>
1711  inline
1712  bool update_ghost_values_start (const std::vector<VectorStruct *> &vec)
1713  {
1714  bool return_value = false;
1715  for (unsigned int comp=0; comp<vec.size(); comp++)
1716  return_value = update_ghost_values_start(*vec[comp], comp);
1717  return return_value;
1718  }
1719 
1720 
1721 
1722  template<typename VectorStruct>
1723  inline
1724  bool update_ghost_values_start_block (const VectorStruct &vec,
1725  const unsigned int channel,
1727  {
1728  bool return_value = false;
1729  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1730  return_value = update_ghost_values_start(vec.block(i), channel+509*i);
1731  return return_value;
1732  }
1733 
1734 
1735 
1736  // if the input vector did not have ghosts imported, clear them here again
1737  // in order to avoid subsequent operations e.g. in linear solvers to work
1738  // with ghosts all the time
1739  template<typename VectorStruct>
1740  inline
1741  void reset_ghost_values (const VectorStruct &vec,
1742  const bool zero_out_ghosts)
1743  {
1744  reset_ghost_values_block(vec, zero_out_ghosts,
1746  }
1747 
1748 
1749 
1750  template<typename Number>
1751  inline
1752  void reset_ghost_values (const parallel::distributed::Vector<Number> &vec,
1753  const bool zero_out_ghosts)
1754  {
1755  if (zero_out_ghosts)
1756  const_cast<parallel::distributed::Vector<Number>&>(vec).zero_out_ghosts();
1757  }
1758 
1759 
1760 
1761  template <typename VectorStruct>
1762  inline
1763  void reset_ghost_values (const std::vector<VectorStruct> &vec,
1764  const bool zero_out_ghosts)
1765  {
1766  for (unsigned int comp=0; comp<vec.size(); comp++)
1767  reset_ghost_values(vec[comp], zero_out_ghosts);
1768  }
1769 
1770 
1771 
1772  template <typename VectorStruct>
1773  inline
1774  void reset_ghost_values (const std::vector<VectorStruct *> &vec,
1775  const bool zero_out_ghosts)
1776  {
1777  for (unsigned int comp=0; comp<vec.size(); comp++)
1778  reset_ghost_values(*vec[comp], zero_out_ghosts);
1779  }
1780 
1781 
1782 
1783  template<typename VectorStruct>
1784  inline
1785  void reset_ghost_values_block (const VectorStruct &vec,
1786  const bool zero_out_ghosts,
1788  {
1789  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1790  reset_ghost_values(vec.block(i), zero_out_ghosts);
1791  }
1792 
1793 
1794 
1795  template <typename VectorStruct>
1796  inline
1797  void update_ghost_values_finish (const VectorStruct &vec)
1798  {
1799  update_ghost_values_finish_block(vec,
1801  }
1802 
1803 
1804 
1805  template <typename Number>
1806  inline
1807  void update_ghost_values_finish (const parallel::distributed::Vector<Number> &vec)
1808  {
1810  }
1811 
1812 
1813 
1814  template <typename VectorStruct>
1815  inline
1816  void update_ghost_values_finish (const std::vector<VectorStruct> &vec)
1817  {
1818  for (unsigned int comp=0; comp<vec.size(); comp++)
1819  update_ghost_values_finish(vec[comp]);
1820  }
1821 
1822 
1823 
1824  template <typename VectorStruct>
1825  inline
1826  void update_ghost_values_finish (const std::vector<VectorStruct *> &vec)
1827  {
1828  for (unsigned int comp=0; comp<vec.size(); comp++)
1829  update_ghost_values_finish(*vec[comp]);
1830  }
1831 
1832 
1833 
1834  template <typename VectorStruct>
1835  inline
1836  void update_ghost_values_finish_block (const VectorStruct &vec,
1838  {
1839  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1840  update_ghost_values_finish(vec.block(i));
1841  }
1842 
1843 
1844 
1845  template <typename VectorStruct>
1846  inline
1847  void compress_start (VectorStruct &vec,
1848  const unsigned int channel = 0)
1849  {
1850  compress_start_block (vec, channel,
1852  }
1853 
1854 
1855 
1856  template <typename Number>
1857  inline
1858  void compress_start (parallel::distributed::Vector<Number> &vec,
1859  const unsigned int channel = 0)
1860  {
1861  vec.compress_start(channel);
1862  }
1863 
1864 
1865 
1866  template <typename VectorStruct>
1867  inline
1868  void compress_start (std::vector<VectorStruct> &vec)
1869  {
1870  for (unsigned int comp=0; comp<vec.size(); comp++)
1871  compress_start (vec[comp], comp);
1872  }
1873 
1874 
1875 
1876  template <typename VectorStruct>
1877  inline
1878  void compress_start (std::vector<VectorStruct *> &vec)
1879  {
1880  for (unsigned int comp=0; comp<vec.size(); comp++)
1881  compress_start (*vec[comp], comp);
1882  }
1883 
1884 
1885 
1886  template <typename VectorStruct>
1887  inline
1888  void compress_start_block (VectorStruct &vec,
1889  const unsigned int channel,
1891  {
1892  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1893  compress_start(vec.block(i), channel + 500*i);
1894  }
1895 
1896 
1897 
1898  template <typename VectorStruct>
1899  inline
1900  void compress_finish (VectorStruct &vec)
1901  {
1902  compress_finish_block(vec,
1904  }
1905 
1906 
1907 
1908  template <typename Number>
1909  inline
1910  void compress_finish (parallel::distributed::Vector<Number> &vec)
1911  {
1912  vec.compress_finish(::VectorOperation::add);
1913  }
1914 
1915 
1916 
1917  template <typename VectorStruct>
1918  inline
1919  void compress_finish (std::vector<VectorStruct> &vec)
1920  {
1921  for (unsigned int comp=0; comp<vec.size(); comp++)
1922  compress_finish(vec[comp]);
1923  }
1924 
1925 
1926 
1927  template <typename VectorStruct>
1928  inline
1929  void compress_finish (std::vector<VectorStruct *> &vec)
1930  {
1931  for (unsigned int comp=0; comp<vec.size(); comp++)
1932  compress_finish(*vec[comp]);
1933  }
1934 
1935 
1936 
1937  template <typename VectorStruct>
1938  inline
1939  void compress_finish_block (VectorStruct &vec,
1941  {
1942  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1943  compress_finish(vec.block(i));
1944  }
1945 
1946 
1947 
1948 #ifdef DEAL_II_WITH_THREADS
1949 
1950  // This defines the TBB data structures that are needed to schedule the
1951  // partition-partition variant
1952 
1953  namespace partition
1954  {
1955  template<typename Worker>
1956  class CellWork : public tbb::task
1957  {
1958  public:
1959  CellWork (const Worker &worker_in,
1960  const unsigned int partition_in,
1961  const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
1962  const bool is_blocked_in)
1963  :
1964  worker (worker_in),
1965  partition (partition_in),
1966  task_info (task_info_in),
1967  is_blocked (is_blocked_in)
1968  {};
1969  tbb::task *execute ()
1970  {
1971  std::pair<unsigned int, unsigned int> cell_range
1972  (task_info.partition_color_blocks_data[partition],
1973  task_info.partition_color_blocks_data[partition+1]);
1974  worker(cell_range);
1975  if (is_blocked==true)
1976  dummy->spawn (*dummy);
1977  return NULL;
1978  }
1979 
1980  tbb::empty_task *dummy;
1981 
1982  private:
1983  const Worker &worker;
1984  const unsigned int partition;
1985  const internal::MatrixFreeFunctions::TaskInfo &task_info;
1986  const bool is_blocked;
1987  };
1988 
1989 
1990 
1991  template<typename Worker>
1992  class PartitionWork : public tbb::task
1993  {
1994  public:
1995  PartitionWork (const Worker &function_in,
1996  const unsigned int partition_in,
1997  const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
1998  const bool is_blocked_in = false)
1999  :
2000  function (function_in),
2001  partition (partition_in),
2002  task_info (task_info_in),
2003  is_blocked (is_blocked_in)
2004  {};
2005  tbb::task *execute ()
2006  {
2007  tbb::empty_task *root = new( tbb::task::allocate_root() )
2008  tbb::empty_task;
2009  unsigned int evens = task_info.partition_evens[partition];
2010  unsigned int odds = task_info.partition_odds[partition];
2011  unsigned int n_blocked_workers =
2012  task_info.partition_n_blocked_workers[partition];
2013  unsigned int n_workers = task_info.partition_n_workers[partition];
2014  std::vector<CellWork<Worker>*> worker(n_workers);
2015  std::vector<CellWork<Worker>*> blocked_worker(n_blocked_workers);
2016 
2017  root->set_ref_count(evens+1);
2018  for (unsigned int j=0; j<evens; j++)
2019  {
2020  worker[j] = new(root->allocate_child())
2021  CellWork<Worker>(function, task_info.
2022  partition_color_blocks_row_index[partition]+2*j,
2023  task_info, false);
2024  if (j>0)
2025  {
2026  worker[j]->set_ref_count(2);
2027  blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
2028  tbb::empty_task;
2029  worker[j-1]->spawn(*blocked_worker[j-1]);
2030  }
2031  else
2032  worker[j]->set_ref_count(1);
2033  if (j<evens-1)
2034  {
2035  blocked_worker[j] = new(worker[j]->allocate_child())
2036  CellWork<Worker>(function, task_info.
2037  partition_color_blocks_row_index
2038  [partition] + 2*j+1, task_info, true);
2039  }
2040  else
2041  {
2042  if (odds==evens)
2043  {
2044  worker[evens] = new(worker[j]->allocate_child())
2045  CellWork<Worker>(function, task_info.
2046  partition_color_blocks_row_index[partition]+2*j+1,
2047  task_info, false);
2048  worker[j]->spawn(*worker[evens]);
2049  }
2050  else
2051  {
2052  tbb::empty_task *child = new(worker[j]->allocate_child())
2053  tbb::empty_task();
2054  worker[j]->spawn(*child);
2055  }
2056  }
2057  }
2058 
2059  root->wait_for_all();
2060  root->destroy(*root);
2061  if (is_blocked==true)
2062  dummy->spawn (*dummy);
2063  return NULL;
2064  }
2065 
2066  tbb::empty_task *dummy;
2067 
2068  private:
2069  const Worker &function;
2070  const unsigned int partition;
2071  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2072  const bool is_blocked;
2073  };
2074 
2075  } // end of namespace partition
2076 
2077 
2078 
2079  namespace color
2080  {
2081  template <typename Worker>
2082  class CellWork
2083  {
2084  public:
2085  CellWork (const Worker &worker_in,
2086  const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
2087  :
2088  worker (worker_in),
2089  task_info (task_info_in)
2090  {};
2091  void operator()(const tbb::blocked_range<unsigned int> &r) const
2092  {
2093  for (unsigned int block=r.begin(); block<r.end(); block++)
2094  {
2095  std::pair<unsigned int,unsigned int> cell_range;
2096  if (task_info.position_short_block<block)
2097  {
2098  cell_range.first = (block-1)*task_info.block_size+
2099  task_info.block_size_last;
2100  cell_range.second = cell_range.first + task_info.block_size;
2101  }
2102  else
2103  {
2104  cell_range.first = block*task_info.block_size;
2105  cell_range.second = cell_range.first +
2106  ((block == task_info.position_short_block)?
2107  (task_info.block_size_last):(task_info.block_size));
2108  }
2109  worker (cell_range);
2110  }
2111  }
2112  private:
2113  const Worker &worker;
2114  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2115  };
2116 
2117 
2118  template<typename Worker>
2119  class PartitionWork : public tbb::task
2120  {
2121  public:
2122  PartitionWork (const Worker &worker_in,
2123  const unsigned int partition_in,
2124  const internal::MatrixFreeFunctions::TaskInfo &task_info_in,
2125  const bool is_blocked_in)
2126  :
2127  worker (worker_in),
2128  partition (partition_in),
2129  task_info (task_info_in),
2130  is_blocked (is_blocked_in)
2131  {};
2132  tbb::task *execute ()
2133  {
2134  unsigned int lower = task_info.partition_color_blocks_data[partition],
2135  upper = task_info.partition_color_blocks_data[partition+1];
2136  parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
2137  CellWork<Worker> (worker,task_info));
2138  if (is_blocked==true)
2139  dummy->spawn (*dummy);
2140  return NULL;
2141  }
2142 
2143  tbb::empty_task *dummy;
2144 
2145  private:
2146  const Worker &worker;
2147  const unsigned int partition;
2148  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2149  const bool is_blocked;
2150  };
2151 
2152  } // end of namespace color
2153 
2154 
2155  template<typename VectorStruct>
2156  class MPIComDistribute : public tbb::task
2157  {
2158  public:
2159  MPIComDistribute (const VectorStruct &src_in)
2160  :
2161  src(src_in)
2162  {};
2163 
2164  tbb::task *execute ()
2165  {
2166  internal::update_ghost_values_finish(src);
2167  return 0;
2168  }
2169 
2170  private:
2171  const VectorStruct &src;
2172  };
2173 
2174 
2175 
2176  template<typename VectorStruct>
2177  class MPIComCompress : public tbb::task
2178  {
2179  public:
2180  MPIComCompress (VectorStruct &dst_in)
2181  :
2182  dst(dst_in)
2183  {};
2184 
2185  tbb::task *execute ()
2186  {
2187  internal::compress_start(dst);
2188  return 0;
2189  }
2190 
2191  private:
2192  VectorStruct &dst;
2193  };
2194 
2195 #endif // DEAL_II_WITH_THREADS
2196 
2197 } // end of namespace internal
2198 
2199 
2200 
2201 template <int dim, typename Number>
2202 template <typename OutVector, typename InVector>
2203 inline
2204 void
2206 (const std_cxx11::function<void (const MatrixFree<dim,Number> &,
2207  OutVector &,
2208  const InVector &,
2209  const std::pair<unsigned int,
2210  unsigned int> &)> &cell_operation,
2211  OutVector &dst,
2212  const InVector &src) const
2213 {
2214  // in any case, need to start the ghost import at the beginning
2215  bool ghosts_were_not_set = internal::update_ghost_values_start (src);
2216 
2217 #ifdef DEAL_II_WITH_THREADS
2218 
2219  // Use multithreading if so requested and if there is enough work to do in
2220  // parallel (the code might hang if there are less than two chunks!)
2221  if (task_info.use_multithreading == true && task_info.n_blocks > 3)
2222  {
2223  // to simplify the function calls, bind away all arguments except the
2224  // cell range
2225  typedef
2226  std_cxx11::function<void (const std::pair<unsigned int,unsigned int> &range)>
2227  Worker;
2228 
2229  const Worker func = std_cxx11::bind (std_cxx11::ref(cell_operation),
2230  std_cxx11::cref(*this),
2231  std_cxx11::ref(dst),
2232  std_cxx11::cref(src),
2233  std_cxx11::_1);
2234 
2235  if (task_info.use_partition_partition == true)
2236  {
2237  tbb::empty_task *root = new( tbb::task::allocate_root() )
2238  tbb::empty_task;
2239  unsigned int evens = task_info.evens;
2240  unsigned int odds = task_info.odds;
2241  root->set_ref_count(evens+1);
2242  unsigned int n_blocked_workers = task_info.n_blocked_workers;
2243  unsigned int n_workers = task_info.n_workers;
2244  std::vector<internal::partition::PartitionWork<Worker>*>
2245  worker(n_workers);
2246  std::vector<internal::partition::PartitionWork<Worker>*>
2247  blocked_worker(n_blocked_workers);
2248  internal::MPIComCompress<OutVector> *worker_compr =
2249  new(root->allocate_child())
2250  internal::MPIComCompress<OutVector>(dst);
2251  worker_compr->set_ref_count(1);
2252  for (unsigned int j=0; j<evens; j++)
2253  {
2254  if (j>0)
2255  {
2256  worker[j] = new(root->allocate_child())
2257  internal::partition::PartitionWork<Worker>
2258  (func,2*j,task_info,false);
2259  worker[j]->set_ref_count(2);
2260  blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
2261  tbb::empty_task;
2262  if (j>1)
2263  worker[j-1]->spawn(*blocked_worker[j-1]);
2264  else
2265  worker_compr->spawn(*blocked_worker[j-1]);
2266  }
2267  else
2268  {
2269  worker[j] = new(worker_compr->allocate_child())
2270  internal::partition::PartitionWork<Worker>
2271  (func,2*j,task_info,false);
2272  worker[j]->set_ref_count(2);
2273  internal::MPIComDistribute<InVector> *worker_dist =
2274  new (worker[j]->allocate_child())
2275  internal::MPIComDistribute<InVector>(src);
2276  worker_dist->spawn(*worker_dist);
2277  }
2278  if (j<evens-1)
2279  {
2280  blocked_worker[j] = new(worker[j]->allocate_child())
2281  internal::partition::PartitionWork<Worker>
2282  (func,2*j+1,task_info,true);
2283  }
2284  else
2285  {
2286  if (odds==evens)
2287  {
2288  worker[evens] = new(worker[j]->allocate_child())
2289  internal::partition::PartitionWork<Worker>
2290  (func,2*j+1,task_info,false);
2291  worker[j]->spawn(*worker[evens]);
2292  }
2293  else
2294  {
2295  tbb::empty_task *child = new(worker[j]->allocate_child())
2296  tbb::empty_task();
2297  worker[j]->spawn(*child);
2298  }
2299  }
2300  }
2301 
2302  root->wait_for_all();
2303  root->destroy(*root);
2304  }
2305  else // end of partition-partition, start of partition-color
2306  {
2307  unsigned int evens = task_info.evens;
2308  unsigned int odds = task_info.odds;
2309 
2310  // check whether there is only one partition. if not, build up the
2311  // tree of partitions
2312  if (odds > 0)
2313  {
2314  tbb::empty_task *root = new( tbb::task::allocate_root() ) tbb::empty_task;
2315  root->set_ref_count(evens+1);
2316  unsigned int n_blocked_workers = odds-(odds+evens+1)%2;
2317  unsigned int n_workers = task_info.partition_color_blocks_data.size()-1-
2318  n_blocked_workers;
2319  std::vector<internal::color::PartitionWork<Worker>*> worker(n_workers);
2320  std::vector<internal::color::PartitionWork<Worker>*> blocked_worker(n_blocked_workers);
2321  unsigned int worker_index = 0, slice_index = 0;
2322  unsigned int spawn_index = 0, spawn_index_new = 0;
2323  int spawn_index_child = -2;
2324  internal::MPIComCompress<OutVector> *worker_compr = new(root->allocate_child())
2325  internal::MPIComCompress<OutVector>(dst);
2326  worker_compr->set_ref_count(1);
2327  for (unsigned int part=0;
2328  part<task_info.partition_color_blocks_row_index.size()-1; part++)
2329  {
2330  spawn_index_new = worker_index;
2331  if (part == 0)
2332  worker[worker_index] = new(worker_compr->allocate_child())
2333  internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
2334  else
2335  worker[worker_index] = new(root->allocate_child())
2336  internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
2337  slice_index++;
2338  for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
2339  slice_index++)
2340  {
2341  worker[worker_index]->set_ref_count(1);
2342  worker_index++;
2343  worker[worker_index] = new (worker[worker_index-1]->allocate_child())
2344  internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
2345  }
2346  worker[worker_index]->set_ref_count(2);
2347  if (part>0)
2348  {
2349  blocked_worker[(part-1)/2]->dummy =
2350  new (worker[worker_index]->allocate_child()) tbb::empty_task;
2351  worker_index++;
2352  if (spawn_index_child == -1)
2353  worker[spawn_index]->spawn(*blocked_worker[(part-1)/2]);
2354  else
2355  worker[spawn_index]->spawn(*worker[spawn_index_child]);
2356  spawn_index = spawn_index_new;
2357  spawn_index_child = -2;
2358  }
2359  else
2360  {
2361  internal::MPIComDistribute<InVector> *worker_dist =
2362  new (worker[worker_index]->allocate_child())
2363  internal::MPIComDistribute<InVector>(src);
2364  worker_dist->spawn(*worker_dist);
2365  worker_index++;
2366  }
2367  part += 1;
2368  if (part<task_info.partition_color_blocks_row_index.size()-1)
2369  {
2370  if (part<task_info.partition_color_blocks_row_index.size()-2)
2371  {
2372  blocked_worker[part/2] = new(worker[worker_index-1]->allocate_child())
2373  internal::color::PartitionWork<Worker>(func,slice_index,task_info,true);
2374  slice_index++;
2375  if (slice_index<
2376  task_info.partition_color_blocks_row_index[part+1])
2377  {
2378  blocked_worker[part/2]->set_ref_count(1);
2379  worker[worker_index] = new(blocked_worker[part/2]->allocate_child())
2380  internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
2381  slice_index++;
2382  }
2383  else
2384  {
2385  spawn_index_child = -1;
2386  continue;
2387  }
2388  }
2389  for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
2390  slice_index++)
2391  {
2392  if (slice_index>
2393  task_info.partition_color_blocks_row_index[part])
2394  {
2395  worker[worker_index]->set_ref_count(1);
2396  worker_index++;
2397  }
2398  worker[worker_index] = new (worker[worker_index-1]->allocate_child())
2399  internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
2400  }
2401  spawn_index_child = worker_index;
2402  worker_index++;
2403  }
2404  else
2405  {
2406  tbb::empty_task *final = new (worker[worker_index-1]->allocate_child())
2407  tbb::empty_task;
2408  worker[spawn_index]->spawn(*final);
2409  spawn_index_child = worker_index-1;
2410  }
2411  }
2412  if (evens==odds)
2413  worker[spawn_index]->spawn(*worker[spawn_index_child]);
2414  root->wait_for_all();
2415  root->destroy(*root);
2416  }
2417  // case when we only have one partition: this is the usual coloring
2418  // scheme, and we just schedule a parallel for loop for each color
2419  else
2420  {
2421  Assert(evens==1,ExcInternalError());
2422  internal::update_ghost_values_finish(src);
2423 
2424  for (unsigned int color=0;
2425  color < task_info.partition_color_blocks_row_index[1];
2426  ++color)
2427  {
2428  unsigned int lower = task_info.partition_color_blocks_data[color],
2429  upper = task_info.partition_color_blocks_data[color+1];
2430  parallel_for(tbb::blocked_range<unsigned int>(lower,upper,1),
2431  internal::color::CellWork<Worker>
2432  (func,task_info));
2433  }
2434 
2435  internal::compress_start(dst);
2436  }
2437  }
2438  }
2439  else
2440 #endif
2441  // serial loop
2442  {
2443  std::pair<unsigned int,unsigned int> cell_range;
2444 
2445  // First operate on cells where no ghost data is needed (inner cells)
2446  {
2447  cell_range.first = 0;
2448  cell_range.second = size_info.boundary_cells_start;
2449  cell_operation (*this, dst, src, cell_range);
2450  }
2451 
2452  // before starting operations on cells that contain ghost nodes (outer
2453  // cells), wait for the MPI commands to finish
2454  internal::update_ghost_values_finish(src);
2455 
2456  // For the outer cells, do the same procedure as for inner cells.
2457  if (size_info.boundary_cells_end > size_info.boundary_cells_start)
2458  {
2459  cell_range.first = size_info.boundary_cells_start;
2460  cell_range.second = size_info.boundary_cells_end;
2461  cell_operation (*this, dst, src, cell_range);
2462  }
2463 
2464  internal::compress_start(dst);
2465 
2466  // Finally operate on cells where no ghost data is needed (inner cells)
2467  if (size_info.n_macro_cells > size_info.boundary_cells_end)
2468  {
2469  cell_range.first = size_info.boundary_cells_end;
2470  cell_range.second = size_info.n_macro_cells;
2471  cell_operation (*this, dst, src, cell_range);
2472  }
2473  }
2474 
2475  // In every case, we need to finish transfers at the very end
2476  internal::compress_finish(dst);
2477  internal::reset_ghost_values(src, ghosts_were_not_set);
2478 }
2479 
2480 
2481 
2482 template <int dim, typename Number>
2483 template <typename CLASS, typename OutVector, typename InVector>
2484 inline
2485 void
2487 (void (CLASS::*function_pointer)(const MatrixFree<dim,Number> &,
2488  OutVector &,
2489  const InVector &,
2490  const std::pair<unsigned int,
2491  unsigned int> &)const,
2492  const CLASS *owning_class,
2493  OutVector &dst,
2494  const InVector &src) const
2495 {
2496  // here, use std_cxx11::bind to hand a function handler with the appropriate
2497  // argument to the other loop function
2498  std_cxx11::function<void (const MatrixFree<dim,Number> &,
2499  OutVector &,
2500  const InVector &,
2501  const std::pair<unsigned int,
2502  unsigned int> &)>
2503  function = std_cxx11::bind<void>(function_pointer,
2504  owning_class,
2505  std_cxx11::_1,
2506  std_cxx11::_2,
2507  std_cxx11::_3,
2508  std_cxx11::_4);
2509  cell_loop (function, dst, src);
2510 }
2511 
2512 
2513 
2514 template <int dim, typename Number>
2515 template <typename CLASS, typename OutVector, typename InVector>
2516 inline
2517 void
2519 (void(CLASS::*function_pointer)(const MatrixFree<dim,Number> &,
2520  OutVector &,
2521  const InVector &,
2522  const std::pair<unsigned int,
2523  unsigned int> &),
2524  CLASS *owning_class,
2525  OutVector &dst,
2526  const InVector &src) const
2527 {
2528  // here, use std_cxx11::bind to hand a function handler with the appropriate
2529  // argument to the other loop function
2530  std_cxx11::function<void (const MatrixFree<dim,Number> &,
2531  OutVector &,
2532  const InVector &,
2533  const std::pair<unsigned int,
2534  unsigned int> &)>
2535  function = std_cxx11::bind<void>(function_pointer,
2536  owning_class,
2537  std_cxx11::_1,
2538  std_cxx11::_2,
2539  std_cxx11::_3,
2540  std_cxx11::_4);
2541  cell_loop (function, dst, src);
2542 }
2543 
2544 
2545 #endif // ifndef DOXYGEN
2546 
2547 
2548 
2549 DEAL_II_NAMESPACE_CLOSE
2550 
2551 #endif
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)
void reinit(const Mapping< dim > &mapping, const DoFHandlerType &dof_handler, const ConstraintMatrix &constraint, const IndexSet &locally_owned_dofs, const QuadratureType &quad, const AdditionalData additional_data=AdditionalData())
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
Definition: matrix_free.h:919
const IndexSet & get_ghost_set(const unsigned int fe_component=0) const
bool indices_are_initialized
Definition: matrix_free.h:924
static const unsigned int invalid_unsigned_int
Definition: types.h:164
bool at_irregular_cell(const unsigned int macro_cell_number) const
const Triangulation< dim, spacedim > & get_triangulation() const
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:889
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1052
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)
Definition: matrix_free.h:146
internal::MatrixFreeFunctions::MappingInfo< dim, Number > mapping_info
Definition: matrix_free.h:895
const DoFHandler< dim > & get_dof_handler(const unsigned int fe_component=0) const
DoFHandlers dof_handlers
Definition: matrix_free.h:869
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)
::ExceptionBase & ExcMessage(std::string arg1)
std::vector< Number > constraint_pool_data
Definition: matrix_free.h:883
#define AssertIndexRange(index, range)
Definition: exceptions.h:1081
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
Definition: matrix_free.h:900
void copy_from(const MatrixFree< dim, Number > &matrix_free_base)
const Triangulation< dim, spacedim > & get_triangulation() const
unsigned int n_constraint_pool_entries() const
std::vector< std::pair< unsigned int, unsigned int > > cell_level_index
Definition: matrix_free.h:908
void print_memory_consumption(StreamType &out) const
const IndexSet & get_locally_owned_set(const unsigned int fe_component=0) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:358
void compress_finish(::VectorOperation::values operation)
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
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:157
bool mapping_is_initialized
Definition: matrix_free.h:929
unsigned int n_macro_cells() const
internal::MatrixFreeFunctions::SizeInfo size_info
Definition: matrix_free.h:914
void initialize_dof_vector(VectorType &vec, const unsigned int vector_component=0) const
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices)
void clear()
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:249
#define Assert(cond, exc)
Definition: exceptions.h:294
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int fe_component=0) const
UpdateFlags
unsigned int n_components() const
bool indices_initialized() 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
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:191
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:201
const Quadrature< dim-1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
Definition: hp.h:102
void reinit(const size_type size, const bool omit_zeroing_entries=false)
const Number * constraint_pool_end(const unsigned int pool_index) const
std::size_t memory_consumption() const
unsigned int tasks_block_size
Definition: matrix_free.h:212
unsigned int get_dofs_per_cell(const unsigned int fe_component=0, const unsigned int hp_active_fe_index=0) 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 cell_loop(const std_cxx11::function< void(const MatrixFree< dim, Number > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src) const
const std_cxx11::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int vector_component=0) const
unsigned int level_mg_handler
Definition: matrix_free.h:235
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
Definition: matrix_free.h:225
void compress_start(const unsigned int communication_channel=0,::VectorOperation::values operation=VectorOperation::add)
Definition: table.h:33
void update_ghost_values_start(const unsigned int communication_channel=0) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
std::vector< internal::MatrixFreeFunctions::DoFInfo > dof_info
Definition: matrix_free.h:875
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
void update_ghost_values_finish() const
unsigned int n_components_filled(const unsigned int macro_cell_number) const