Reference documentation for deal.II version Git 5a2787e538 2021-09-21 14:55:10 -0600
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mapping_q_cache.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 - 2021 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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/utilities.h>
19 
20 #include <deal.II/dofs/dof_tools.h>
21 
22 #include <deal.II/fe/fe_dgq.h>
23 #include <deal.II/fe/fe_nothing.h>
24 #include <deal.II/fe/fe_q.h>
25 #include <deal.II/fe/fe_tools.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/fe/mapping_q1.h>
29 
31 #include <deal.II/lac/la_vector.h>
34 
35 #include <functional>
36 
38 
39 template <int dim, int spacedim>
41  const unsigned int polynomial_degree)
42  : MappingQ<dim, spacedim>(polynomial_degree)
43  , uses_level_info(false)
44 {}
45 
46 
47 
48 template <int dim, int spacedim>
50  const MappingQCache<dim, spacedim> &mapping)
51  : MappingQ<dim, spacedim>(mapping)
54 {}
55 
56 
57 
58 template <int dim, int spacedim>
60 {
61  // When this object goes out of scope, we want the cache to get cleared and
62  // free its memory before the signal is disconnected in order to not work on
63  // invalid memory that has been left back by freeing an object of this
64  // class.
65  support_point_cache.reset();
66  clear_signal.disconnect();
67 }
68 
69 
70 
71 template <int dim, int spacedim>
72 std::unique_ptr<Mapping<dim, spacedim>>
74 {
75  return std::make_unique<MappingQCache<dim, spacedim>>(*this);
76 }
77 
78 
79 
80 template <int dim, int spacedim>
81 bool
83 {
84  return false;
85 }
86 
87 
88 
89 template <int dim, int spacedim>
90 void
92  const Mapping<dim, spacedim> & mapping,
94 {
95  // FE and FEValues in the case they are needed
98  fe_values_all;
99 
100  this->initialize(
101  triangulation,
102  [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell) {
103  const auto mapping_q =
104  dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping);
105  if (mapping_q != nullptr && this->get_degree() == mapping_q->get_degree())
106  {
107  return mapping_q->compute_mapping_support_points(cell);
108  }
109  else
110  {
111  // get FEValues (thread-safe); in the case that this thread has not
112  // created a an FEValues object yet, this helper-function also
113  // creates one with the right quadrature rule
114  auto &fe_values = fe_values_all.get();
115  if (fe_values.get() == nullptr)
116  {
117  QGaussLobatto<dim> quadrature_gl(this->polynomial_degree + 1);
118 
119  std::vector<Point<dim>> quadrature_points;
120  for (const auto i :
121  FETools::hierarchic_to_lexicographic_numbering<dim>(
122  this->polynomial_degree))
123  quadrature_points.push_back(quadrature_gl.point(i));
124  Quadrature<dim> quadrature(quadrature_points);
125 
126  fe_values = std::make_unique<FEValues<dim, spacedim>>(
127  mapping, fe, quadrature, update_quadrature_points);
128  }
129 
130  fe_values->reinit(cell);
131  return fe_values->get_quadrature_points();
132  }
133  });
134 }
135 
136 
137 
138 template <int dim, int spacedim>
139 void
142  const MappingQ<dim, spacedim> & mapping)
143 {
144  this->initialize(mapping, triangulation);
145 }
146 
147 
148 
149 template <int dim, int spacedim>
150 void
153  const std::function<std::vector<Point<spacedim>>(
155  &compute_points_on_cell)
156 {
157  clear_signal.disconnect();
158  clear_signal = triangulation.signals.any_change.connect(
159  [&]() -> void { this->support_point_cache.reset(); });
160 
162  std::make_shared<std::vector<std::vector<std::vector<Point<spacedim>>>>>(
163  triangulation.n_levels());
164  for (unsigned int l = 0; l < triangulation.n_levels(); ++l)
165  (*support_point_cache)[l].resize(triangulation.n_raw_cells(l));
166 
168  triangulation.begin(),
169  triangulation.end(),
170  [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
171  void *,
172  void *) {
173  (*support_point_cache)[cell->level()][cell->index()] =
174  compute_points_on_cell(cell);
176  (*support_point_cache)[cell->level()][cell->index()].size(),
177  Utilities::pow(this->get_degree() + 1, dim));
178  },
179  /* copier */ std::function<void(void *)>(),
180  /* scratch_data */ nullptr,
181  /* copy_data */ nullptr,
183  /* chunk_size = */ 1);
184 
185  uses_level_info = true;
186 }
187 
188 
189 
190 template <int dim, int spacedim>
191 void
193  const Mapping<dim, spacedim> & mapping,
195  const std::function<Point<spacedim>(
197  const Point<spacedim> &)> & transformation_function,
198  const bool function_describes_relative_displacement)
199 {
200  // FE and FEValues in the case they are needed
203  fe_values_all;
204 
205  this->initialize(
206  tria,
207  [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell) {
208  std::vector<Point<spacedim>> points;
209 
210  const auto mapping_q =
211  dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping);
212 
213  if (mapping_q != nullptr && this->get_degree() == mapping_q->get_degree())
214  {
215  points = mapping_q->compute_mapping_support_points(cell);
216  }
217  else
218  {
219  // get FEValues (thread-safe); in the case that this thread has not
220  // created a an FEValues object yet, this helper-function also
221  // creates one with the right quadrature rule
222  auto &fe_values = fe_values_all.get();
223  if (fe_values.get() == nullptr)
224  {
225  QGaussLobatto<dim> quadrature_gl(this->polynomial_degree + 1);
226 
227  std::vector<Point<dim>> quadrature_points;
228  for (const auto i :
229  FETools::hierarchic_to_lexicographic_numbering<dim>(
230  this->polynomial_degree))
231  quadrature_points.push_back(quadrature_gl.point(i));
232  Quadrature<dim> quadrature(quadrature_points);
233 
234  fe_values = std::make_unique<FEValues<dim, spacedim>>(
235  mapping, fe, quadrature, update_quadrature_points);
236  }
237 
238  fe_values->reinit(cell);
239  points = fe_values->get_quadrature_points();
240  }
241 
242  for (auto &p : points)
243  if (function_describes_relative_displacement)
244  p += transformation_function(cell, p);
245  else
246  p = transformation_function(cell, p);
247 
248  return points;
249  });
250 
251  uses_level_info = true;
252 }
253 
254 
255 
256 template <int dim, int spacedim>
257 void
259  const Mapping<dim, spacedim> & mapping,
260  const Triangulation<dim, spacedim> &tria,
261  const Function<spacedim> & transformation_function,
262  const bool function_describes_relative_displacement)
263 {
264  AssertDimension(transformation_function.n_components, spacedim);
265 
266  this->initialize(
267  mapping,
268  tria,
269  [&](const auto &, const auto &point) {
270  Point<spacedim> new_point;
271  for (int c = 0; c < spacedim; ++c)
272  new_point[c] = transformation_function.value(point, c);
273  return new_point;
274  },
275  function_describes_relative_displacement);
276 
277  uses_level_info = true;
278 }
279 
280 
281 
282 namespace
283 {
284  template <typename VectorType>
285  void
286  copy_locally_owned_data_from(
287  const VectorType &vector,
289  &vector_ghosted)
290  {
292  temp.reinit(vector.locally_owned_elements());
293  temp.import(vector, VectorOperation::insert);
294  vector_ghosted.import(temp, VectorOperation::insert);
295  }
296 } // namespace
297 
298 
299 
300 template <int dim, int spacedim>
301 template <typename VectorType>
302 void
304  const Mapping<dim, spacedim> & mapping,
305  const DoFHandler<dim, spacedim> &dof_handler,
306  const VectorType & vector,
307  const bool vector_describes_relative_displacement)
308 {
309  AssertDimension(dof_handler.get_fe_collection().size(), 1);
310  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
311  AssertDimension(fe.n_base_elements(), 1);
312  AssertDimension(fe.element_multiplicity(0), spacedim);
313 
314  const unsigned int is_fe_q =
315  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe.base_element(0)) != nullptr;
316  const unsigned int is_fe_dgq =
317  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe.base_element(0)) != nullptr;
318 
321  FETools::hierarchic_to_lexicographic_numbering<spacedim>(
322  this->get_degree()));
323 
324  // Step 1: copy global vector so that the ghost values are such that the
325  // cache can be set up for all ghost cells
327  vector_ghosted;
328  IndexSet locally_relevant_dofs;
329  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
330  vector_ghosted.reinit(dof_handler.locally_owned_dofs(),
331  locally_relevant_dofs,
332  dof_handler.get_communicator());
333  copy_locally_owned_data_from(vector, vector_ghosted);
334  vector_ghosted.update_ghost_values();
335 
336  // FE and FEValues in the case they are needed
337  FE_Nothing<dim, spacedim> fe_nothing;
339  fe_values_all;
340 
341  // Interpolation of values is needed if we cannot just read off locations
342  // from the solution vectors (as in the case of FE_Q and FE_DGQ with the
343  // same polynomial degree as this class has).
344  const bool interpolation_of_values_is_needed =
345  ((is_fe_q || is_fe_dgq) && fe.degree == this->get_degree()) == false;
346 
347  // Step 2: loop over all cells
348  this->initialize(
349  dof_handler.get_triangulation(),
350  [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell_tria)
351  -> std::vector<Point<spacedim>> {
352  const bool is_active_non_artificial_cell =
353  (cell_tria->is_active() == true) &&
354  (cell_tria->is_artificial() == false);
355 
356  const typename DoFHandler<dim, spacedim>::cell_iterator cell_dofs(
357  &cell_tria->get_triangulation(),
358  cell_tria->level(),
359  cell_tria->index(),
360  &dof_handler);
361 
362  const auto mapping_q =
363  dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping);
364 
365  // Step 2a) set up and reinit FEValues (if needed)
366  if (
367  ((vector_describes_relative_displacement ||
368  (is_active_non_artificial_cell == false)) &&
369  ((mapping_q != nullptr &&
370  this->get_degree() == mapping_q->get_degree()) ==
371  false)) /*condition 1: points need to be computed via FEValues*/
372  ||
373  (is_active_non_artificial_cell && interpolation_of_values_is_needed) /*condition 2: interpolation of values is needed*/)
374  {
375  // get FEValues (thread-safe); in the case that this thread has
376  // not created a an FEValues object yet, this helper-function also
377  // creates one with the right quadrature rule
378  auto &fe_values = fe_values_all.get();
379  if (fe_values.get() == nullptr)
380  {
381  QGaussLobatto<dim> quadrature_gl(this->polynomial_degree + 1);
382 
383  std::vector<Point<dim>> quadrature_points;
384  for (const auto i :
385  FETools::hierarchic_to_lexicographic_numbering<dim>(
386  this->polynomial_degree))
387  quadrature_points.push_back(quadrature_gl.point(i));
388  Quadrature<dim> quadrature(quadrature_points);
389 
390  fe_values = std::make_unique<FEValues<dim, spacedim>>(
391  mapping,
392  interpolation_of_values_is_needed ?
393  fe :
394  static_cast<const FiniteElement<dim, spacedim> &>(fe_nothing),
395  quadrature,
397  }
398 
399  if (interpolation_of_values_is_needed)
400  fe_values->reinit(cell_dofs);
401  else
402  fe_values->reinit(cell_tria);
403  }
404 
405  std::vector<Point<spacedim>> result;
406 
407  // Step 2b) read of quadrature points in the relative displacement case
408  // note: we also take this path for non-active or artificial cells so that
409  // these cells are filled with some useful data
410  if (vector_describes_relative_displacement ||
411  is_active_non_artificial_cell == false)
412  {
413  if (mapping_q != nullptr &&
414  this->get_degree() == mapping_q->get_degree())
415  result = mapping_q->compute_mapping_support_points(cell_tria);
416  else
417  result = fe_values_all.get()->get_quadrature_points();
418 
419  // for non-active or artificial cells we are done here and return
420  // the absolute positions, since the provided vector cannot contain
421  // any useful information for these cells
422  if (is_active_non_artificial_cell == false)
423  return result;
424  }
425  else
426  {
427  result.resize(
428  Utilities::pow<unsigned int>(this->get_degree() + 1, dim));
429  }
430 
431  // Step 2c) read global vector and adjust points accordingly
432  if (interpolation_of_values_is_needed == false)
433  {
434  // case 1: FE_Q or FE_DGQ with same degree as this class has; this
435  // is the simple case since no interpolation is needed
436  std::vector<types::global_dof_index> dof_indices(
437  fe.n_dofs_per_cell());
438  cell_dofs->get_dof_indices(dof_indices);
439 
440  for (unsigned int i = 0; i < dof_indices.size(); ++i)
441  {
442  const auto id = fe.system_to_component_index(i);
443 
444  if (is_fe_q)
445  {
446  // case 1a: FE_Q
447  if (vector_describes_relative_displacement)
448  result[id.second][id.first] +=
449  vector_ghosted(dof_indices[i]);
450  else
451  result[id.second][id.first] =
452  vector_ghosted(dof_indices[i]);
453  }
454  else
455  {
456  // case 1b: FE_DGQ
457  if (vector_describes_relative_displacement)
458  result[lexicographic_to_hierarchic_numbering[id.second]]
459  [id.first] += vector_ghosted(dof_indices[i]);
460  else
461  result[lexicographic_to_hierarchic_numbering[id.second]]
462  [id.first] = vector_ghosted(dof_indices[i]);
463  }
464  }
465  }
466  else
467  {
468  // case 2: general case; interpolation is needed
469  // note: the following code could be optimized for tensor-product
470  // elements via application of sum factorization as is done on
471  // MatrixFree/FEEvaluation
472  auto &fe_values = fe_values_all.get();
473 
474  std::vector<Vector<typename VectorType::value_type>> values(
475  fe_values->n_quadrature_points,
477 
478  fe_values->get_function_values(vector_ghosted, values);
479 
480  for (unsigned int q = 0; q < fe_values->n_quadrature_points; ++q)
481  for (unsigned int c = 0; c < spacedim; ++c)
482  if (vector_describes_relative_displacement)
483  result[q][c] += values[q][c];
484  else
485  result[q][c] = values[q][c];
486  }
487 
488  return result;
489  });
490 
491  uses_level_info = false;
492 }
493 
494 
495 
496 template <int dim, int spacedim>
497 template <typename VectorType>
498 void
500  const Mapping<dim, spacedim> & mapping,
501  const DoFHandler<dim, spacedim> &dof_handler,
502  const MGLevelObject<VectorType> &vectors,
503  const bool vector_describes_relative_displacement)
504 {
505  AssertDimension(dof_handler.get_fe_collection().size(), 1);
506  const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
507  AssertDimension(fe.n_base_elements(), 1);
508  AssertDimension(fe.element_multiplicity(0), spacedim);
509  AssertDimension(0, vectors.min_level());
510  AssertDimension(dof_handler.get_triangulation().n_global_levels() - 1,
511  vectors.max_level());
512 
513  const unsigned int is_fe_q =
514  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe.base_element(0)) != nullptr;
515  const unsigned int is_fe_dgq =
516  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe.base_element(0)) != nullptr;
517 
520  FETools::hierarchic_to_lexicographic_numbering<spacedim>(
521  this->get_degree()));
522 
523  // Step 1: copy global vector so that the ghost values are such that the
524  // cache can be set up for all ghost cells
527  vectors_ghosted(vectors.min_level(), vectors.max_level());
528 
529  for (unsigned int l = vectors.min_level(); l <= vectors.max_level(); ++l)
530  {
531  IndexSet locally_relevant_dofs;
533  l,
534  locally_relevant_dofs);
535  vectors_ghosted[l].reinit(dof_handler.locally_owned_mg_dofs(l),
536  locally_relevant_dofs,
537  dof_handler.get_communicator());
538  copy_locally_owned_data_from(vectors[l], vectors_ghosted[l]);
539  vectors_ghosted[l].update_ghost_values();
540  }
541 
542  // FE and FEValues in the case they are needed
543  FE_Nothing<dim, spacedim> fe_nothing;
545  fe_values_all;
546 
547  // Interpolation of values is needed if we cannot just read off locations
548  // from the solution vectors (as in the case of FE_Q and FE_DGQ with the
549  // same polynomial degree as this class has).
550  const bool interpolation_of_values_is_needed =
551  ((is_fe_q || is_fe_dgq) && fe.degree == this->get_degree()) == false;
552 
553  // Step 2: loop over all cells
554  this->initialize(
555  dof_handler.get_triangulation(),
556  [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell_tria)
557  -> std::vector<Point<spacedim>> {
558  const bool is_non_artificial_cell =
559  cell_tria->level_subdomain_id() != numbers::artificial_subdomain_id;
560 
561  const typename DoFHandler<dim, spacedim>::level_cell_iterator cell_dofs(
562  &cell_tria->get_triangulation(),
563  cell_tria->level(),
564  cell_tria->index(),
565  &dof_handler);
566 
567  const auto mapping_q =
568  dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping);
569 
570  // Step 2a) set up and reinit FEValues (if needed)
571  if (
572  ((vector_describes_relative_displacement ||
573  (is_non_artificial_cell == false)) &&
574  ((mapping_q != nullptr &&
575  this->get_degree() == mapping_q->get_degree()) ==
576  false)) /*condition 1: points need to be computed via FEValues*/
577  ||
578  (is_non_artificial_cell == true && interpolation_of_values_is_needed) /*condition 2: interpolation of values is needed*/)
579  {
580  // get FEValues (thread-safe); in the case that this thread has
581  // not created a an FEValues object yet, this helper-function also
582  // creates one with the right quadrature rule
583  auto &fe_values = fe_values_all.get();
584  if (fe_values.get() == nullptr)
585  {
586  QGaussLobatto<dim> quadrature_gl(this->polynomial_degree + 1);
587 
588  std::vector<Point<dim>> quadrature_points;
589  for (const auto i :
590  FETools::hierarchic_to_lexicographic_numbering<dim>(
591  this->polynomial_degree))
592  quadrature_points.push_back(quadrature_gl.point(i));
593  Quadrature<dim> quadrature(quadrature_points);
594 
595  fe_values = std::make_unique<FEValues<dim, spacedim>>(
596  mapping,
597  interpolation_of_values_is_needed ?
598  fe :
599  static_cast<const FiniteElement<dim, spacedim> &>(fe_nothing),
600  quadrature,
602  }
603 
604  if (interpolation_of_values_is_needed)
605  fe_values->reinit(cell_dofs);
606  else
607  fe_values->reinit(cell_tria);
608  }
609 
610  std::vector<Point<spacedim>> result;
611 
612  // Step 2b) read of quadrature points in the relative displacement case
613  // note: we also take this path for non-active or artificial cells so that
614  // these cells are filled with some useful data
615  if (vector_describes_relative_displacement ||
616  (is_non_artificial_cell == false))
617  {
618  if (mapping_q != nullptr &&
619  this->get_degree() == mapping_q->get_degree())
620  result = mapping_q->compute_mapping_support_points(cell_tria);
621  else
622  result = fe_values_all.get()->get_quadrature_points();
623 
624  // for non-active or artificial cells we are done here and return
625  // the absolute positions, since the provided vector cannot contain
626  // any useful information for these cells
627  if (is_non_artificial_cell == false)
628  return result;
629  }
630  else
631  {
632  result.resize(
633  Utilities::pow<unsigned int>(this->get_degree() + 1, dim));
634  }
635 
636  // Step 2c) read global vector and adjust points accordingly
637  if (interpolation_of_values_is_needed == false)
638  {
639  // case 1: FE_Q or FE_DGQ with same degree as this class has; this
640  // is the simple case since no interpolation is needed
641  std::vector<types::global_dof_index> dof_indices(
642  fe.n_dofs_per_cell());
643  cell_dofs->get_mg_dof_indices(dof_indices);
644 
645  for (unsigned int i = 0; i < dof_indices.size(); ++i)
646  {
647  const auto id = fe.system_to_component_index(i);
648 
649  if (is_fe_q)
650  {
651  // case 1a: FE_Q
652  if (vector_describes_relative_displacement)
653  result[id.second][id.first] +=
654  vectors_ghosted[cell_tria->level()](dof_indices[i]);
655  else
656  result[id.second][id.first] =
657  vectors_ghosted[cell_tria->level()](dof_indices[i]);
658  }
659  else
660  {
661  // case 1b: FE_DGQ
662  if (vector_describes_relative_displacement)
663  result[lexicographic_to_hierarchic_numbering[id.second]]
664  [id.first] +=
665  vectors_ghosted[cell_tria->level()](dof_indices[i]);
666  else
667  result[lexicographic_to_hierarchic_numbering[id.second]]
668  [id.first] =
669  vectors_ghosted[cell_tria->level()](dof_indices[i]);
670  }
671  }
672  }
673  else
674  {
675  // case 2: general case; interpolation is needed
676  // note: the following code could be optimized for tensor-product
677  // elements via application of sum factorization as is done on
678  // MatrixFree/FEEvaluation
679  auto &fe_values = fe_values_all.get();
680 
681  std::vector<types::global_dof_index> dof_indices(
682  fe.n_dofs_per_cell());
683  cell_dofs->get_mg_dof_indices(dof_indices);
684 
685  std::vector<typename VectorType::value_type> dof_values(
686  fe.n_dofs_per_cell());
687 
688  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
689  dof_values[i] = vectors_ghosted[cell_tria->level()](dof_indices[i]);
690 
691  for (unsigned int c = 0; c < spacedim; ++c)
692  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
693  for (unsigned int q = 0; q < fe_values->n_quadrature_points; ++q)
694  if (vector_describes_relative_displacement == false && i == 0)
695  result[q][c] =
696  dof_values[i] * fe_values->shape_value_component(i, q, c);
697  else
698  result[q][c] +=
699  dof_values[i] * fe_values->shape_value_component(i, q, c);
700  }
701 
702  return result;
703  });
704 
705  uses_level_info = true;
706 }
707 
708 
709 
710 template <int dim, int spacedim>
711 std::size_t
713 {
714  if (support_point_cache.get() != nullptr)
715  return sizeof(*this) +
717  else
718  return sizeof(*this);
719 }
720 
721 
722 
723 template <int dim, int spacedim>
724 std::vector<Point<spacedim>>
726  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
727 {
728  Assert(support_point_cache.get() != nullptr,
729  ExcMessage("Must call MappingQCache::initialize() before "
730  "using it or after mesh has changed!"));
731 
732  Assert(uses_level_info || cell->is_active(), ExcInternalError());
733 
734  AssertIndexRange(cell->level(), support_point_cache->size());
735  AssertIndexRange(cell->index(), (*support_point_cache)[cell->level()].size());
736  return (*support_point_cache)[cell->level()][cell->index()];
737 }
738 
739 
740 
741 template <int dim, int spacedim>
742 boost::container::small_vector<Point<spacedim>,
745  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
746 {
747  Assert(support_point_cache.get() != nullptr,
748  ExcMessage("Must call MappingQCache::initialize() before "
749  "using it or after mesh has changed!"));
750 
751  Assert(uses_level_info || cell->is_active(), ExcInternalError());
752 
753  AssertIndexRange(cell->level(), support_point_cache->size());
754  AssertIndexRange(cell->index(), (*support_point_cache)[cell->level()].size());
755  const auto ptr = (*support_point_cache)[cell->level()][cell->index()].begin();
756  return boost::container::small_vector<Point<spacedim>,
758  ptr, ptr + cell->n_vertices());
759 }
760 
761 
762 
763 //--------------------------- Explicit instantiations -----------------------
764 #include "mapping_q_cache.inst"
765 
766 
unsigned int get_degree() const
Definition: mapping_q.cc:444
unsigned int max_level() const
Shape function values.
virtual bool preserves_vertex_locations() const override
const Triangulation< dim, spacedim > & get_triangulation() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1653
const unsigned int n_components
Definition: function.h:164
A class that provides a separate storage location on each thread that accesses the object...
unsigned int n_raw_cells(const unsigned int level) const
Definition: tria.cc:13832
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:451
std::vector< unsigned int > lexicographic_to_hierarchic_numbering(unsigned int degree)
const ::Triangulation< dim, spacedim > & tria
Definition: fe_dgq.h:110
#define AssertIndexRange(index, range)
Definition: exceptions.h:1718
MPI_Comm get_communicator() const
const Point< dim > & point(const unsigned int i) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
void import(const Vector< Number, MemorySpace2 > &src, VectorOperation::values operation)
Transformed quadrature points.
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:13047
unsigned int min_level() const
unsigned int n_levels() const
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
virtual void reinit(const size_type size, const bool omit_zeroing_entries=false)
void extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1175
std::size_t memory_consumption() const
cell_iterator end() const
Definition: tria.cc:13158
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
void initialize(const Mapping< dim, spacedim > &mapping, const Triangulation< dim, spacedim > &triangulation)
static ::ExceptionBase & ExcMessage(std::string arg1)
Definition: fe_q.h:548
#define Assert(cond, exc)
Definition: exceptions.h:1461
Signals signals
Definition: tria.h:2289
Abstract base class for mapping classes.
Definition: mapping.h:303
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
const unsigned int polynomial_degree
Definition: mapping_q.h:628
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:502
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1486
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:185
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping_q.cc:1680
boost::signals2::connection clear_signal
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
MappingQCache(const unsigned int polynomial_degree)
void import(const ::Vector< Number > &vec, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1337
static unsigned int n_threads()
virtual std::vector< Point< spacedim > > compute_mapping_support_points(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
void extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1133
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
const IndexSet & locally_owned_dofs() const
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
std::shared_ptr< std::vector< std::vector< std::vector< Point< spacedim > > > > > support_point_cache
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()