Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mapping_q.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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 
16 #include <deal.II/base/array_view.h>
17 #include <deal.II/base/memory_consumption.h>
18 #include <deal.II/base/polynomial.h>
19 #include <deal.II/base/quadrature.h>
20 #include <deal.II/base/quadrature_lib.h>
21 #include <deal.II/base/std_cxx14/memory.h>
22 #include <deal.II/base/tensor_product_polynomials.h>
23 #include <deal.II/base/utilities.h>
24 
25 #include <deal.II/dofs/dof_accessor.h>
26 
27 #include <deal.II/fe/fe_q.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping_q.h>
30 
31 #include <deal.II/grid/tria_iterator.h>
32 
33 #include <deal.II/lac/full_matrix.h>
34 
35 #include <memory>
36 #include <numeric>
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 
41 template <int dim, int spacedim>
43  : use_mapping_q1_on_current_cell(false)
44 {}
45 
46 
47 
48 template <int dim, int spacedim>
49 std::size_t
51 {
52  return (
54  MemoryConsumption::memory_consumption(use_mapping_q1_on_current_cell) +
55  MemoryConsumption::memory_consumption(mapping_q1_data) +
56  MemoryConsumption::memory_consumption(mapping_qp_data));
57 }
58 
59 
60 
61 template <int dim, int spacedim>
62 MappingQ<dim, spacedim>::MappingQ(const unsigned int degree,
63  const bool use_mapping_q_on_all_cells)
64  : polynomial_degree(degree)
65  ,
66 
67  // see whether we want to use *this* mapping objects on *all* cells,
68  // or defer to an explicit Q1 mapping on interior cells. if
69  // degree==1, then we are already that Q1 mapping, so we don't need
70  // it; if dim!=spacedim, there is also no need for anything because
71  // we're most likely on a curved manifold
73  (dim != spacedim))
74  ,
75  // create a Q1 mapping for use on interior cells (if necessary)
76  // or to create a good initial guess in transform_real_to_unit_cell()
77  q1_mapping(std::make_shared<MappingQGeneric<dim, spacedim>>(1))
78  ,
79 
80  // create a Q_p mapping; if p=1, simply share the Q_1 mapping already
81  // created via the shared_ptr objects
82  qp_mapping(this->polynomial_degree > 1 ?
83  std::make_shared<const MappingQGeneric<dim, spacedim>>(degree) :
84  q1_mapping)
85 {}
86 
87 
88 
89 template <int dim, int spacedim>
91  : polynomial_degree(mapping.polynomial_degree)
92  , use_mapping_q_on_all_cells(mapping.use_mapping_q_on_all_cells)
93 {
94  // Note that we really do have to use clone() here, since mapping.q1_mapping
95  // may be MappingQ1Eulerian and mapping.qp_mapping may be MappingQEulerian.
96  std::shared_ptr<const Mapping<dim, spacedim>> other_q1_map =
97  mapping.q1_mapping->clone();
98  q1_mapping = std::dynamic_pointer_cast<const MappingQGeneric<dim, spacedim>>(
99  other_q1_map);
100  Assert(q1_mapping != nullptr, ExcInternalError());
101  Assert(q1_mapping->get_degree() == 1, ExcInternalError());
102 
103  // Same as the other constructor: if possible reuse the Q1 mapping
104  if (this->polynomial_degree == 1)
105  {
107  }
108  else
109  {
110  std::shared_ptr<const Mapping<dim, spacedim>> other_qp_map =
111  mapping.qp_mapping->clone();
112  qp_mapping =
113  std::dynamic_pointer_cast<const MappingQGeneric<dim, spacedim>>(
114  other_qp_map);
115  Assert(qp_mapping != nullptr, ExcInternalError());
116  }
117 }
118 
119 
120 
121 template <int dim, int spacedim>
122 unsigned int
124 {
125  return polynomial_degree;
126 }
127 
128 
129 
130 template <int dim, int spacedim>
131 inline bool
133 {
134  return true;
135 }
136 
137 
138 
139 template <int dim, int spacedim>
142 {
143  return (q1_mapping->requires_update_flags(in) |
144  qp_mapping->requires_update_flags(in));
145 }
146 
147 
148 
149 template <int dim, int spacedim>
150 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
152  const Quadrature<dim> &quadrature) const
153 {
154  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
155  std_cxx14::make_unique<InternalData>();
156  auto &data = dynamic_cast<InternalData &>(*data_ptr);
157 
158  // build the Q1 and Qp internal data objects in parallel
160  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
162  *qp_mapping,
163  update_flags,
164  quadrature);
165 
166  if (!use_mapping_q_on_all_cells)
167  data.mapping_q1_data = Utilities::dynamic_unique_cast<
169  std::move(q1_mapping->get_data(update_flags, quadrature)));
170 
171  // wait for the task above to finish and use returned value
172  data.mapping_qp_data = Utilities::dynamic_unique_cast<
174  std::move(do_get_data.return_value()));
175  return data_ptr;
176 }
177 
178 
179 
180 template <int dim, int spacedim>
181 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
183  const UpdateFlags update_flags,
184  const Quadrature<dim - 1> &quadrature) const
185 {
186  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
187  std_cxx14::make_unique<InternalData>();
188  auto &data = dynamic_cast<InternalData &>(*data_ptr);
189 
190  // build the Q1 and Qp internal data objects in parallel
192  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
193  do_get_data =
195  *qp_mapping,
196  update_flags,
197  quadrature);
198 
199  if (!use_mapping_q_on_all_cells)
200  data.mapping_q1_data = Utilities::dynamic_unique_cast<
202  std::move(q1_mapping->get_face_data(update_flags, quadrature)));
203 
204  // wait for the task above to finish and use returned value
205  data.mapping_qp_data = Utilities::dynamic_unique_cast<
207  std::move(do_get_data.return_value()));
208  return data_ptr;
209 }
210 
211 
212 
213 template <int dim, int spacedim>
214 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
216  const UpdateFlags update_flags,
217  const Quadrature<dim - 1> &quadrature) const
218 {
219  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
220  std_cxx14::make_unique<InternalData>();
221  auto &data = dynamic_cast<InternalData &>(*data_ptr);
222 
223  // build the Q1 and Qp internal data objects in parallel
225  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
226  do_get_data =
228  *qp_mapping,
229  update_flags,
230  quadrature);
231 
232  if (!use_mapping_q_on_all_cells)
233  data.mapping_q1_data = Utilities::dynamic_unique_cast<
235  std::move(q1_mapping->get_subface_data(update_flags, quadrature)));
236 
237  // wait for the task above to finish and use returned value
238  data.mapping_qp_data = Utilities::dynamic_unique_cast<
240  std::move(do_get_data.return_value()));
241  return data_ptr;
242 }
243 
244 
245 // Note that the CellSimilarity flag is modifiable, since MappingQ can need to
246 // recalculate data even when cells are similar.
247 template <int dim, int spacedim>
250  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
251  const CellSimilarity::Similarity cell_similarity,
252  const Quadrature<dim> & quadrature,
253  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
255  &output_data) const
256 {
257  // convert data object to internal data for this class. fails with an
258  // exception if that is not possible
259  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
260  ExcInternalError());
261  const InternalData &data = static_cast<const InternalData &>(internal_data);
262 
263  // check whether this cell needs the full mapping or can be treated by a
264  // reduced Q1 mapping, e.g. if the cell is in the interior of the domain
265  data.use_mapping_q1_on_current_cell =
266  !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
267 
268 
269  // call the base class. we need to ensure that the flag indicating whether
270  // we can use some similarity has to be modified - for a general MappingQ,
271  // the data needs to be recomputed anyway since then the mapping changes the
272  // data. this needs to be known also for later operations, so modify the
273  // variable here. this also affects the calculation of the next cell -- if
274  // we use Q1 data on the next cell, the data will still be invalid.
275  const CellSimilarity::Similarity updated_cell_similarity =
276  ((data.use_mapping_q1_on_current_cell == false) &&
277  (this->polynomial_degree > 1) ?
279  cell_similarity);
280 
281  // depending on the results above, decide whether the Q1 mapping or
282  // the Qp mapping needs to handle this cell
283  if (data.use_mapping_q1_on_current_cell)
284  q1_mapping->fill_fe_values(cell,
285  updated_cell_similarity,
286  quadrature,
287  *data.mapping_q1_data,
288  output_data);
289  else
290  qp_mapping->fill_fe_values(cell,
291  updated_cell_similarity,
292  quadrature,
293  *data.mapping_qp_data,
294  output_data);
295 
296  return updated_cell_similarity;
297 }
298 
299 
300 
301 template <int dim, int spacedim>
302 void
304  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
305  const unsigned int face_no,
306  const Quadrature<dim - 1> & quadrature,
307  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
309  &output_data) const
310 {
311  // convert data object to internal data for this class. fails with an
312  // exception if that is not possible
313  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
314  ExcInternalError());
315  const InternalData &data = static_cast<const InternalData &>(internal_data);
316 
317  // check whether this cell needs the full mapping or can be treated by a
318  // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
319  // domain. note that it is not sufficient to ask whether the present _face_
320  // is in the interior, as the mapping on the face depends on the mapping of
321  // the cell, which in turn depends on the fact whether _any_ of the faces of
322  // this cell is at the boundary, not only the present face
323  data.use_mapping_q1_on_current_cell =
324  !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
325 
326  // depending on the results above, decide whether the Q1 mapping or
327  // the Qp mapping needs to handle this cell
328  if (data.use_mapping_q1_on_current_cell)
329  q1_mapping->fill_fe_face_values(
330  cell, face_no, quadrature, *data.mapping_q1_data, output_data);
331  else
332  qp_mapping->fill_fe_face_values(
333  cell, face_no, quadrature, *data.mapping_qp_data, output_data);
334 }
335 
336 
337 template <int dim, int spacedim>
338 void
340  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
341  const unsigned int face_no,
342  const unsigned int subface_no,
343  const Quadrature<dim - 1> & quadrature,
344  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
346  &output_data) const
347 {
348  // convert data object to internal data for this class. fails with an
349  // exception if that is not possible
350  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
351  ExcInternalError());
352  const InternalData &data = static_cast<const InternalData &>(internal_data);
353 
354  // check whether this cell needs the full mapping or can be treated by a
355  // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
356  // domain. note that it is not sufficient to ask whether the present _face_
357  // is in the interior, as the mapping on the face depends on the mapping of
358  // the cell, which in turn depends on the fact whether _any_ of the faces of
359  // this cell is at the boundary, not only the present face
360  data.use_mapping_q1_on_current_cell =
361  !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
362 
363  // depending on the results above, decide whether the Q1 mapping or
364  // the Qp mapping needs to handle this cell
365  if (data.use_mapping_q1_on_current_cell)
366  q1_mapping->fill_fe_subface_values(cell,
367  face_no,
368  subface_no,
369  quadrature,
370  *data.mapping_q1_data,
371  output_data);
372  else
373  qp_mapping->fill_fe_subface_values(cell,
374  face_no,
375  subface_no,
376  quadrature,
377  *data.mapping_qp_data,
378  output_data);
379 }
380 
381 
382 
383 template <int dim, int spacedim>
384 void
386  const ArrayView<const Tensor<1, dim>> & input,
387  const MappingType mapping_type,
388  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
389  const ArrayView<Tensor<1, spacedim>> & output) const
390 {
391  AssertDimension(input.size(), output.size());
392 
393  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
394  Assert(data != nullptr, ExcInternalError());
395 
396  // check whether we should in fact work on the Q1 portion of it
397  if (data->use_mapping_q1_on_current_cell)
398  q1_mapping->transform(input, mapping_type, *data->mapping_q1_data, output);
399  else
400  // otherwise use the full mapping
401  qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output);
402 }
403 
404 
405 
406 template <int dim, int spacedim>
407 void
409  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
410  const MappingType mapping_type,
411  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
412  const ArrayView<Tensor<2, spacedim>> & output) const
413 {
414  AssertDimension(input.size(), output.size());
415  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
416  &mapping_data) != nullptr),
417  ExcInternalError());
418  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
419 
420  // check whether we should in fact work on the Q1 portion of it
422  q1_mapping->transform(input, mapping_type, *data->mapping_q1_data, output);
423  else
424  // otherwise use the full mapping
425  qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output);
426 }
427 
428 
429 
430 template <int dim, int spacedim>
431 void
433  const ArrayView<const Tensor<2, dim>> & input,
434  const MappingType mapping_type,
435  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
436  const ArrayView<Tensor<2, spacedim>> & output) const
437 {
438  AssertDimension(input.size(), output.size());
439  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
440  &mapping_data) != nullptr),
441  ExcInternalError());
442  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
443 
444  // check whether we should in fact work on the Q1 portion of it
446  q1_mapping->transform(input, mapping_type, *data->mapping_q1_data, output);
447  else
448  // otherwise use the full mapping
449  qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output);
450 }
451 
452 
453 
454 template <int dim, int spacedim>
455 void
457  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
458  const MappingType mapping_type,
459  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
460  const ArrayView<Tensor<3, spacedim>> & output) const
461 {
462  AssertDimension(input.size(), output.size());
463  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
464  &mapping_data) != nullptr),
465  ExcInternalError());
466  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
467 
468  // check whether we should in fact work on the Q1 portion of it
470  q1_mapping->transform(input, mapping_type, *data->mapping_q1_data, output);
471  else
472  // otherwise use the full mapping
473  qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output);
474 }
475 
476 
477 
478 template <int dim, int spacedim>
479 void
481  const ArrayView<const Tensor<3, dim>> & input,
482  const MappingType mapping_type,
483  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
484  const ArrayView<Tensor<3, spacedim>> & output) const
485 {
486  AssertDimension(input.size(), output.size());
487  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
488  &mapping_data) != nullptr),
489  ExcInternalError());
490  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
491 
492  // check whether we should in fact work on the Q1 portion of it
494  q1_mapping->transform(input, mapping_type, *data->mapping_q1_data, output);
495  else
496  // otherwise use the full mapping
497  qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output);
498 }
499 
500 
501 
502 template <int dim, int spacedim>
505  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
506  const Point<dim> & p) const
507 {
508  // first see, whether we want to use a linear or a higher order
509  // mapping, then either use our own facilities or that of the Q1
510  // mapping we store
511  if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
512  return qp_mapping->transform_unit_to_real_cell(cell, p);
513  else
514  return q1_mapping->transform_unit_to_real_cell(cell, p);
515 }
516 
517 
518 
519 template <int dim, int spacedim>
522  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
523  const Point<spacedim> & p) const
524 {
525  if (cell->has_boundary_lines() || use_mapping_q_on_all_cells ||
526  (dim != spacedim))
527  return qp_mapping->transform_real_to_unit_cell(cell, p);
528  else
529  return q1_mapping->transform_real_to_unit_cell(cell, p);
530 }
531 
532 
533 
534 template <int dim, int spacedim>
535 std::unique_ptr<Mapping<dim, spacedim>>
537 {
538  return std_cxx14::make_unique<MappingQ<dim, spacedim>>(
539  this->polynomial_degree, this->use_mapping_q_on_all_cells);
540 }
541 
542 
543 
544 // explicit instantiations
545 #include "mapping_q.inst"
546 
547 
548 DEAL_II_NAMESPACE_CLOSE
unsigned int get_degree() const
Definition: mapping_q.cc:123
virtual std::size_t memory_consumption() const override
Definition: mapping_q.cc:50
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::shared_ptr< const MappingQGeneric< dim, spacedim > > qp_mapping
Definition: mapping_q.h:368
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_q.cc:151
Task< RT > new_task(const std::function< RT()> &function)
MappingType
Definition: mapping.h:61
bool use_mapping_q1_on_current_cell
Definition: mapping_q.h:251
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_face_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition: mapping_q.cc:182
STL namespace.
const bool use_mapping_q_on_all_cells
Definition: mapping_q.h:333
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:738
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
Abstract base class for mapping classes.
Definition: dof_tools.h:57
const unsigned int polynomial_degree
Definition: mapping_q.h:327
std::unique_ptr< typename MappingQGeneric< dim, spacedim >::InternalData > mapping_qp_data
Definition: mapping_q.h:265
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_subface_data(const UpdateFlags flags, const Quadrature< dim - 1 > &quadrature) const override
Definition: mapping_q.cc:215
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const override
Definition: mapping_q.cc:504
virtual bool preserves_vertex_locations() const override
Definition: mapping_q.cc:132
virtual void transform(const ArrayView< const Tensor< 1, dim >> &input, const MappingType type, const typename Mapping< dim, spacedim >::InternalDataBase &internal, const ArrayView< Tensor< 1, spacedim >> &output) const override
Definition: mapping_q.cc:385
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_q.cc:536
std::unique_ptr< typename MappingQGeneric< dim, spacedim >::InternalData > mapping_q1_data
Definition: mapping_q.h:258
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_q.cc:141
MappingQ(const unsigned int polynomial_degree, const bool use_mapping_q_on_all_cells=false)
Definition: mapping_q.cc:62
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const override
Definition: mapping_q.cc:521
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()
std::shared_ptr< const MappingQGeneric< dim, spacedim > > q1_mapping
Definition: mapping_q.h:352