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