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