Reference documentation for deal.II version 8.1.0
matrix_free.h
1 // ---------------------------------------------------------------------
2 // @f$Id: matrix_free.h 31370 2013-10-21 08:49:41Z kronbichler @f$
3 //
4 // Copyright (C) 2011 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #ifndef __deal2__matrix_free_h
19 #define __deal2__matrix_free_h
20 
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>
41 
42 #ifdef DEAL_II_WITH_THREADS
43 #include <tbb/task.h>
44 #include <tbb/task_scheduler_init.h>
45 #include <tbb/parallel_for.h>
46 #include <tbb/blocked_range.h>
47 #endif
48 
49 #include <stdlib.h>
50 #include <memory>
51 #include <limits>
52 
53 
55 
56 
57 
106 template <int dim, typename Number=double>
108 {
109 public:
110 
140  {
144  enum TasksParallelScheme {none, partition_partition, partition_color, color};
145 
149  AdditionalData (const MPI_Comm mpi_communicator = MPI_COMM_SELF,
150  const TasksParallelScheme tasks_parallel_scheme = partition_partition,
151  const unsigned int tasks_block_size = 0,
153  const unsigned int level_mg_handler = numbers::invalid_unsigned_int,
154  const bool store_plain_indices = true,
155  const bool initialize_indices = true,
156  const bool initialize_mapping = true)
157  :
166  {};
167 
174  MPI_Comm mpi_communicator;
175 
206 
217  unsigned int tasks_block_size;
218 
231 
240  unsigned int level_mg_handler;
241 
249 
259 
267  };
268 
276  MatrixFree ();
277 
281  ~MatrixFree();
282 
299  template <typename DH, typename Quadrature>
300  void reinit (const Mapping<dim> &mapping,
301  const DH &dof_handler,
302  const ConstraintMatrix &constraint,
303  const IndexSet &locally_owned_dofs,
304  const Quadrature &quad,
305  const AdditionalData additional_data = AdditionalData());
306 
311  template <typename DH, typename Quadrature>
312  void reinit (const Mapping<dim> &mapping,
313  const DH &dof_handler,
314  const ConstraintMatrix &constraint,
315  const Quadrature &quad,
316  const AdditionalData additional_data = AdditionalData());
317 
322  template <typename DH, typename Quadrature>
323  void reinit (const DH &dof_handler,
324  const ConstraintMatrix &constraint,
325  const Quadrature &quad,
326  const AdditionalData additional_data = AdditionalData());
327 
354  template <typename DH, typename Quadrature>
355  void reinit (const Mapping<dim> &mapping,
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,
360  const AdditionalData additional_data = AdditionalData());
361 
367  template <typename DH, typename Quadrature>
368  void reinit (const Mapping<dim> &mapping,
369  const std::vector<const DH *> &dof_handler,
370  const std::vector<const ConstraintMatrix *> &constraint,
371  const std::vector<Quadrature> &quad,
372  const AdditionalData additional_data = AdditionalData());
373 
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,
382  const AdditionalData additional_data = AdditionalData());
383 
391  template <typename DH, typename Quadrature>
392  void reinit (const Mapping<dim> &mapping,
393  const std::vector<const DH *> &dof_handler,
394  const std::vector<const ConstraintMatrix *> &constraint,
395  const Quadrature &quad,
396  const AdditionalData additional_data = AdditionalData());
397 
402  template <typename DH, typename Quadrature>
403  void reinit (const std::vector<const DH *> &dof_handler,
404  const std::vector<const ConstraintMatrix *> &constraint,
405  const Quadrature &quad,
406  const AdditionalData additional_data = AdditionalData());
407 
413  void copy_from (const MatrixFree<dim,Number> &matrix_free_base);
414 
419  void clear();
420 
422 
440  template <typename OutVector, typename InVector>
441  void cell_loop (const std_cxx1x::function<void (const MatrixFree<dim,Number> &,
442  OutVector &,
443  const InVector &,
444  const std::pair<unsigned int,
445  unsigned int> &)> &cell_operation,
446  OutVector &dst,
447  const InVector &src) const;
448 
458  template <typename CLASS, typename OutVector, typename InVector>
459  void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
460  OutVector &,
461  const InVector &,
462  const std::pair<unsigned int,
463  unsigned int> &)const,
464  const CLASS *owning_class,
465  OutVector &dst,
466  const InVector &src) const;
467 
471  template <typename CLASS, typename OutVector, typename InVector>
472  void cell_loop (void (CLASS::*function_pointer)(const MatrixFree &,
473  OutVector &,
474  const InVector &,
475  const std::pair<unsigned int,
476  unsigned int> &),
477  CLASS *owning_class,
478  OutVector &dst,
479  const InVector &src) const;
480 
488  std::pair<unsigned int,unsigned int>
489  create_cell_subrange_hp (const std::pair<unsigned int,unsigned int> &range,
490  const unsigned int fe_degree,
491  const unsigned int vector_component = 0) const;
492 
499  std::pair<unsigned int,unsigned int>
500  create_cell_subrange_hp_by_index (const std::pair<unsigned int,unsigned int> &range,
501  const unsigned int fe_index,
502  const unsigned int vector_component = 0) const;
503 
505 
518  template <typename VectorType>
519  void initialize_dof_vector(VectorType &vec,
520  const unsigned int vector_component=0) const;
521 
530  template <typename Number2>
532  const unsigned int vector_component=0) const;
533 
544  const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &
545  get_vector_partitioner (const unsigned int vector_component=0) const;
546 
550  const IndexSet &
551  get_locally_owned_set (const unsigned int fe_component = 0) const;
552 
556  const IndexSet &
557  get_ghost_set (const unsigned int fe_component = 0) const;
558 
568  const std::vector<unsigned int> &
569  get_constrained_dofs (const unsigned int fe_component = 0) const;
570 
575  void renumber_dofs (std::vector<types::global_dof_index> &renumbering,
576  const unsigned int vector_component = 0);
577 
579 
587  unsigned int n_components () const;
588 
597  unsigned int n_physical_cells () const;
598 
609  unsigned int n_macro_cells () const;
610 
615  const DoFHandler<dim> &
616  get_dof_handler (const unsigned int fe_component = 0) const;
617 
629  get_cell_iterator (const unsigned int macro_cell_number,
630  const unsigned int vector_number,
631  const unsigned int fe_component = 0) const;
632 
644  typename hp::DoFHandler<dim>::active_cell_iterator
645  get_hp_cell_iterator (const unsigned int macro_cell_number,
646  const unsigned int vector_number,
647  const unsigned int fe_component = 0) const;
648 
660  bool
661  at_irregular_cell (const unsigned int macro_cell_number) const;
662 
671  unsigned int
672  n_components_filled (const unsigned int macro_cell_number) const;
673 
677  unsigned int
678  get_dofs_per_cell (const unsigned int fe_component = 0,
679  const unsigned int hp_active_fe_index = 0) const;
680 
684  unsigned int
685  get_n_q_points (const unsigned int quad_index = 0,
686  const unsigned int hp_active_fe_index = 0) const;
687 
692  unsigned int
693  get_dofs_per_face (const unsigned int fe_component = 0,
694  const unsigned int hp_active_fe_index = 0) const;
695 
700  unsigned int
701  get_n_q_points_face (const unsigned int quad_index = 0,
702  const unsigned int hp_active_fe_index = 0) const;
703 
707  const Quadrature<dim> &
708  get_quadrature (const unsigned int quad_index = 0,
709  const unsigned int hp_active_fe_index = 0) const;
710 
714  const Quadrature<dim-1> &
715  get_face_quadrature (const unsigned int quad_index = 0,
716  const unsigned int hp_active_fe_index = 0) const;
717 
721  bool indices_initialized () const;
722 
728  bool mapping_initialized () const;
729 
734  std::size_t memory_consumption() const;
735 
740  template <typename STREAM>
741  void print_memory_consumption(STREAM &out) const;
742 
747  void print (std::ostream &out) const;
748 
750 
759  get_task_info () const;
760 
765  get_size_info () const;
766 
767  /*
768  * Returns geometry-dependent information on the cells.
769  */
771  get_mapping_info () const;
772 
776  const internal::MatrixFreeFunctions::DoFInfo &
777  get_dof_info (const unsigned int fe_component = 0) const;
778 
782  unsigned int n_constraint_pool_entries() const;
783 
788  const Number *
789  constraint_pool_begin (const unsigned int pool_index) const;
790 
796  const Number *
797  constraint_pool_end (const unsigned int pool_index) const;
798 
803  get_shape_info (const unsigned int fe_component = 0,
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;
807 
809 
810 private:
811 
816  void internal_reinit (const Mapping<dim> &mapping,
817  const std::vector<const DoFHandler<dim> *> &dof_handler,
818  const std::vector<const ConstraintMatrix *> &constraint,
819  const std::vector<IndexSet> &locally_owned_set,
820  const std::vector<hp::QCollection<1> > &quad,
821  const AdditionalData additional_data);
822 
826  void internal_reinit (const Mapping<dim> &mapping,
827  const std::vector<const hp::DoFHandler<dim>*> &dof_handler,
828  const std::vector<const ConstraintMatrix *> &constraint,
829  const std::vector<IndexSet> &locally_owned_set,
830  const std::vector<hp::QCollection<1> > &quad,
831  const AdditionalData additional_data);
832 
839  void
840  initialize_indices (const std::vector<const ConstraintMatrix *> &constraint,
841  const std::vector<IndexSet> &locally_owned_set);
842 
846  void initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handlers,
847  const unsigned int level);
848 
852  void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
853  const unsigned int level);
854 
860  struct DoFHandlers
861  {
862  DoFHandlers () : n_dof_handlers (0), level (numbers::invalid_unsigned_int) {};
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;
867  unsigned int level;
868  };
869 
874 
879  std::vector<internal::MatrixFreeFunctions::DoFInfo> dof_info;
880 
887  std::vector<Number> constraint_pool_data;
888 
893  std::vector<unsigned int> constraint_pool_row_index;
894 
900 
905 
912  std::vector<std::pair<unsigned int,unsigned int> > cell_level_index;
913 
919 
924 
929 
934 };
935 
936 
937 
938 /*----------------------- Inline functions ----------------------------------*/
939 
940 #ifndef DOXYGEN
941 
942 
943 template <int dim, typename Number>
944 template <typename VectorType>
945 inline
946 void
948  const unsigned int comp) const
949 {
950  AssertIndexRange (comp, n_components());
951  vec.reinit(dof_info[comp].vector_partitioner->size());
952 }
953 
954 
955 
956 template <int dim, typename Number>
957 template <typename Number2>
958 inline
959 void
961  const unsigned int comp) const
962 {
963  AssertIndexRange (comp, n_components());
964  vec.reinit(dof_info[comp].vector_partitioner);
965 }
966 
967 
968 
969 template <int dim, typename Number>
970 inline
971 const std_cxx1x::shared_ptr<const Utilities::MPI::Partitioner> &
972 MatrixFree<dim,Number>::get_vector_partitioner (const unsigned int comp) const
973 {
974  AssertIndexRange (comp, n_components());
975  return dof_info[comp].vector_partitioner;
976 }
977 
978 
979 
980 template <int dim, typename Number>
981 inline
982 const std::vector<unsigned int> &
983 MatrixFree<dim,Number>::get_constrained_dofs (const unsigned int comp) const
984 {
985  AssertIndexRange (comp, n_components());
986  return dof_info[comp].constrained_dofs;
987 }
988 
989 
990 
991 template <int dim, typename Number>
992 inline
993 unsigned int
995 {
996  AssertDimension (dof_handlers.n_dof_handlers, dof_info.size());
997  return dof_handlers.n_dof_handlers;
998 }
999 
1000 
1001 
1002 template <int dim, typename Number>
1003 inline
1006 {
1007  return task_info;
1008 }
1009 
1010 
1011 
1012 template <int dim, typename Number>
1013 inline
1016 {
1017  return size_info;
1018 }
1019 
1020 
1021 
1022 template <int dim, typename Number>
1023 inline
1024 unsigned int
1026 {
1027  return size_info.n_macro_cells;
1028 }
1029 
1030 
1031 
1032 template <int dim, typename Number>
1033 inline
1034 unsigned int
1036 {
1037  return size_info.n_active_cells;
1038 }
1039 
1040 
1041 
1042 template <int dim, typename Number>
1043 inline
1046 {
1047  return mapping_info;
1048 }
1049 
1050 
1051 
1052 template <int dim, typename Number>
1053 inline
1054 const internal::MatrixFreeFunctions::DoFInfo &
1055 MatrixFree<dim,Number>::get_dof_info (unsigned int dof_index) const
1056 {
1057  AssertIndexRange (dof_index, n_components());
1058  return dof_info[dof_index];
1059 }
1060 
1061 
1062 
1063 template <int dim, typename Number>
1064 inline
1065 unsigned int
1067 {
1068  return constraint_pool_row_index.size()-1;
1069 }
1070 
1071 
1072 
1073 template <int dim, typename Number>
1074 inline
1075 const Number *
1076 MatrixFree<dim,Number>::constraint_pool_begin (const unsigned int row) const
1077 {
1078  AssertIndexRange (row, constraint_pool_row_index.size()-1);
1079  return constraint_pool_data.empty() ? 0 :
1080  &constraint_pool_data[0] + constraint_pool_row_index[row];
1081 }
1082 
1083 
1084 
1085 template <int dim, typename Number>
1086 inline
1087 const Number *
1088 MatrixFree<dim,Number>::constraint_pool_end (const unsigned int row) const
1089 {
1090  AssertIndexRange (row, constraint_pool_row_index.size()-1);
1091  return constraint_pool_data.empty() ? 0 :
1092  &constraint_pool_data[0] + constraint_pool_row_index[row+1];
1093 }
1094 
1095 
1096 
1097 template <int dim, typename Number>
1098 inline
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
1104 {
1105  AssertIndexRange (vector_component, dof_info.size());
1106  if (dof_info[vector_component].cell_active_fe_index.empty())
1107  {
1108  AssertDimension (dof_info[vector_component].fe_index_conversion.size(),1);
1109  if (dof_info[vector_component].fe_index_conversion[0].first == degree)
1110  return range;
1111  else
1112  return std::pair<unsigned int,unsigned int> (range.second,range.second);
1113  }
1114 
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);
1119  else
1120  return create_cell_subrange_hp_by_index (range, fe_index, vector_component);
1121 }
1122 
1123 
1124 
1125 template <int dim, typename Number>
1126 inline
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
1132 {
1133  AssertIndexRange (fe_index, dof_info[vector_component].max_fe_index);
1134  const std::vector<unsigned int> &fe_indices =
1135  dof_info[vector_component].cell_active_fe_index;
1136  if (fe_indices.size() == 0)
1137  return range;
1138  else
1139  {
1140  // the range over which we are searching must be ordered, otherwise we
1141  // got a range that spans over too many cells
1142 #ifdef DEBUG
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!"));
1146  AssertIndexRange(range.first,fe_indices.size()+1);
1147  AssertIndexRange(range.second,fe_indices.size()+1);
1148 #endif
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)
1153  -&fe_indices[0] ;
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 &&
1159  return_range.second <= range.second, ExcInternalError());
1160  return return_range;
1161  }
1162 }
1163 
1164 
1165 
1166 template <int dim, typename Number>
1167 inline
1168 void
1169 MatrixFree<dim,Number>::renumber_dofs (std::vector<types::global_dof_index> &renumbering,
1170  const unsigned int vector_component)
1171 {
1172  AssertIndexRange(vector_component, dof_info.size());
1173  dof_info[vector_component].renumber_dofs (renumbering);
1174 }
1175 
1176 
1177 
1178 template <int dim, typename Number>
1179 inline
1180 const DoFHandler<dim> &
1181 MatrixFree<dim,Number>::get_dof_handler (const unsigned int dof_index) const
1182 {
1183  AssertIndexRange (dof_index, n_components());
1184  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
1185  {
1186  AssertDimension (dof_handlers.dof_handler.size(),
1187  dof_handlers.n_dof_handlers);
1188  return *dof_handlers.dof_handler[dof_index];
1189  }
1190  else
1191  {
1192  Assert (false, ExcNotImplemented());
1193  // put pseudo return argument to avoid compiler error, but trigger a
1194  // segfault in case this is only run in optimized mode
1195  return *dof_handlers.dof_handler[numbers::invalid_unsigned_int];
1196  }
1197 }
1198 
1199 
1200 
1201 template <int dim, typename Number>
1202 inline
1204 MatrixFree<dim,Number>::get_cell_iterator(const unsigned int macro_cell_number,
1205  const unsigned int vector_number,
1206  const unsigned int dof_index) const
1207 {
1208  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
1209 #ifdef DEBUG
1210  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
1211  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
1212  AssertIndexRange (vector_number, vectorization_length);
1213  const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
1214  if (irreg_filled > 0)
1215  AssertIndexRange (vector_number, irreg_filled);
1216 #endif
1217 
1218  const DoFHandler<dim> *dofh = 0;
1219  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
1220  {
1221  AssertDimension (dof_handlers.dof_handler.size(),
1222  dof_handlers.n_dof_handlers);
1223  dofh = dof_handlers.dof_handler[dof_index];
1224  }
1225  else
1226  {
1227  Assert (false, ExcMessage ("Cannot return DoFHandler<dim>::cell_iterator "
1228  "for underlying DoFHandler!"));
1229  }
1230 
1231  std::pair<unsigned int,unsigned int> index =
1232  cell_level_index[macro_cell_number*vectorization_length+vector_number];
1233  return typename DoFHandler<dim>::cell_iterator
1234  (&dofh->get_tria(), index.first, index.second, dofh);
1235 }
1236 
1237 
1238 
1239 template <int dim, typename Number>
1240 inline
1241 typename hp::DoFHandler<dim>::active_cell_iterator
1242 MatrixFree<dim,Number>::get_hp_cell_iterator(const unsigned int macro_cell_number,
1243  const unsigned int vector_number,
1244  const unsigned int dof_index) const
1245 {
1246  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
1247 #ifdef DEBUG
1248  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
1249  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
1250  AssertIndexRange (vector_number, vectorization_length);
1251  const unsigned int irreg_filled = dof_info[dof_index].row_starts[macro_cell_number][2];
1252  if (irreg_filled > 0)
1253  AssertIndexRange (vector_number, irreg_filled);
1254 #endif
1255 
1256  Assert (dof_handlers.active_dof_handler == DoFHandlers::hp,
1257  ExcNotImplemented());
1258  const hp::DoFHandler<dim> *dofh = dof_handlers.hp_dof_handler[dof_index];
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);
1263 }
1264 
1265 
1266 
1267 template <int dim, typename Number>
1268 inline
1269 bool
1270 MatrixFree<dim,Number>::at_irregular_cell (const unsigned int macro_cell) const
1271 {
1272  AssertIndexRange (macro_cell, size_info.n_macro_cells);
1273  return dof_info[0].row_starts[macro_cell][2] > 0;
1274 }
1275 
1276 
1277 
1278 template <int dim, typename Number>
1279 inline
1280 unsigned int
1281 MatrixFree<dim,Number>::n_components_filled (const unsigned int macro_cell) const
1282 {
1283  AssertIndexRange (macro_cell, size_info.n_macro_cells);
1284  const unsigned int n_filled = dof_info[0].row_starts[macro_cell][2];
1285  if (n_filled == 0)
1287  else
1288  return n_filled;
1289 }
1290 
1291 
1292 
1293 template <int dim, typename Number>
1294 inline
1295 unsigned int
1296 MatrixFree<dim,Number>::get_dofs_per_cell(const unsigned int dof_index,
1297  const unsigned int active_fe_index) const
1298 {
1299  AssertIndexRange (dof_index, dof_info.size());
1300  return dof_info[dof_index].dofs_per_cell[active_fe_index];
1301 }
1302 
1303 
1304 
1305 template <int dim, typename Number>
1306 inline
1307 unsigned int
1308 MatrixFree<dim,Number>::get_n_q_points(const unsigned int quad_index,
1309  const unsigned int active_fe_index) const
1310 {
1311  AssertIndexRange (quad_index,
1312  mapping_info.mapping_data_gen.size());
1313  return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index];
1314 }
1315 
1316 
1317 
1318 template <int dim, typename Number>
1319 inline
1320 unsigned int
1321 MatrixFree<dim,Number>::get_dofs_per_face(const unsigned int dof_index,
1322  const unsigned int active_fe_index) const
1323 {
1324  AssertIndexRange (dof_index, dof_info.size());
1325  return dof_info[dof_index].dofs_per_face[active_fe_index];
1326 }
1327 
1328 
1329 
1330 template <int dim, typename Number>
1331 inline
1332 unsigned int
1333 MatrixFree<dim,Number>::get_n_q_points_face(const unsigned int quad_index,
1334  const unsigned int active_fe_index) const
1335 {
1336  AssertIndexRange (quad_index,
1337  mapping_info.mapping_data_gen.size());
1338  return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index];
1339 }
1340 
1341 
1342 
1343 template <int dim, typename Number>
1344 inline
1345 const IndexSet &
1346 MatrixFree<dim,Number>::get_locally_owned_set(const unsigned int dof_index) const
1347 {
1348  AssertIndexRange (dof_index, dof_info.size());
1349  return dof_info[dof_index].vector_partitioner->locally_owned_range();
1350 }
1351 
1352 
1353 
1354 template <int dim, typename Number>
1355 inline
1356 const IndexSet &
1357 MatrixFree<dim,Number>::get_ghost_set(const unsigned int dof_index) const
1358 {
1359  AssertIndexRange (dof_index, dof_info.size());
1360  return dof_info[dof_index].vector_partitioner->ghost_indices();
1361 }
1362 
1363 
1364 
1365 template <int dim, typename Number>
1366 inline
1368 MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
1369  const unsigned int index_quad,
1370  const unsigned int active_fe_index,
1371  const unsigned int active_quad_index) const
1372 {
1373  AssertIndexRange (index_fe, shape_info.size(0));
1374  AssertIndexRange (index_quad, shape_info.size(1));
1375  AssertIndexRange (active_fe_index, shape_info.size(2));
1376  AssertIndexRange (active_quad_index, shape_info.size(3));
1377  return shape_info(index_fe, index_quad,
1378  active_fe_index, active_quad_index);
1379 }
1380 
1381 
1382 
1383 template <int dim, typename Number>
1384 inline
1385 const Quadrature<dim> &
1386 MatrixFree<dim,Number>::get_quadrature (const unsigned int quad_index,
1387  const unsigned int active_fe_index) const
1388 {
1389  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
1390  return mapping_info.mapping_data_gen[quad_index].
1391  quadrature[active_fe_index];
1392 }
1393 
1394 
1395 
1396 template <int dim, typename Number>
1397 inline
1398 const Quadrature<dim-1> &
1399 MatrixFree<dim,Number>::get_face_quadrature (const unsigned int quad_index,
1400  const unsigned int active_fe_index) const
1401 {
1402  AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
1403  return mapping_info.mapping_data_gen[quad_index].
1404  face_quadrature[active_fe_index];
1405 }
1406 
1407 
1408 
1409 template <int dim, typename Number>
1410 inline
1411 bool
1413 {
1414  return indices_are_initialized;
1415 }
1416 
1417 
1418 
1419 template <int dim, typename Number>
1420 inline
1421 bool
1423 {
1424  return mapping_is_initialized;
1425 }
1426 
1427 
1428 
1429 // ------------------------------ reinit functions ---------------------------
1430 
1431 namespace internal
1432 {
1433  namespace MatrixFree
1434  {
1435  template <typename DH>
1436  inline
1437  std::vector<IndexSet>
1438  extract_locally_owned_index_sets (const std::vector<const DH *> &dofh,
1439  const unsigned int level)
1440  {
1441  std::vector<IndexSet> locally_owned_set;
1442  locally_owned_set.reserve (dofh.size());
1443  for (unsigned int j=0; j<dofh.size(); j++)
1444  if (level == numbers::invalid_unsigned_int)
1445  locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
1446  else
1447  {
1448  // TODO: not distributed yet
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);
1452  }
1453  return locally_owned_set;
1454  }
1455  }
1456 }
1457 
1458 
1459 
1460 template <int dim, typename Number>
1461 template <typename DH, typename Quad>
1463 reinit(const DH &dof_handler,
1464  const ConstraintMatrix &constraints_in,
1465  const Quad &quad,
1466  const MatrixFree<dim,Number>::AdditionalData additional_data)
1467 {
1468  MappingQ1<dim> mapping;
1469  std::vector<const DH *> dof_handlers;
1470  std::vector<const ConstraintMatrix *> constraints;
1471  std::vector<Quad> quads;
1472 
1473  dof_handlers.push_back(&dof_handler);
1474  constraints.push_back (&constraints_in);
1475  quads.push_back (quad);
1476 
1477  std::vector<IndexSet> locally_owned_sets =
1478  internal::MatrixFree::extract_locally_owned_index_sets
1479  (dof_handlers, additional_data.level_mg_handler);
1480  reinit(mapping, dof_handlers,constraints, locally_owned_sets, quads,
1481  additional_data);
1482 }
1483 
1484 
1485 
1486 template <int dim, typename Number>
1487 template <typename DH, typename Quad>
1489 reinit(const Mapping<dim> &mapping,
1490  const DH &dof_handler,
1491  const ConstraintMatrix &constraints_in,
1492  const Quad &quad,
1493  const MatrixFree<dim,Number>::AdditionalData additional_data)
1494 {
1495  std::vector<const DH *> dof_handlers;
1496  std::vector<const ConstraintMatrix *> constraints;
1497  std::vector<Quad> quads;
1498 
1499  dof_handlers.push_back(&dof_handler);
1500  constraints.push_back (&constraints_in);
1501  quads.push_back (quad);
1502 
1503  std::vector<IndexSet> locally_owned_sets =
1504  internal::MatrixFree::extract_locally_owned_index_sets
1505  (dof_handlers, additional_data.level_mg_handler);
1506  reinit(mapping, dof_handlers,constraints,locally_owned_sets, quads,
1507  additional_data);
1508 }
1509 
1510 
1511 
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,
1518  const MatrixFree<dim,Number>::AdditionalData additional_data)
1519 {
1520  MappingQ1<dim> mapping;
1521  std::vector<IndexSet> locally_owned_set =
1522  internal::MatrixFree::extract_locally_owned_index_sets
1523  (dof_handler, additional_data.level_mg_handler);
1524  reinit(mapping, dof_handler,constraint,locally_owned_set,
1525  static_cast<const std::vector<Quadrature<1> >&>(quad),
1526  additional_data);
1527 }
1528 
1529 
1530 
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,
1536  const Quad &quad,
1537  const MatrixFree<dim,Number>::AdditionalData additional_data)
1538 {
1539  MappingQ1<dim> mapping;
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
1544  (dof_handler, additional_data.level_mg_handler);
1545  reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
1546  additional_data);
1547 }
1548 
1549 
1550 
1551 template <int dim, typename Number>
1552 template <typename DH, typename Quad>
1554 reinit(const Mapping<dim> &mapping,
1555  const std::vector<const DH *> &dof_handler,
1556  const std::vector<const ConstraintMatrix *> &constraint,
1557  const Quad &quad,
1558  const MatrixFree<dim,Number>::AdditionalData additional_data)
1559 {
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
1564  (dof_handler, additional_data.level_mg_handler);
1565  reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
1566  additional_data);
1567 }
1568 
1569 
1570 
1571 template <int dim, typename Number>
1572 template <typename DH, typename Quad>
1574 reinit(const Mapping<dim> &mapping,
1575  const std::vector<const DH *> &dof_handler,
1576  const std::vector<const ConstraintMatrix *> &constraint,
1577  const std::vector<Quad> &quad,
1578  const MatrixFree<dim,Number>::AdditionalData additional_data)
1579 {
1580  std::vector<IndexSet> locally_owned_set =
1581  internal::MatrixFree::extract_locally_owned_index_sets
1582  (dof_handler, additional_data.level_mg_handler);
1583  reinit(mapping, dof_handler,constraint,locally_owned_set,
1584  quad, additional_data);
1585 }
1586 
1587 
1588 
1589 namespace internal
1590 {
1591  namespace MatrixFree
1592  {
1593  // resolve DoFHandler types
1594 
1595  // MGDoFHandler is deprecated in deal.II but might still be present in
1596  // user code, so we need to resolve its type (fortunately, it is derived
1597  // from DoFHandler, so we can static_cast it to a DoFHandler<dim>)
1598  template <typename DH>
1599  inline
1600  std::vector<const ::DoFHandler<DH::dimension> *>
1601  resolve_dof_handler (const std::vector<const DH *> &dof_handler)
1602  {
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]);
1606  return conversion;
1607  }
1608 
1609  template <int dim>
1610  inline
1611  std::vector<const ::hp::DoFHandler<dim> *>
1612  resolve_dof_handler (const std::vector<const ::hp::DoFHandler<dim> *> &dof_handler)
1613  {
1614  return dof_handler;
1615  }
1616  }
1617 }
1618 
1619 
1620 
1621 template <int dim, typename Number>
1622 template <typename DH, typename Quad>
1624 reinit(const Mapping<dim> &mapping,
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,
1629  const MatrixFree<dim,Number>::AdditionalData additional_data)
1630 {
1631  // find out whether we use a hp Quadrature or a standard quadrature
1632  std::vector<hp::QCollection<1> > quad_hp;
1633  for (unsigned int q=0; q<quad.size(); ++q)
1634  quad_hp.push_back (hp::QCollection<1>(quad[q]));
1635  internal_reinit (mapping,
1636  internal::MatrixFree::resolve_dof_handler(dof_handler),
1637  constraint, locally_owned_set, quad_hp, additional_data);
1638 }
1639 
1640 
1641 
1642 // ------------------------------ implementation of cell_loop ---------------
1643 
1644 // internal helper functions that define how to call MPI data exchange
1645 // functions: for generic vectors, do nothing at all. For distributed vectors,
1646 // call update_ghost_values_start function and so on. If we have collections
1647 // of vectors, just do the individual functions of the components. In order to
1648 // keep ghost values consistent (whether we are in read or write mode). the whole situation is a bit complicated by the fact
1649 // that we need to treat block vectors differently, which use some additional
1650 // helper functions to select the blocks and template magic.
1651 namespace internal
1652 {
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,
1671 
1672  template<typename VectorStruct>
1673  bool update_ghost_values_start_block (const VectorStruct &,
1674  const unsigned int,
1676  {
1677  return false;
1678  }
1679  template<typename VectorStruct>
1680  void reset_ghost_values_block (const VectorStruct &,
1681  const bool,
1683  {}
1684  template<typename VectorStruct>
1685  void update_ghost_values_finish_block (const VectorStruct &,
1687  {}
1688  template<typename VectorStruct>
1689  void compress_start_block (const VectorStruct &,
1690  const unsigned int,
1692  {}
1693  template<typename VectorStruct>
1694  void compress_finish_block (const VectorStruct &,
1696  {}
1697 
1698 
1699 
1700  // returns true if the vector was in a state without ghost values before,
1701  // i.e., we need to zero out ghosts in the very end
1702  template<typename VectorStruct>
1703  inline
1704  bool update_ghost_values_start (const VectorStruct &vec,
1705  const unsigned int channel = 0)
1706  {
1707  return
1708  update_ghost_values_start_block(vec, channel,
1710  }
1711 
1712 
1713 
1714  template<typename Number>
1715  inline
1716  bool update_ghost_values_start (const parallel::distributed::Vector<Number> &vec,
1717  const unsigned int channel = 0)
1718  {
1719  bool return_value = !vec.has_ghost_elements();
1720  vec.update_ghost_values_start(channel);
1721  return return_value;
1722  }
1723 
1724 
1725 
1726  template <typename VectorStruct>
1727  inline
1728  bool update_ghost_values_start (const std::vector<VectorStruct> &vec)
1729  {
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;
1734  }
1735 
1736 
1737 
1738  template <typename VectorStruct>
1739  inline
1740  bool update_ghost_values_start (const std::vector<VectorStruct *> &vec)
1741  {
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;
1746  }
1747 
1748 
1749 
1750  template<typename VectorStruct>
1751  inline
1752  bool update_ghost_values_start_block (const VectorStruct &vec,
1753  const unsigned int channel,
1755  {
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;
1760  }
1761 
1762 
1763 
1764  // if the input vector did not have ghosts imported, clear them here again
1765  // in order to avoid subsequent operations e.g. in linear solvers to work
1766  // with ghosts all the time
1767  template<typename VectorStruct>
1768  inline
1769  void reset_ghost_values (const VectorStruct &vec,
1770  const bool zero_out_ghosts)
1771  {
1772  reset_ghost_values_block(vec, zero_out_ghosts,
1774  }
1775 
1776 
1777 
1778  template<typename Number>
1779  inline
1780  void reset_ghost_values (const parallel::distributed::Vector<Number> &vec,
1781  const bool zero_out_ghosts)
1782  {
1783  if (zero_out_ghosts)
1784  const_cast<parallel::distributed::Vector<Number>&>(vec).zero_out_ghosts();
1785  }
1786 
1787 
1788 
1789  template <typename VectorStruct>
1790  inline
1791  void reset_ghost_values (const std::vector<VectorStruct> &vec,
1792  const bool zero_out_ghosts)
1793  {
1794  for (unsigned int comp=0; comp<vec.size(); comp++)
1795  reset_ghost_values(vec[comp], zero_out_ghosts);
1796  }
1797 
1798 
1799 
1800  template <typename VectorStruct>
1801  inline
1802  void reset_ghost_values (const std::vector<VectorStruct *> &vec,
1803  const bool zero_out_ghosts)
1804  {
1805  for (unsigned int comp=0; comp<vec.size(); comp++)
1806  reset_ghost_values(*vec[comp], zero_out_ghosts);
1807  }
1808 
1809 
1810 
1811  template<typename VectorStruct>
1812  inline
1813  void reset_ghost_values_block (const VectorStruct &vec,
1814  const bool zero_out_ghosts,
1816  {
1817  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1818  reset_ghost_values(vec.block(i), zero_out_ghosts);
1819  }
1820 
1821 
1822 
1823  template <typename VectorStruct>
1824  inline
1825  void update_ghost_values_finish (const VectorStruct &vec)
1826  {
1827  update_ghost_values_finish_block(vec,
1829  }
1830 
1831 
1832 
1833  template <typename Number>
1834  inline
1835  void update_ghost_values_finish (const parallel::distributed::Vector<Number> &vec)
1836  {
1838  }
1839 
1840 
1841 
1842  template <typename VectorStruct>
1843  inline
1844  void update_ghost_values_finish (const std::vector<VectorStruct> &vec)
1845  {
1846  for (unsigned int comp=0; comp<vec.size(); comp++)
1847  update_ghost_values_finish(vec[comp]);
1848  }
1849 
1850 
1851 
1852  template <typename VectorStruct>
1853  inline
1854  void update_ghost_values_finish (const std::vector<VectorStruct *> &vec)
1855  {
1856  for (unsigned int comp=0; comp<vec.size(); comp++)
1857  update_ghost_values_finish(*vec[comp]);
1858  }
1859 
1860 
1861 
1862  template <typename VectorStruct>
1863  inline
1864  void update_ghost_values_finish_block (const VectorStruct &vec,
1866  {
1867  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1868  update_ghost_values_finish(vec.block(i));
1869  }
1870 
1871 
1872 
1873  template <typename VectorStruct>
1874  inline
1875  void compress_start (VectorStruct &vec,
1876  const unsigned int channel = 0)
1877  {
1878  compress_start_block (vec, channel,
1880  }
1881 
1882 
1883 
1884  template <typename Number>
1885  inline
1886  void compress_start (parallel::distributed::Vector<Number> &vec,
1887  const unsigned int channel = 0)
1888  {
1889  vec.compress_start(channel);
1890  }
1891 
1892 
1893 
1894  template <typename VectorStruct>
1895  inline
1896  void compress_start (std::vector<VectorStruct> &vec)
1897  {
1898  for (unsigned int comp=0; comp<vec.size(); comp++)
1899  compress_start (vec[comp], comp);
1900  }
1901 
1902 
1903 
1904  template <typename VectorStruct>
1905  inline
1906  void compress_start (std::vector<VectorStruct *> &vec)
1907  {
1908  for (unsigned int comp=0; comp<vec.size(); comp++)
1909  compress_start (*vec[comp], comp);
1910  }
1911 
1912 
1913 
1914  template <typename VectorStruct>
1915  inline
1916  void compress_start_block (VectorStruct &vec,
1917  const unsigned int channel,
1919  {
1920  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1921  compress_start(vec.block(i), channel + 500*i);
1922  }
1923 
1924 
1925 
1926  template <typename VectorStruct>
1927  inline
1928  void compress_finish (VectorStruct &vec)
1929  {
1930  compress_finish_block(vec,
1932  }
1933 
1934 
1935 
1936  template <typename Number>
1937  inline
1938  void compress_finish (parallel::distributed::Vector<Number> &vec)
1939  {
1940  vec.compress_finish(::VectorOperation::add);
1941  }
1942 
1943 
1944 
1945  template <typename VectorStruct>
1946  inline
1947  void compress_finish (std::vector<VectorStruct> &vec)
1948  {
1949  for (unsigned int comp=0; comp<vec.size(); comp++)
1950  compress_finish(vec[comp]);
1951  }
1952 
1953 
1954 
1955  template <typename VectorStruct>
1956  inline
1957  void compress_finish (std::vector<VectorStruct *> &vec)
1958  {
1959  for (unsigned int comp=0; comp<vec.size(); comp++)
1960  compress_finish(*vec[comp]);
1961  }
1962 
1963 
1964 
1965  template <typename VectorStruct>
1966  inline
1967  void compress_finish_block (VectorStruct &vec,
1969  {
1970  for (unsigned int i=0; i<vec.n_blocks(); ++i)
1971  compress_finish(vec.block(i));
1972  }
1973 
1974 
1975 
1976 #ifdef DEAL_II_WITH_THREADS
1977 
1978  // This defines the TBB data structures that are needed to schedule the
1979  // partition-partition variant
1980 
1981  namespace partition
1982  {
1983  template<typename Worker, bool blocked=false>
1984  class CellWork : public tbb::task
1985  {
1986  public:
1987  CellWork (const Worker &worker_in,
1988  const unsigned int partition_in,
1989  const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
1990  :
1991  worker (worker_in),
1992  partition (partition_in),
1993  task_info (task_info_in)
1994  {};
1995  tbb::task *execute ()
1996  {
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]);
2000  worker(cell_range);
2001  if (blocked==true)
2002  dummy->spawn (*dummy);
2003  return NULL;
2004  }
2005 
2006  tbb::empty_task *dummy;
2007 
2008  private:
2009  const Worker &worker;
2010  const unsigned int partition;
2011  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2012  };
2013 
2014 
2015 
2016  template<typename Worker, bool blocked=false>
2017  class PartitionWork : public tbb::task
2018  {
2019  public:
2020  PartitionWork (const Worker &function_in,
2021  const unsigned int partition_in,
2022  const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
2023  :
2024  function (function_in),
2025  partition (partition_in),
2026  task_info (task_info_in)
2027  {};
2028  tbb::task *execute ()
2029  {
2030  if (false)
2031  {
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);
2038  }
2039  else
2040  {
2041  tbb::empty_task *root = new( tbb::task::allocate_root() )
2042  tbb::empty_task;
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);
2050 
2051  root->set_ref_count(evens+1);
2052  for (unsigned int j=0; j<evens; j++)
2053  {
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);
2058  if (j>0)
2059  {
2060  worker[j]->set_ref_count(2);
2061  blocked_worker[j-1]->dummy = new(worker[j]->allocate_child())
2062  tbb::empty_task;
2063  worker[j-1]->spawn(*blocked_worker[j-1]);
2064  }
2065  else
2066  worker[j]->set_ref_count(1);
2067  if (j<evens-1)
2068  {
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);
2073  }
2074  else
2075  {
2076  if (odds==evens)
2077  {
2078  worker[evens] = new(worker[j]->allocate_child())
2079  CellWork<Worker,false>(function,
2080  task_info.
2081  partition_color_blocks_row_index
2082  [partition]+2*j+1,task_info);
2083  worker[j]->spawn(*worker[evens]);
2084  }
2085  else
2086  {
2087  tbb::empty_task *child = new(worker[j]->allocate_child())
2088  tbb::empty_task();
2089  worker[j]->spawn(*child);
2090  }
2091  }
2092  }
2093 
2094  root->wait_for_all();
2095  root->destroy(*root);
2096  }
2097  if (blocked==true)
2098  dummy->spawn (*dummy);
2099  return NULL;
2100  }
2101 
2102  tbb::empty_task *dummy;
2103 
2104  private:
2105  const Worker &function;
2106  const unsigned int partition;
2107  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2108  };
2109 
2110  } // end of namespace partition
2111 
2112 
2113 
2114  namespace color
2115  {
2116  template <typename Worker>
2117  class CellWork
2118  {
2119  public:
2120  CellWork (const Worker &worker_in,
2121  const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
2122  :
2123  worker (worker_in),
2124  task_info (task_info_in)
2125  {};
2126  void operator()(const tbb::blocked_range<unsigned int> &r) const
2127  {
2128  for (unsigned int block=r.begin(); block<r.end(); block++)
2129  {
2130  std::pair<unsigned int,unsigned int> cell_range;
2131  if (task_info.position_short_block<block)
2132  {
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;
2136  }
2137  else
2138  {
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));
2143  }
2144  worker (cell_range);
2145  }
2146  }
2147  private:
2148  const Worker &worker;
2149  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2150  };
2151 
2152 
2153  template<typename Worker, bool blocked=false>
2154  class PartitionWork : public tbb::task
2155  {
2156  public:
2157  PartitionWork (const Worker &worker_in,
2158  const unsigned int partition_in,
2159  const internal::MatrixFreeFunctions::TaskInfo &task_info_in)
2160  :
2161  worker (worker_in),
2162  partition (partition_in),
2163  task_info (task_info_in)
2164  {};
2165  tbb::task *execute ()
2166  {
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));
2171  if (blocked==true)
2172  dummy->spawn (*dummy);
2173  return NULL;
2174  }
2175 
2176  tbb::empty_task *dummy;
2177 
2178  private:
2179  const Worker &worker;
2180  const unsigned int partition;
2181  const internal::MatrixFreeFunctions::TaskInfo &task_info;
2182  };
2183 
2184  } // end of namespace color
2185 
2186 
2187  template<typename VectorStruct>
2188  class MPIComDistribute : public tbb::task
2189  {
2190  public:
2191  MPIComDistribute (const VectorStruct &src_in)
2192  :
2193  src(src_in)
2194  {};
2195 
2196  tbb::task *execute ()
2197  {
2198  internal::update_ghost_values_finish(src);
2199  return 0;
2200  }
2201 
2202  private:
2203  const VectorStruct &src;
2204  };
2205 
2206 
2207 
2208  template<typename VectorStruct>
2209  class MPIComCompress : public tbb::task
2210  {
2211  public:
2212  MPIComCompress (VectorStruct &dst_in)
2213  :
2214  dst(dst_in)
2215  {};
2216 
2217  tbb::task *execute ()
2218  {
2219  internal::compress_start(dst);
2220  return 0;
2221  }
2222 
2223  private:
2224  VectorStruct &dst;
2225  };
2226 
2227 #endif // DEAL_II_WITH_THREADS
2228 
2229 } // end of namespace internal
2230 
2231 
2232 
2233 template <int dim, typename Number>
2234 template <typename OutVector, typename InVector>
2235 inline
2236 void
2238 (const std_cxx1x::function<void (const MatrixFree<dim,Number> &,
2239  OutVector &,
2240  const InVector &,
2241  const std::pair<unsigned int,
2242  unsigned int> &)> &cell_operation,
2243  OutVector &dst,
2244  const InVector &src) const
2245 {
2246  // in any case, need to start the ghost import at the beginning
2247  bool ghosts_were_not_set = internal::update_ghost_values_start (src);
2248 
2249 #ifdef DEAL_II_WITH_THREADS
2250 
2251  // Use multithreading if so requested and if there is enough work to do in
2252  // parallel (the code might hang if there are less than two chunks!)
2253  if (task_info.use_multithreading == true && task_info.n_blocks > 3)
2254  {
2255  // to simplify the function calls, bind away all arguments except the
2256  // cell range
2257  typedef
2258  std_cxx1x::function<void (const std::pair<unsigned int,unsigned int> &range)>
2259  Worker;
2260 
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),
2265  std_cxx1x::_1);
2266 
2267  if (task_info.use_partition_partition == true)
2268  {
2269  tbb::empty_task *root = new( tbb::task::allocate_root() )
2270  tbb::empty_task;
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>*>
2277  worker(n_workers);
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++)
2285  {
2286  if (j>0)
2287  {
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())
2293  tbb::empty_task;
2294  if (j>1)
2295  worker[j-1]->spawn(*blocked_worker[j-1]);
2296  else
2297  worker_compr->spawn(*blocked_worker[j-1]);
2298  }
2299  else
2300  {
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);
2309  }
2310  if (j<evens-1)
2311  {
2312  blocked_worker[j] = new(worker[j]->allocate_child())
2313  internal::partition::PartitionWork<Worker,true>
2314  (func,2*j+1,task_info);
2315  }
2316  else
2317  {
2318  if (odds==evens)
2319  {
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]);
2324  }
2325  else
2326  {
2327  tbb::empty_task *child = new(worker[j]->allocate_child())
2328  tbb::empty_task();
2329  worker[j]->spawn(*child);
2330  }
2331  }
2332  }
2333 
2334  root->wait_for_all();
2335  root->destroy(*root);
2336  }
2337  else // end of partition-partition, start of partition-color
2338  {
2339  unsigned int evens = task_info.evens;
2340  unsigned int odds = task_info.odds;
2341 
2342  // check whether there is only one partition. if not, build up the
2343  // tree of partitions
2344  if (odds > 0)
2345  {
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-
2350  n_blocked_workers;
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++)
2361  {
2362  spawn_index_new = worker_index;
2363  if (part == 0)
2364  worker[worker_index] = new(worker_compr->allocate_child())
2365  internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2366  else
2367  worker[worker_index] = new(root->allocate_child())
2368  internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2369  slice_index++;
2370  for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
2371  slice_index++)
2372  {
2373  worker[worker_index]->set_ref_count(1);
2374  worker_index++;
2375  worker[worker_index] = new (worker[worker_index-1]->allocate_child())
2376  internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2377  }
2378  worker[worker_index]->set_ref_count(2);
2379  if (part>0)
2380  {
2381  blocked_worker[(part-1)/2]->dummy =
2382  new (worker[worker_index]->allocate_child()) tbb::empty_task;
2383  worker_index++;
2384  if (spawn_index_child == -1)
2385  worker[spawn_index]->spawn(*blocked_worker[(part-1)/2]);
2386  else
2387  worker[spawn_index]->spawn(*worker[spawn_index_child]);
2388  spawn_index = spawn_index_new;
2389  spawn_index_child = -2;
2390  }
2391  else
2392  {
2393  internal::MPIComDistribute<InVector> *worker_dist =
2394  new (worker[worker_index]->allocate_child())
2395  internal::MPIComDistribute<InVector>(src);
2396  worker_dist->spawn(*worker_dist);
2397  worker_index++;
2398  }
2399  part += 1;
2400  if (part<task_info.partition_color_blocks_row_index.size()-1)
2401  {
2402  if (part<task_info.partition_color_blocks_row_index.size()-2)
2403  {
2404  blocked_worker[part/2] = new(worker[worker_index-1]->allocate_child())
2405  internal::color::PartitionWork<Worker,true>(func,slice_index,task_info);
2406  slice_index++;
2407  if (slice_index<
2408  task_info.partition_color_blocks_row_index[part+1])
2409  {
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);
2413  slice_index++;
2414  }
2415  else
2416  {
2417  spawn_index_child = -1;
2418  continue;
2419  }
2420  }
2421  for (; slice_index<task_info.partition_color_blocks_row_index[part+1];
2422  slice_index++)
2423  {
2424  if (slice_index>
2425  task_info.partition_color_blocks_row_index[part])
2426  {
2427  worker[worker_index]->set_ref_count(1);
2428  worker_index++;
2429  }
2430  worker[worker_index] = new (worker[worker_index-1]->allocate_child())
2431  internal::color::PartitionWork<Worker,false>(func,slice_index,task_info);
2432  }
2433  spawn_index_child = worker_index;
2434  worker_index++;
2435  }
2436  else
2437  {
2438  tbb::empty_task *final = new (worker[worker_index-1]->allocate_child())
2439  tbb::empty_task;
2440  worker[spawn_index]->spawn(*final);
2441  spawn_index_child = worker_index-1;
2442  }
2443  }
2444  if (evens==odds)
2445  worker[spawn_index]->spawn(*worker[spawn_index_child]);
2446  root->wait_for_all();
2447  root->destroy(*root);
2448  }
2449  // case when we only have one partition: this is the usual coloring
2450  // scheme, and we just schedule a parallel for loop for each color
2451  else
2452  {
2453  Assert(evens==1,ExcInternalError());
2454  internal::update_ghost_values_finish(src);
2455 
2456  for (unsigned int color=0;
2457  color < task_info.partition_color_blocks_row_index[1];
2458  ++color)
2459  {
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>
2464  (func,task_info));
2465  }
2466 
2467  internal::compress_start(dst);
2468  }
2469  }
2470  }
2471  else
2472 #endif
2473  // serial loop
2474  {
2475  std::pair<unsigned int,unsigned int> cell_range;
2476 
2477  // First operate on cells where no ghost data is needed (inner cells)
2478  {
2479  cell_range.first = 0;
2480  cell_range.second = size_info.boundary_cells_start;
2481  cell_operation (*this, dst, src, cell_range);
2482  }
2483 
2484  // before starting operations on cells that contain ghost nodes (outer
2485  // cells), wait for the MPI commands to finish
2486  internal::update_ghost_values_finish(src);
2487 
2488  // For the outer cells, do the same procedure as for inner cells.
2489  if (size_info.boundary_cells_end > size_info.boundary_cells_start)
2490  {
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);
2494  }
2495 
2496  internal::compress_start(dst);
2497 
2498  // Finally operate on cells where no ghost data is needed (inner cells)
2499  if (size_info.n_macro_cells > size_info.boundary_cells_end)
2500  {
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);
2504  }
2505  }
2506 
2507  // In every case, we need to finish transfers at the very end
2508  internal::compress_finish(dst);
2509  internal::reset_ghost_values(src, ghosts_were_not_set);
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> &)const,
2524  const CLASS *owning_class,
2525  OutVector &dst,
2526  const InVector &src) const
2527 {
2528  // here, use std_cxx1x::bind to hand a function handler with the appropriate
2529  // argument to the other loop function
2530  std_cxx1x::function<void (const MatrixFree<dim,Number> &,
2531  OutVector &,
2532  const InVector &,
2533  const std::pair<unsigned int,
2534  unsigned int> &)>
2535  function = std_cxx1x::bind<void>(function_pointer,
2536  std_cxx1x::cref(*owning_class),
2537  std_cxx1x::_1,
2538  std_cxx1x::_2,
2539  std_cxx1x::_3,
2540  std_cxx1x::_4);
2541  cell_loop (function, dst, src);
2542 }
2543 
2544 
2545 
2546 template <int dim, typename Number>
2547 template <typename CLASS, typename OutVector, typename InVector>
2548 inline
2549 void
2551 (void(CLASS::*function_pointer)(const MatrixFree<dim,Number> &,
2552  OutVector &,
2553  const InVector &,
2554  const std::pair<unsigned int,
2555  unsigned int> &),
2556  CLASS *owning_class,
2557  OutVector &dst,
2558  const InVector &src) const
2559 {
2560  // here, use std_cxx1x::bind to hand a function handler with the appropriate
2561  // argument to the other loop function
2562  std_cxx1x::function<void (const MatrixFree<dim,Number> &,
2563  OutVector &,
2564  const InVector &,
2565  const std::pair<unsigned int,
2566  unsigned int> &)>
2567  function = std_cxx1x::bind<void>(function_pointer,
2568  std_cxx1x::ref(*owning_class),
2569  std_cxx1x::_1,
2570  std_cxx1x::_2,
2571  std_cxx1x::_3,
2572  std_cxx1x::_4);
2573  cell_loop (function, dst, src);
2574 }
2575 
2576 
2577 #endif // ifndef DOXYGEN
2578 
2579 
2580 
2581 DEAL_II_NAMESPACE_CLOSE
2582 
2583 #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)
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:923
const IndexSet & get_ghost_set(const unsigned int fe_component=0) const
bool indices_are_initialized
Definition: matrix_free.h:928
static const unsigned int invalid_unsigned_int
Definition: types.h:191
bool at_irregular_cell(const unsigned int macro_cell_number) const
std::vector< unsigned int > constraint_pool_row_index
Definition: matrix_free.h:893
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
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:149
internal::MatrixFreeFunctions::MappingInfo< dim, Number > mapping_info
Definition: matrix_free.h:899
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)
DoFHandlers dof_handlers
Definition: matrix_free.h:873
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
Definition: matrix_free.h:887
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
Table< 4, internal::MatrixFreeFunctions::ShapeInfo< Number > > shape_info
Definition: matrix_free.h:904
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
Definition: matrix_free.h:912
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
Definition: matrix_free.h:933
unsigned int n_macro_cells() const
internal::MatrixFreeFunctions::SizeInfo size_info
Definition: matrix_free.h:918
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 print_memory_consumption(STREAM &out) const
#define Assert(cond, exc)
Definition: exceptions.h:299
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
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:205
const Quadrature< dim-1 > & get_face_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void update_ghost_values_start(const unsigned int communication_channel=0) const
Definition: hp.h:103
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:217
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
Definition: matrix_free.h:240
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:230
::ExceptionBase & ExcNotImplemented()
Definition: table.h:33
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
Definition: matrix_free.h:879
::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