Reference documentation for deal.II version Git 53084d8c6a 2020-10-21 00:34:33 -0400
\(\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.cc
Go to the documentation of this file.
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 
22 #include <deal.II/base/utilities.h>
23 
25 
26 #include <deal.II/fe/fe_q.h>
27 #include <deal.II/fe/fe_values.h>
28 #include <deal.II/fe/mapping_q.h>
29 
31 
33 
34 #include <memory>
35 #include <numeric>
36 
38 
39 
40 template <int dim, int spacedim>
42  : use_mapping_q1_on_current_cell(false)
43 {}
44 
45 
46 
47 template <int dim, int spacedim>
48 std::size_t
50 {
51  return (
56 }
57 
58 
59 
60 template <int dim, int spacedim>
61 MappingQ<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 
88 template <int dim, int spacedim>
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  {
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 
120 template <int dim, int spacedim>
121 unsigned int
123 {
124  return polynomial_degree;
125 }
126 
127 
128 
129 template <int dim, int spacedim>
130 inline bool
132 {
133  return true;
134 }
135 
136 
137 
138 template <int dim, int spacedim>
141 {
142  return (q1_mapping->requires_update_flags(in) |
143  qp_mapping->requires_update_flags(in));
144 }
145 
146 
147 
148 template <int dim, int spacedim>
149 std::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<
172  typename MappingQGeneric<dim, spacedim>::InternalData>(
173  std::move(do_get_data.return_value()));
174  return data_ptr;
175 }
176 
177 
178 
179 template <int dim, int spacedim>
180 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
182  const UpdateFlags update_flags,
183  const Quadrature<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  // build the Q1 and Qp internal data objects in parallel
191  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
192  do_get_data =
194  *qp_mapping,
195  update_flags,
196  quadrature);
197 
199  data.mapping_q1_data = Utilities::dynamic_unique_cast<
201  std::move(q1_mapping->get_face_data(update_flags, quadrature)));
202 
203  // wait for the task above to finish and use returned value
204  data.mapping_qp_data = Utilities::dynamic_unique_cast<
205  typename MappingQGeneric<dim, spacedim>::InternalData>(
206  std::move(do_get_data.return_value()));
207  return data_ptr;
208 }
209 
210 
211 
212 template <int dim, int spacedim>
213 std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
215  const UpdateFlags update_flags,
216  const Quadrature<dim - 1> &quadrature) const
217 {
218  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
219  std::make_unique<InternalData>();
220  auto &data = dynamic_cast<InternalData &>(*data_ptr);
221 
222  // build the Q1 and Qp internal data objects in parallel
224  std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
225  do_get_data =
227  *qp_mapping,
228  update_flags,
229  quadrature);
230 
232  data.mapping_q1_data = Utilities::dynamic_unique_cast<
234  std::move(q1_mapping->get_subface_data(update_flags, quadrature)));
235 
236  // wait for the task above to finish and use returned value
237  data.mapping_qp_data = Utilities::dynamic_unique_cast<
238  typename MappingQGeneric<dim, spacedim>::InternalData>(
239  std::move(do_get_data.return_value()));
240  return data_ptr;
241 }
242 
243 
244 // Note that the CellSimilarity flag is modifiable, since MappingQ can need to
245 // recalculate data even when cells are similar.
246 template <int dim, int spacedim>
249  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
250  const CellSimilarity::Similarity cell_similarity,
251  const Quadrature<dim> & quadrature,
252  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
254  &output_data) const
255 {
256  // convert data object to internal data for this class. fails with an
257  // exception if that is not possible
258  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
259  ExcInternalError());
260  const InternalData &data = static_cast<const InternalData &>(internal_data);
261 
262  // check whether this cell needs the full mapping or can be treated by a
263  // reduced Q1 mapping, e.g. if the cell is in the interior of the domain
265  !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
266 
267 
268  // call the base class. we need to ensure that the flag indicating whether
269  // we can use some similarity has to be modified - for a general MappingQ,
270  // the data needs to be recomputed anyway since then the mapping changes the
271  // data. this needs to be known also for later operations, so modify the
272  // variable here. this also affects the calculation of the next cell -- if
273  // we use Q1 data on the next cell, the data will still be invalid.
274  const CellSimilarity::Similarity updated_cell_similarity =
275  ((data.use_mapping_q1_on_current_cell == false) &&
276  (this->polynomial_degree > 1) ?
278  cell_similarity);
279 
280  // depending on the results above, decide whether the Q1 mapping or
281  // the Qp mapping needs to handle this cell
283  q1_mapping->fill_fe_values(cell,
284  updated_cell_similarity,
285  quadrature,
286  *data.mapping_q1_data,
287  output_data);
288  else
289  qp_mapping->fill_fe_values(cell,
290  updated_cell_similarity,
291  quadrature,
292  *data.mapping_qp_data,
293  output_data);
294 
295  return updated_cell_similarity;
296 }
297 
298 
299 
300 template <int dim, int spacedim>
301 void
303  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
304  const unsigned int face_no,
305  const Quadrature<dim - 1> & quadrature,
306  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
308  &output_data) const
309 {
310  // convert data object to internal data for this class. fails with an
311  // exception if that is not possible
312  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
313  ExcInternalError());
314  const InternalData &data = static_cast<const InternalData &>(internal_data);
315 
316  // check whether this cell needs the full mapping or can be treated by a
317  // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
318  // domain. note that it is not sufficient to ask whether the present _face_
319  // is in the interior, as the mapping on the face depends on the mapping of
320  // the cell, which in turn depends on the fact whether _any_ of the faces of
321  // this cell is at the boundary, not only the present face
323  !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
324 
325  // depending on the results above, decide whether the Q1 mapping or
326  // the Qp mapping needs to handle this cell
328  q1_mapping->fill_fe_face_values(
329  cell, face_no, quadrature, *data.mapping_q1_data, output_data);
330  else
331  qp_mapping->fill_fe_face_values(
332  cell, face_no, quadrature, *data.mapping_qp_data, output_data);
333 }
334 
335 
336 template <int dim, int spacedim>
337 void
339  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
340  const unsigned int face_no,
341  const unsigned int subface_no,
342  const Quadrature<dim - 1> & quadrature,
343  const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
345  &output_data) const
346 {
347  // convert data object to internal data for this class. fails with an
348  // exception if that is not possible
349  Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
350  ExcInternalError());
351  const InternalData &data = static_cast<const InternalData &>(internal_data);
352 
353  // check whether this cell needs the full mapping or can be treated by a
354  // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
355  // domain. note that it is not sufficient to ask whether the present _face_
356  // is in the interior, as the mapping on the face depends on the mapping of
357  // the cell, which in turn depends on the fact whether _any_ of the faces of
358  // this cell is at the boundary, not only the present face
360  !(use_mapping_q_on_all_cells || cell->has_boundary_lines());
361 
362  // depending on the results above, decide whether the Q1 mapping or
363  // the Qp mapping needs to handle this cell
365  q1_mapping->fill_fe_subface_values(cell,
366  face_no,
367  subface_no,
368  quadrature,
369  *data.mapping_q1_data,
370  output_data);
371  else
372  qp_mapping->fill_fe_subface_values(cell,
373  face_no,
374  subface_no,
375  quadrature,
376  *data.mapping_qp_data,
377  output_data);
378 }
379 
380 
381 
382 template <int dim, int spacedim>
383 void
385  const ArrayView<const Tensor<1, dim>> & input,
386  const MappingKind mapping_kind,
387  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
388  const ArrayView<Tensor<1, spacedim>> & output) const
389 {
390  AssertDimension(input.size(), output.size());
391 
392  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
393  Assert(data != nullptr, ExcInternalError());
394 
395  // check whether we should in fact work on the Q1 portion of it
396  if (data->use_mapping_q1_on_current_cell)
397  q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
398  else
399  // otherwise use the full mapping
400  qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
401 }
402 
403 
404 
405 template <int dim, int spacedim>
406 void
408  const ArrayView<const DerivativeForm<1, dim, spacedim>> &input,
409  const MappingKind mapping_kind,
410  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
411  const ArrayView<Tensor<2, spacedim>> & output) const
412 {
413  AssertDimension(input.size(), output.size());
414  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
415  &mapping_data) != nullptr),
416  ExcInternalError());
417  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
418 
419  // check whether we should in fact work on the Q1 portion of it
421  q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
422  else
423  // otherwise use the full mapping
424  qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
425 }
426 
427 
428 
429 template <int dim, int spacedim>
430 void
432  const ArrayView<const Tensor<2, dim>> & input,
433  const MappingKind mapping_kind,
434  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
435  const ArrayView<Tensor<2, spacedim>> & output) const
436 {
437  AssertDimension(input.size(), output.size());
438  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
439  &mapping_data) != nullptr),
440  ExcInternalError());
441  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
442 
443  // check whether we should in fact work on the Q1 portion of it
445  q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
446  else
447  // otherwise use the full mapping
448  qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
449 }
450 
451 
452 
453 template <int dim, int spacedim>
454 void
456  const ArrayView<const DerivativeForm<2, dim, spacedim>> &input,
457  const MappingKind mapping_kind,
458  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
459  const ArrayView<Tensor<3, spacedim>> & output) const
460 {
461  AssertDimension(input.size(), output.size());
462  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
463  &mapping_data) != nullptr),
464  ExcInternalError());
465  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
466 
467  // check whether we should in fact work on the Q1 portion of it
469  q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
470  else
471  // otherwise use the full mapping
472  qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
473 }
474 
475 
476 
477 template <int dim, int spacedim>
478 void
480  const ArrayView<const Tensor<3, dim>> & input,
481  const MappingKind mapping_kind,
482  const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
483  const ArrayView<Tensor<3, spacedim>> & output) const
484 {
485  AssertDimension(input.size(), output.size());
486  Assert((dynamic_cast<const typename MappingQ<dim, spacedim>::InternalData *>(
487  &mapping_data) != nullptr),
488  ExcInternalError());
489  const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data);
490 
491  // check whether we should in fact work on the Q1 portion of it
493  q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output);
494  else
495  // otherwise use the full mapping
496  qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output);
497 }
498 
499 
500 
501 template <int dim, int spacedim>
504  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
505  const Point<dim> & p) const
506 {
507  // first see, whether we want to use a linear or a higher order
508  // mapping, then either use our own facilities or that of the Q1
509  // mapping we store
510  if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
511  return qp_mapping->transform_unit_to_real_cell(cell, p);
512  else
513  return q1_mapping->transform_unit_to_real_cell(cell, p);
514 }
515 
516 
517 
518 template <int dim, int spacedim>
521  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
522  const Point<spacedim> & p) const
523 {
524  if (cell->has_boundary_lines() || use_mapping_q_on_all_cells ||
525  (dim != spacedim))
526  return qp_mapping->transform_real_to_unit_cell(cell, p);
527  else
528  return q1_mapping->transform_real_to_unit_cell(cell, p);
529 }
530 
531 
532 
533 template <int dim, int spacedim>
534 std::unique_ptr<Mapping<dim, spacedim>>
536 {
537  return std::make_unique<MappingQ<dim, spacedim>>(
539 }
540 
541 
542 
543 // explicit instantiations
544 #include "mapping_q.inst"
545 
546 
unsigned int get_degree() const
Definition: mapping_q.cc:122
virtual std::size_t memory_consumption() const override
Definition: mapping_q.cc:49
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
std::shared_ptr< const MappingQGeneric< dim, spacedim > > qp_mapping
Definition: mapping_q.h:365
virtual std::unique_ptr< typename Mapping< dim, spacedim >::InternalDataBase > get_data(const UpdateFlags, const Quadrature< dim > &quadrature) const override
Definition: mapping_q.cc:150
Task< RT > new_task(const std::function< RT()> &function)
bool use_mapping_q1_on_current_cell
Definition: mapping_q.h:248
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:181
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
Definition: mapping_q.cc:384
STL namespace.
const bool use_mapping_q_on_all_cells
Definition: mapping_q.h:330
MappingKind
Definition: mapping.h:62
std::unique_ptr< To > dynamic_unique_cast(std::unique_ptr< From > &&p)
Definition: utilities.h:752
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
Definition: mapping_q.cc:338
#define Assert(cond, exc)
Definition: exceptions.h:1423
UpdateFlags
Abstract base class for mapping classes.
Definition: mapping.h:301
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
const unsigned int polynomial_degree
Definition: mapping_q.h:324
std::unique_ptr< typename MappingQGeneric< dim, spacedim >::InternalData > mapping_qp_data
Definition: mapping_q.h:262
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:214
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:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
virtual bool preserves_vertex_locations() const override
Definition: mapping_q.cc:131
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_q.cc:535
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
Definition: mapping_q.cc:248
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const Quadrature< dim - 1 > &quadrature, const typename Mapping< dim, spacedim >::InternalDataBase &internal_data, internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &output_data) const override
Definition: mapping_q.cc:302
std::unique_ptr< typename MappingQGeneric< dim, spacedim >::InternalData > mapping_q1_data
Definition: mapping_q.h:255
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
Definition: mapping_q.cc:140
MappingQ(const unsigned int polynomial_degree, const bool use_mapping_q_on_all_cells=false)
Definition: mapping_q.cc:61
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:520
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:349