Reference documentation for deal.II version GIT f6a5d312c9 2023-10-04 08:50:02+00:00
\(\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\}}\)
fe_values.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 #include <deal.II/base/numbers.h>
18 
28 #include <deal.II/lac/vector.h>
29 
31 
33 
34 namespace NonMatching
35 {
39  , surface(update_default)
40  {}
41 
42 
43 
44  template <int dim>
45  template <typename VectorType>
47  const Quadrature<1> &quadrature,
48  const RegionUpdateFlags region_update_flags,
49  const MeshClassifier<dim> &mesh_classifier,
50  const DoFHandler<dim> &dof_handler,
51  const VectorType &level_set,
52  const AdditionalData &additional_data)
53  : mapping_collection(&::hp::StaticMappingQ1<dim>::mapping_collection)
54  , fe_collection(&fe_collection)
55  , q_collection_1D(quadrature)
56  , region_update_flags(region_update_flags)
57  , mesh_classifier(&mesh_classifier)
58  , quadrature_generator(q_collection_1D,
59  dof_handler,
60  level_set,
61  additional_data)
62  {
63  // Tensor products of each quadrature in q_collection_1d. Used on the
64  // non-intersected cells.
65  hp::QCollection<dim> q_collection;
66  for (const auto &quadrature : q_collection_1D)
67  q_collection.push_back(Quadrature<dim>(quadrature));
68 
69  initialize(q_collection);
70  }
71 
72 
73 
74  template <int dim>
75  template <typename VectorType>
77  const hp::FECollection<dim> &fe_collection,
78  const hp::QCollection<dim> &q_collection,
79  const hp::QCollection<1> &q_collection_1D,
80  const RegionUpdateFlags region_update_flags,
81  const MeshClassifier<dim> &mesh_classifier,
82  const DoFHandler<dim> &dof_handler,
83  const VectorType &level_set,
84  const AdditionalData &additional_data)
85  : mapping_collection(&mapping_collection)
86  , fe_collection(&fe_collection)
87  , q_collection_1D(q_collection_1D)
88  , region_update_flags(region_update_flags)
89  , mesh_classifier(&mesh_classifier)
90  , quadrature_generator(q_collection_1D,
91  dof_handler,
92  level_set,
93  additional_data)
94  {
95  initialize(q_collection);
96  }
97 
98 
99 
100  template <int dim>
101  void
103  {
104  current_cell_location = LocationToLevelSet::unassigned;
105  active_fe_index = numbers::invalid_unsigned_int;
106 
107  Assert(fe_collection->size() > 0,
108  ExcMessage("Incoming hp::FECollection can not be empty."));
109  Assert(mapping_collection->size() == fe_collection->size() ||
110  mapping_collection->size() == 1,
111  ExcMessage("Size of hp::MappingCollection must be "
112  "the same as hp::FECollection or 1."));
113  Assert(q_collection.size() == fe_collection->size() ||
114  q_collection.size() == 1,
115  ExcMessage("Size of hp::QCollection<dim> must be the "
116  "same as hp::FECollection or 1."));
117  Assert(q_collection_1D.size() == fe_collection->size() ||
118  q_collection_1D.size() == 1,
119  ExcMessage("Size of hp::QCollection<1> must be the "
120  "same as hp::FECollection or 1."));
121 
122  // For each element in fe_collection, create ::FEValues objects to use
123  // on the non-intersected cells.
124  fe_values_inside_full_quadrature.resize(fe_collection->size());
125  fe_values_outside_full_quadrature.resize(fe_collection->size());
126  for (unsigned int fe_index = 0; fe_index < fe_collection->size();
127  ++fe_index)
128  {
129  const unsigned int mapping_index =
130  mapping_collection->size() > 1 ? fe_index : 0;
131  const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
132 
133  fe_values_inside_full_quadrature[fe_index].emplace(
134  (*mapping_collection)[mapping_index],
135  (*fe_collection)[fe_index],
136  q_collection[q_index],
137  region_update_flags.inside);
138  fe_values_outside_full_quadrature[fe_index].emplace(
139  (*mapping_collection)[mapping_index],
140  (*fe_collection)[fe_index],
141  q_collection[q_index],
142  region_update_flags.outside);
143  }
144  }
145 
146 
147 
148  template <int dim>
149  template <bool level_dof_access>
150  void
153  {
154  current_cell_location = mesh_classifier->location_to_level_set(cell);
155  active_fe_index = cell->active_fe_index();
156 
157  // These objects were created with a quadrature based on the previous cell
158  // and are thus no longer valid.
159  fe_values_inside.reset();
160  fe_values_surface.reset();
161  fe_values_outside.reset();
162 
163  switch (current_cell_location)
164  {
166  {
167  fe_values_inside_full_quadrature.at(active_fe_index)->reinit(cell);
168  break;
169  }
171  {
172  fe_values_outside_full_quadrature.at(active_fe_index)->reinit(cell);
173  break;
174  }
176  {
177  const unsigned int mapping_index =
178  mapping_collection->size() > 1 ? active_fe_index : 0;
179 
180  const unsigned int q1D_index =
181  q_collection_1D.size() > 1 ? active_fe_index : 0;
182  quadrature_generator.set_1D_quadrature(q1D_index);
183  quadrature_generator.generate(cell);
184 
185  const Quadrature<dim> &inside_quadrature =
186  quadrature_generator.get_inside_quadrature();
187  const Quadrature<dim> &outside_quadrature =
188  quadrature_generator.get_outside_quadrature();
189  const ImmersedSurfaceQuadrature<dim> &surface_quadrature =
190  quadrature_generator.get_surface_quadrature();
191 
192  // Even if a cell is formally intersected the number of created
193  // quadrature points can be 0. Avoid creating an FEValues object
194  // if that is the case.
195  if (inside_quadrature.size() > 0)
196  {
197  fe_values_inside.emplace((*mapping_collection)[mapping_index],
198  (*fe_collection)[active_fe_index],
199  inside_quadrature,
200  region_update_flags.inside);
201 
202  fe_values_inside->reinit(cell);
203  }
204 
205  if (outside_quadrature.size() > 0)
206  {
207  fe_values_outside.emplace((*mapping_collection)[mapping_index],
208  (*fe_collection)[active_fe_index],
209  outside_quadrature,
210  region_update_flags.outside);
211 
212  fe_values_outside->reinit(cell);
213  }
214 
215  if (surface_quadrature.size() > 0)
216  {
217  fe_values_surface.emplace((*mapping_collection)[mapping_index],
218  (*fe_collection)[active_fe_index],
219  surface_quadrature,
220  region_update_flags.surface);
221  fe_values_surface->reinit(cell);
222  }
223 
224  break;
225  }
226  default:
227  {
228  Assert(false, ExcInternalError());
229  break;
230  }
231  }
232  }
233 
234 
235 
236  template <int dim>
237  const std::optional<::FEValues<dim>> &
239  {
240  if (current_cell_location == LocationToLevelSet::inside)
241  return fe_values_inside_full_quadrature.at(active_fe_index);
242  else
243  return fe_values_inside;
244  }
245 
246 
247 
248  template <int dim>
249  const std::optional<::FEValues<dim>> &
251  {
252  if (current_cell_location == LocationToLevelSet::outside)
253  return fe_values_outside_full_quadrature.at(active_fe_index);
254  else
255  return fe_values_outside;
256  }
257 
258 
259 
260  template <int dim>
261  const std::optional<FEImmersedSurfaceValues<dim>> &
263  {
264  return fe_values_surface;
265  }
266 
267 
268 
269  template <int dim>
270  template <typename VectorType>
272  const hp::FECollection<dim> &fe_collection,
273  const Quadrature<1> &quadrature,
274  const RegionUpdateFlags region_update_flags,
275  const MeshClassifier<dim> &mesh_classifier,
276  const DoFHandler<dim> &dof_handler,
277  const VectorType &level_set,
278  const AdditionalData &additional_data)
279  : mapping_collection(&::hp::StaticMappingQ1<dim>::mapping_collection)
280  , fe_collection(&fe_collection)
281  , q_collection_1D(quadrature)
282  , region_update_flags(region_update_flags)
283  , mesh_classifier(&mesh_classifier)
284  , face_quadrature_generator(q_collection_1D,
285  dof_handler,
286  level_set,
287  additional_data)
288  {
289  // Tensor products of each quadrature in q_collection_1d. Used on the
290  // non-intersected cells.
291  hp::QCollection<dim - 1> q_collection;
292  for (const auto &quadrature : q_collection_1D)
293  q_collection.push_back(quadrature);
294 
295  initialize(q_collection);
296  }
297 
298 
299 
300  template <int dim>
301  template <typename VectorType>
303  const hp::MappingCollection<dim> &mapping_collection,
304  const hp::FECollection<dim> &fe_collection,
305  const hp::QCollection<dim - 1> &q_collection,
306  const hp::QCollection<1> &q_collection_1D,
307  const RegionUpdateFlags region_update_flags,
308  const MeshClassifier<dim> &mesh_classifier,
309  const DoFHandler<dim> &dof_handler,
310  const VectorType &level_set,
311  const AdditionalData &additional_data)
312  : mapping_collection(&mapping_collection)
313  , fe_collection(&fe_collection)
314  , q_collection_1D(q_collection_1D)
315  , region_update_flags(region_update_flags)
316  , mesh_classifier(&mesh_classifier)
317  , face_quadrature_generator(q_collection_1D,
318  dof_handler,
319  level_set,
320  additional_data)
321  {
322  initialize(q_collection);
323  }
324 
325 
326 
327  template <int dim>
328  void
330  const hp::QCollection<dim - 1> &q_collection)
331  {
332  current_face_location = LocationToLevelSet::unassigned;
333  active_fe_index = numbers::invalid_unsigned_int;
334 
335  Assert(fe_collection->size() > 0,
336  ExcMessage("Incoming hp::FECollection can not be empty."));
337  Assert(
338  mapping_collection->size() == fe_collection->size() ||
339  mapping_collection->size() == 1,
340  ExcMessage(
341  "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
342  Assert(
343  q_collection.size() == fe_collection->size() || q_collection.size() == 1,
344  ExcMessage(
345  "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
346  Assert(
347  q_collection_1D.size() == fe_collection->size() ||
348  q_collection_1D.size() == 1,
349  ExcMessage(
350  "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
351 
352  // For each element in fe_collection, create ::FEInterfaceValues
353  // objects to use on the non-intersected cells.
354  fe_values_inside_full_quadrature.resize(fe_collection->size());
355  fe_values_outside_full_quadrature.resize(fe_collection->size());
356  for (unsigned int fe_index = 0; fe_index < fe_collection->size();
357  ++fe_index)
358  {
359  const unsigned int mapping_index =
360  mapping_collection->size() > 1 ? fe_index : 0;
361  const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
362 
363  fe_values_inside_full_quadrature[fe_index].emplace(
364  (*mapping_collection)[mapping_index],
365  (*fe_collection)[fe_index],
366  q_collection[q_index],
367  region_update_flags.inside);
368  fe_values_outside_full_quadrature[fe_index].emplace(
369  (*mapping_collection)[mapping_index],
370  (*fe_collection)[fe_index],
371  q_collection[q_index],
372  region_update_flags.outside);
373  }
374  }
375 
376 
377 
378  template <int dim>
379  template <bool level_dof_access>
380  void
383  const unsigned int face_no,
384  const std::function<void(::FEInterfaceValues<dim> &)> &call_reinit)
385  {
386  current_face_location =
387  mesh_classifier->location_to_level_set(cell, face_no);
388  active_fe_index = cell->active_fe_index();
389 
390  // These objects were created with a quadrature based on the previous cell
391  // and are thus no longer valid.
392  fe_values_inside.reset();
393  fe_values_outside.reset();
394 
395  switch (current_face_location)
396  {
398  {
399  call_reinit(*fe_values_inside_full_quadrature.at(active_fe_index));
400  break;
401  }
403  {
404  call_reinit(*fe_values_outside_full_quadrature.at(active_fe_index));
405  break;
406  }
408  {
409  const unsigned int mapping_index =
410  mapping_collection->size() > 1 ? active_fe_index : 0;
411  const unsigned int q1D_index =
412  q_collection_1D.size() > 1 ? active_fe_index : 0;
413 
414  face_quadrature_generator.set_1D_quadrature(q1D_index);
415  face_quadrature_generator.generate(cell, face_no);
416 
417  const Quadrature<dim - 1> &inside_quadrature =
418  face_quadrature_generator.get_inside_quadrature();
419  const Quadrature<dim - 1> &outside_quadrature =
420  face_quadrature_generator.get_outside_quadrature();
421 
422  // Even if a cell is formally intersected the number of created
423  // quadrature points can be 0. Avoid creating an FEInterfaceValues
424  // object if that is the case.
425  if (inside_quadrature.size() > 0)
426  {
427  fe_values_inside.emplace((*mapping_collection)[mapping_index],
428  (*fe_collection)[active_fe_index],
429  inside_quadrature,
430  region_update_flags.inside);
431 
432  call_reinit(*fe_values_inside);
433  }
434 
435  if (outside_quadrature.size() > 0)
436  {
437  fe_values_outside.emplace((*mapping_collection)[mapping_index],
438  (*fe_collection)[active_fe_index],
439  outside_quadrature,
440  region_update_flags.outside);
441 
442  call_reinit(*fe_values_outside);
443  }
444  break;
445  }
446  default:
447  {
448  Assert(false, ExcInternalError());
449  break;
450  }
451  }
452  }
453 
454 
455 
456  template <int dim>
457  const std::optional<::FEInterfaceValues<dim>> &
459  {
460  if (current_face_location == LocationToLevelSet::inside)
461  return fe_values_inside_full_quadrature.at(active_fe_index);
462  else
463  return fe_values_inside;
464  }
465 
466 
467 
468  template <int dim>
469  const std::optional<::FEInterfaceValues<dim>> &
471  {
472  if (current_face_location == LocationToLevelSet::outside)
473  return fe_values_outside_full_quadrature.at(active_fe_index);
474  else
475  return fe_values_outside;
476  }
477 
478 
479 #include "fe_values.inst"
480 
481 } // namespace NonMatching
const std::optional<::FEInterfaceValues< dim > > & get_outside_fe_values() const
Definition: fe_values.cc:470
void initialize(const hp::QCollection< dim - 1 > &q_collection)
Definition: fe_values.cc:329
const hp::QCollection< 1 > q_collection_1D
Definition: fe_values.h:603
typename FaceQuadratureGenerator< dim >::AdditionalData AdditionalData
Definition: fe_values.h:440
void do_reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access >> &cell, const unsigned int face_no, const std::function< void(::FEInterfaceValues< dim > &)> &call_reinit)
Definition: fe_values.cc:381
const std::optional<::FEInterfaceValues< dim > > & get_inside_fe_values() const
Definition: fe_values.cc:458
FEInterfaceValues(const hp::FECollection< dim > &fe_collection, const Quadrature< 1 > &quadrature, const RegionUpdateFlags region_update_flags, const MeshClassifier< dim > &mesh_classifier, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
Definition: fe_values.cc:271
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access >> &cell)
Definition: fe_values.cc:151
const std::optional<::FEValues< dim > > & get_outside_fe_values() const
Definition: fe_values.cc:250
void initialize(const hp::QCollection< dim > &q_collection)
Definition: fe_values.cc:102
const std::optional< FEImmersedSurfaceValues< dim > > & get_surface_fe_values() const
Definition: fe_values.cc:262
const hp::QCollection< 1 > q_collection_1D
Definition: fe_values.h:284
typename QuadratureGenerator< dim >::AdditionalData AdditionalData
Definition: fe_values.h:147
FEValues(const hp::FECollection< dim > &fe_collection, const Quadrature< 1 > &quadrature, const RegionUpdateFlags region_update_flags, const MeshClassifier< dim > &mesh_classifier, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
Definition: fe_values.cc:46
const std::optional<::FEValues< dim > > & get_inside_fe_values() const
Definition: fe_values.cc:238
unsigned int size() const
unsigned int size() const
Definition: collection.h:265
void push_back(const Quadrature< dim_in > &new_quadrature)
Definition: q_collection.h:216
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcMessage(std::string arg1)
@ update_default
No update.
Definition: hp.h:118
static const unsigned int invalid_unsigned_int
Definition: types.h:213
unsigned short int fe_index
Definition: types.h:60