Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
mapping_q.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2001 - 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
23
25
26#include <deal.II/fe/fe_q.h>
29
31
33
34#include <memory>
35#include <numeric>
36
38
39
40template <int dim, int spacedim>
42 : use_mapping_q1_on_current_cell(false)
43{}
44
45
46
47template <int dim, int spacedim>
48std::size_t
50{
51 return (
53 MemoryConsumption::memory_consumption(use_mapping_q1_on_current_cell) +
56}
57
58
59
60template <int dim, int spacedim>
61MappingQ<dim, spacedim>::MappingQ(const unsigned int degree,
62 const bool use_mapping_q_on_all_cells)
63 : polynomial_degree(degree)
64 ,
65
66 // see whether we want to use *this* mapping objects on *all* cells,
67 // or defer to an explicit Q1 mapping on interior cells. if
68 // degree==1, then we are already that Q1 mapping, so we don't need
69 // it; if dim!=spacedim, there is also no need for anything because
70 // we're most likely on a curved manifold
71 use_mapping_q_on_all_cells(degree == 1 || use_mapping_q_on_all_cells ||
72 (dim != spacedim))
73 ,
74 // create a Q1 mapping for use on interior cells (if necessary)
75 // or to create a good initial guess in transform_real_to_unit_cell()
76 q1_mapping(std::make_shared<MappingQGeneric<dim, spacedim>>(1))
77 ,
78
79 // create a Q_p mapping; if p=1, simply share the Q_1 mapping already
80 // created via the shared_ptr objects
81 qp_mapping(this->polynomial_degree > 1 ?
82 std::make_shared<const MappingQGeneric<dim, spacedim>>(degree) :
83 q1_mapping)
84{}
85
86
87
88template <int dim, int spacedim>
90 : polynomial_degree(mapping.polynomial_degree)
91 , use_mapping_q_on_all_cells(mapping.use_mapping_q_on_all_cells)
92{
93 // Note that we really do have to use clone() here, since mapping.q1_mapping
94 // may be MappingQ1Eulerian and mapping.qp_mapping may be MappingQEulerian.
95 std::shared_ptr<const Mapping<dim, spacedim>> other_q1_map =
96 mapping.q1_mapping->clone();
97 q1_mapping = std::dynamic_pointer_cast<const MappingQGeneric<dim, spacedim>>(
98 other_q1_map);
99 Assert(q1_mapping != nullptr, ExcInternalError());
100 Assert(q1_mapping->get_degree() == 1, ExcInternalError());
101
102 // Same as the other constructor: if possible reuse the Q1 mapping
103 if (this->polynomial_degree == 1)
104 {
105 qp_mapping = q1_mapping;
106 }
107 else
108 {
109 std::shared_ptr<const Mapping<dim, spacedim>> other_qp_map =
110 mapping.qp_mapping->clone();
111 qp_mapping =
112 std::dynamic_pointer_cast<const MappingQGeneric<dim, spacedim>>(
113 other_qp_map);
114 Assert(qp_mapping != nullptr, ExcInternalError());
115 }
116}
117
118
119
120template <int dim, int spacedim>
121unsigned int
123{
124 return polynomial_degree;
125}
126
127
128
129template <int dim, int spacedim>
130inline bool
132{
133 return true;
134}
135
136
137
138template <int dim, int spacedim>
141{
142 return (q1_mapping->requires_update_flags(in) |
143 qp_mapping->requires_update_flags(in));
144}
145
146
147
148template <int dim, int spacedim>
149std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
151 const Quadrature<dim> &quadrature) const
152{
153 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
154 std::make_unique<InternalData>();
155 auto &data = dynamic_cast<InternalData &>(*data_ptr);
156
157 // build the Q1 and Qp internal data objects in parallel
159 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
161 *qp_mapping,
162 update_flags,
163 quadrature);
164
166 data.mapping_q1_data = Utilities::dynamic_unique_cast<
168 std::move(q1_mapping->get_data(update_flags, quadrature)));
169
170 // wait for the task above to finish and use returned value
171 data.mapping_qp_data = Utilities::dynamic_unique_cast<
173 std::move(do_get_data.return_value()));
174 return data_ptr;
175}
176
177
178
179template <int dim, int spacedim>
180std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
182 const UpdateFlags update_flags,
183 const hp::QCollection<dim - 1> &quadrature) const
184{
185 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
186 std::make_unique<InternalData>();
187 auto &data = dynamic_cast<InternalData &>(*data_ptr);
188
189 std::unique_ptr<typename MappingQGeneric<dim, spacedim>::InternalDataBase> (
190 MappingQGeneric<dim, spacedim>::*mapping_get_face_data)(
191 const UpdateFlags, const hp::QCollection<dim - 1> &) const =
193
194 // build the Q1 and Qp internal data objects in parallel
196 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
197 do_get_data = Threads::new_task(mapping_get_face_data,
198 *qp_mapping,
199 update_flags,
200 quadrature);
201
203 data.mapping_q1_data = Utilities::dynamic_unique_cast<
205 std::move(q1_mapping->get_face_data(update_flags, quadrature)));
206
207 // wait for the task above to finish and use returned value
208 data.mapping_qp_data = Utilities::dynamic_unique_cast<
210 std::move(do_get_data.return_value()));
211 return data_ptr;
212}
213
214
215
216template <int dim, int spacedim>
217std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
219 const UpdateFlags update_flags,
220 const Quadrature<dim - 1> &quadrature) const
221{
222 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
223 std::make_unique<InternalData>();
224 auto &data = dynamic_cast<InternalData &>(*data_ptr);
225
226 // build the Q1 and Qp internal data objects in parallel
228 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
229 do_get_data =
231 *qp_mapping,
232 update_flags,
233 quadrature);
234
236 data.mapping_q1_data = Utilities::dynamic_unique_cast<
238 std::move(q1_mapping->get_subface_data(update_flags, quadrature)));
239
240 // wait for the task above to finish and use returned value
241 data.mapping_qp_data = Utilities::dynamic_unique_cast<
243 std::move(do_get_data.return_value()));
244 return data_ptr;
245}
246
247
248// Note that the CellSimilarity flag is modifiable, since MappingQ can need to
249// recalculate data even when cells are similar.
250template <int dim, int spacedim>
254 const CellSimilarity::Similarity cell_similarity,
255 const Quadrature<dim> & quadrature,
256 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
258 &output_data) const
259{
260 // convert data object to internal data for this class. fails with an
261 // exception if that is not possible
262 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
264 const InternalData &data = static_cast<const InternalData &>(internal_data);
265
266 // check whether this cell needs the full mapping or can be treated by a
267 // reduced Q1 mapping, e.g. if the cell is in the interior of the domain
268 data.use_mapping_q1_on_current_cell =
269 !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
270
271
272 // call the base class. we need to ensure that the flag indicating whether
273 // we can use some similarity has to be modified - for a general MappingQ,
274 // the data needs to be recomputed anyway since then the mapping changes the
275 // data. this needs to be known also for later operations, so modify the
276 // variable here. this also affects the calculation of the next cell -- if
277 // we use Q1 data on the next cell, the data will still be invalid.
278 const CellSimilarity::Similarity updated_cell_similarity =
279 ((data.use_mapping_q1_on_current_cell == false) &&
280 (this->polynomial_degree > 1) ?
282 cell_similarity);
283
284 // depending on the results above, decide whether the Q1 mapping or
285 // the Qp mapping needs to handle this cell
286 if (data.use_mapping_q1_on_current_cell)
287 q1_mapping->fill_fe_values(cell,
288 updated_cell_similarity,
289 quadrature,
290 *data.mapping_q1_data,
291 output_data);
292 else
293 qp_mapping->fill_fe_values(cell,
294 updated_cell_similarity,
295 quadrature,
296 *data.mapping_qp_data,
297 output_data);
298
299 return updated_cell_similarity;
300}
301
302
303
304template <int dim, int spacedim>
305void
308 const unsigned int face_no,
309 const hp::QCollection<dim - 1> & quadrature,
310 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
312 &output_data) const
313{
314 // convert data object to internal data for this class. fails with an
315 // exception if that is not possible
316 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
318 const InternalData &data = static_cast<const InternalData &>(internal_data);
319
320 // check whether this cell needs the full mapping or can be treated by a
321 // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
322 // domain. note that it is not sufficient to ask whether the present _face_
323 // is in the interior, as the mapping on the face depends on the mapping of
324 // the cell, which in turn depends on the fact whether _any_ of the faces of
325 // this cell is at the boundary, not only the present face
326 data.use_mapping_q1_on_current_cell =
327 !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
328
329 // depending on the results above, decide whether the Q1 mapping or
330 // the Qp mapping needs to handle this cell
331 if (data.use_mapping_q1_on_current_cell)
332 q1_mapping->fill_fe_face_values(
333 cell, face_no, quadrature, *data.mapping_q1_data, output_data);
334 else
335 qp_mapping->fill_fe_face_values(
336 cell, face_no, quadrature, *data.mapping_qp_data, output_data);
337}
338
339
340template <int dim, int spacedim>
341void
344 const unsigned int face_no,
345 const unsigned int subface_no,
346 const Quadrature<dim - 1> & quadrature,
347 const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
349 &output_data) const
350{
351 // convert data object to internal data for this class. fails with an
352 // exception if that is not possible
353 Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
355 const InternalData &data = static_cast<const InternalData &>(internal_data);
356
357 // check whether this cell needs the full mapping or can be treated by a
358 // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
359 // domain. note that it is not sufficient to ask whether the present _face_
360 // is in the interior, as the mapping on the face depends on the mapping of
361 // the cell, which in turn depends on the fact whether _any_ of the faces of
362 // this cell is at the boundary, not only the present face
363 data.use_mapping_q1_on_current_cell =
364 !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
365
366 // depending on the results above, decide whether the Q1 mapping or
367 // the Qp mapping needs to handle this cell
368 if (data.use_mapping_q1_on_current_cell)
369 q1_mapping->fill_fe_subface_values(cell,
370 face_no,
371 subface_no,
372 quadrature,
373 *data.mapping_q1_data,
374 output_data);
375 else
376 qp_mapping->fill_fe_subface_values(cell,
377 face_no,
378 subface_no,
379 quadrature,
380 *data.mapping_qp_data,
381 output_data);
382}
383
384
385
386template <int dim, int spacedim>
387void
389 const ArrayView<const Tensor<1, dim>> & input,
390 const MappingKind mapping_kind,
391 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
392 const ArrayView<Tensor<1, spacedim>> & output) const
393{
394 AssertDimension(input.size(), output.size());
395
396 const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
397 Assert(data != nullptr, ExcInternalError());
398
399 // check whether we should in fact work on the Q1 portion of it
400 if (data->use_mapping_q1_on_current_cell)
401 q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
402 else
403 // otherwise use the full mapping
404 qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
405}
406
407
408
409template <int dim, int spacedim>
410void
413 const MappingKind mapping_kind,
414 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
415 const ArrayView<Tensor<2, spacedim>> & output) const
416{
417 AssertDimension(input.size(), output.size());
418 Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
419 &mapping_data) != nullptr),
421 const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
422
423 // check whether we should in fact work on the Q1 portion of it
424 if (data->use_mapping_q1_on_current_cell)
425 q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
426 else
427 // otherwise use the full mapping
428 qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
429}
430
431
432
433template <int dim, int spacedim>
434void
436 const ArrayView<const Tensor<2, dim>> & input,
437 const MappingKind mapping_kind,
438 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
439 const ArrayView<Tensor<2, spacedim>> & output) const
440{
441 AssertDimension(input.size(), output.size());
442 Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
443 &mapping_data) != nullptr),
445 const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
446
447 // check whether we should in fact work on the Q1 portion of it
448 if (data->use_mapping_q1_on_current_cell)
449 q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
450 else
451 // otherwise use the full mapping
452 qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
453}
454
455
456
457template <int dim, int spacedim>
458void
461 const MappingKind mapping_kind,
462 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
463 const ArrayView<Tensor<3, spacedim>> & output) const
464{
465 AssertDimension(input.size(), output.size());
466 Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
467 &mapping_data) != nullptr),
469 const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
470
471 // check whether we should in fact work on the Q1 portion of it
472 if (data->use_mapping_q1_on_current_cell)
473 q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
474 else
475 // otherwise use the full mapping
476 qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
477}
478
479
480
481template <int dim, int spacedim>
482void
484 const ArrayView<const Tensor<3, dim>> & input,
485 const MappingKind mapping_kind,
486 const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
487 const ArrayView<Tensor<3, spacedim>> & output) const
488{
489 AssertDimension(input.size(), output.size());
490 Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
491 &mapping_data) != nullptr),
493 const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
494
495 // check whether we should in fact work on the Q1 portion of it
496 if (data->use_mapping_q1_on_current_cell)
497 q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
498 else
499 // otherwise use the full mapping
500 qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
501}
502
503
504
505template <int dim, int spacedim>
509 const Point<dim> & p) const
510{
511 // first see, whether we want to use a linear or a higher order
512 // mapping, then either use our own facilities or that of the Q1
513 // mapping we store
514 if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
515 return qp_mapping->transform_unit_to_real_cell(cell, p);
516 else
517 return q1_mapping->transform_unit_to_real_cell(cell, p);
518}
519
520
521
522template <int dim, int spacedim>
526 const Point<spacedim> & p) const
527{
528 if (cell->has_boundary_lines() || use_mapping_q_on_all_cells ||
529 (dim != spacedim))
530 return qp_mapping->transform_real_to_unit_cell(cell, p);
531 else
532 return q1_mapping->transform_real_to_unit_cell(cell, p);
533}
534
535
536
537template <int dim, int spacedim>
538std::unique_ptr<Mapping<dim, spacedim>>
540{
541 return std::make_unique<MappingQ<dim, spacedim>>(
543}
544
545
546
547template <int dim, int spacedim>
550 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
551{
552 if (cell->has_boundary_lines() || use_mapping_q_on_all_cells ||
553 (dim != spacedim))
555 qp_mapping->compute_mapping_support_points(cell));
556 else
557 return BoundingBox<spacedim>(q1_mapping->get_vertices(cell));
558}
559
560
561
562template <int dim, int spacedim>
563bool
565 const ReferenceCell &reference_cell) const
566{
567 Assert(dim == reference_cell.get_dimension(),
568 ExcMessage("The dimension of your mapping (" +
570 ") and the reference cell cell_type (" +
571 Utilities::to_string(reference_cell.get_dimension()) +
572 " ) do not agree."));
573
574 return reference_cell.is_hyper_cube();
575}
576
577
578
579// explicit instantiations
580#include "mapping_q.inst"
581
582
virtual std::size_t memory_consumption() const override
virtual void transform(const ArrayView< const Tensor< 1, dim > > &input, const MappingKind kind, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim > > &output) const override
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
std::shared_ptr< const MappingQGeneric< dim, spacedim > > q1_mapping
Definition: mapping_q.h:361
MappingQ(const unsigned int polynomial_degree, const bool use_mapping_q_on_all_cells=false)
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const override
virtual bool is_compatible_with(const ReferenceCell &reference_cell) const override
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
const unsigned int polynomial_degree
Definition: mapping_q.h:336
virtual bool preserves_vertex_locations() const override
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
const bool use_mapping_q_on_all_cells
Definition: mapping_q.h:342
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
virtual CellSimilarity::Similarity fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
std::shared_ptr< const MappingQGeneric< dim, spacedim > > qp_mapping
Definition: mapping_q.h:377
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const hp::QCollection< dim - 1 > &quadrature) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
unsigned int get_degree() const
Abstract base class for mapping classes.
Definition: mapping.h:304
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
UpdateFlags
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
MappingKind
Definition: mapping.h:65
Task< RT > new_task(const std::function< RT()> &function)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:1403
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
STL namespace.