Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
fe_values.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
27#include <deal.II/lac/vector.h>
28
30
32
33namespace NonMatching
34{
40
41
42
43 template <int dim>
44 template <typename VectorType>
46 const Quadrature<1> &quadrature,
47 const RegionUpdateFlags region_update_flags,
48 const MeshClassifier<dim> &mesh_classifier,
49 const DoFHandler<dim> &dof_handler,
50 const VectorType &level_set,
51 const AdditionalData &additional_data)
52 : mapping_collection(&::hp::StaticMappingQ1<dim>::mapping_collection)
53 , fe_collection(&fe_collection)
54 , q_collection_1D(quadrature)
55 , region_update_flags(region_update_flags)
56 , mesh_classifier(&mesh_classifier)
57 , quadrature_generator(q_collection_1D,
58 dof_handler,
59 level_set,
60 additional_data)
61 {
62 // Tensor products of each quadrature in q_collection_1d. Used on the
63 // non-intersected cells.
64 hp::QCollection<dim> q_collection;
65 for (const auto &quadrature : q_collection_1D)
66 q_collection.push_back(Quadrature<dim>(quadrature));
67
68 initialize(q_collection);
69 }
70
71
72
73 template <int dim>
74 template <typename VectorType>
76 const hp::FECollection<dim> &fe_collection,
77 const hp::QCollection<dim> &q_collection,
78 const hp::QCollection<1> &q_collection_1D,
79 const RegionUpdateFlags region_update_flags,
80 const MeshClassifier<dim> &mesh_classifier,
81 const DoFHandler<dim> &dof_handler,
82 const VectorType &level_set,
83 const AdditionalData &additional_data)
84 : mapping_collection(&mapping_collection)
85 , fe_collection(&fe_collection)
86 , q_collection_1D(q_collection_1D)
87 , region_update_flags(region_update_flags)
88 , mesh_classifier(&mesh_classifier)
89 , quadrature_generator(q_collection_1D,
90 dof_handler,
91 level_set,
92 additional_data)
93 {
94 initialize(q_collection);
95 }
96
97
98
99 template <int dim>
100 void
102 {
103 current_cell_location = LocationToLevelSet::unassigned;
104 active_fe_index = numbers::invalid_unsigned_int;
105
106 Assert(fe_collection->size() > 0,
107 ExcMessage("Incoming hp::FECollection can not be empty."));
108 Assert(mapping_collection->size() == fe_collection->size() ||
109 mapping_collection->size() == 1,
110 ExcMessage("Size of hp::MappingCollection must be "
111 "the same as hp::FECollection or 1."));
112 Assert(q_collection.size() == fe_collection->size() ||
113 q_collection.size() == 1,
114 ExcMessage("Size of hp::QCollection<dim> must be the "
115 "same as hp::FECollection or 1."));
116 Assert(q_collection_1D.size() == fe_collection->size() ||
117 q_collection_1D.size() == 1,
118 ExcMessage("Size of hp::QCollection<1> must be the "
119 "same as hp::FECollection or 1."));
120
121 // For each element in fe_collection, create ::FEValues objects to use
122 // on the non-intersected cells.
123 fe_values_inside_full_quadrature.resize(fe_collection->size());
124 fe_values_outside_full_quadrature.resize(fe_collection->size());
125 for (unsigned int fe_index = 0; fe_index < fe_collection->size();
126 ++fe_index)
127 {
128 const unsigned int mapping_index =
129 mapping_collection->size() > 1 ? fe_index : 0;
130 const unsigned int q_index = q_collection.size() > 1 ? fe_index : 0;
131
132 fe_values_inside_full_quadrature[fe_index].emplace(
133 (*mapping_collection)[mapping_index],
134 (*fe_collection)[fe_index],
135 q_collection[q_index],
136 region_update_flags.inside);
137 fe_values_outside_full_quadrature[fe_index].emplace(
138 (*mapping_collection)[mapping_index],
139 (*fe_collection)[fe_index],
140 q_collection[q_index],
141 region_update_flags.outside);
142 }
143 }
144
145
146
147 template <int dim>
148 template <bool level_dof_access>
149 void
152 const unsigned int q_index,
153 const unsigned int mapping_index)
154 {
155 this->reinit_internal(cell,
156 q_index,
157 mapping_index,
158 cell->active_fe_index());
159 }
160
161
162
163 template <int dim>
164 void
166 const unsigned int q_index,
167 const unsigned int mapping_index,
168 const unsigned int fe_index)
169 {
170 this->reinit_internal(cell, q_index, mapping_index, fe_index);
171 }
172
173
174
175 template <int dim>
176 template <typename CellIteratorType>
177 void
178 FEValues<dim>::reinit_internal(const CellIteratorType &cell,
179 const unsigned int q_index_in,
180 const unsigned int mapping_index_in,
181 const unsigned int fe_index_in)
182 {
183 current_cell_location = mesh_classifier->location_to_level_set(cell);
184
185 if (fe_index_in == numbers::invalid_unsigned_int)
186 this->active_fe_index = 0;
187 else
188 this->active_fe_index = fe_index_in;
189
190 unsigned int mapping_index = mapping_index_in;
191 unsigned int q_index = q_index_in;
192 unsigned int q_index_1D = q_index_in;
193
194 if (mapping_index == numbers::invalid_unsigned_int)
195 {
196 if (mapping_collection->size() > 1)
197 mapping_index = active_fe_index;
198 else
199 mapping_index = 0;
200 }
201
202 if (q_index == numbers::invalid_unsigned_int)
203 {
204 if (fe_values_inside_full_quadrature.size() > 1)
205 q_index = active_fe_index;
206 else
207 q_index = 0;
208 }
209
210 if (q_index_1D == numbers::invalid_unsigned_int)
211 {
212 if (q_collection_1D.size() > 1)
213 q_index_1D = active_fe_index;
214 else
215 q_index_1D = 0;
216 }
217
218 // These objects were created with a quadrature based on the previous cell
219 // and are thus no longer valid.
220 fe_values_inside.reset();
221 fe_values_surface.reset();
222 fe_values_outside.reset();
223
224 switch (current_cell_location)
225 {
227 {
228 Assert((active_fe_index == mapping_index) ||
229 ((mapping_collection->size() == 1) &&
230 (mapping_index == 0)),
232 Assert(active_fe_index == q_index, ExcNotImplemented());
233
234 fe_values_inside_full_quadrature.at(q_index)->reinit(cell);
235 break;
236 }
238 {
239 Assert((active_fe_index == mapping_index) ||
240 ((mapping_collection->size() == 1) &&
241 (mapping_index == 0)),
243 Assert(active_fe_index == q_index, ExcNotImplemented());
244
245 fe_values_outside_full_quadrature.at(q_index)->reinit(cell);
246 break;
247 }
249 {
250 quadrature_generator.set_1D_quadrature(q_index_1D);
251 quadrature_generator.generate(cell);
252
253 const Quadrature<dim> &inside_quadrature =
254 quadrature_generator.get_inside_quadrature();
255 const Quadrature<dim> &outside_quadrature =
256 quadrature_generator.get_outside_quadrature();
257 const ImmersedSurfaceQuadrature<dim> &surface_quadrature =
258 quadrature_generator.get_surface_quadrature();
259
260 // Even if a cell is formally intersected the number of created
261 // quadrature points can be 0. Avoid creating an FEValues object
262 // if that is the case.
263 if (inside_quadrature.size() > 0)
264 {
265 fe_values_inside.emplace((*mapping_collection)[mapping_index],
266 (*fe_collection)[active_fe_index],
267 inside_quadrature,
268 region_update_flags.inside);
269
270 fe_values_inside->reinit(cell);
271 }
272
273 if (outside_quadrature.size() > 0)
274 {
275 fe_values_outside.emplace((*mapping_collection)[mapping_index],
276 (*fe_collection)[active_fe_index],
277 outside_quadrature,
278 region_update_flags.outside);
279
280 fe_values_outside->reinit(cell);
281 }
282
283 if (surface_quadrature.size() > 0)
284 {
285 fe_values_surface.emplace((*mapping_collection)[mapping_index],
286 (*fe_collection)[active_fe_index],
287 surface_quadrature,
288 region_update_flags.surface);
289 fe_values_surface->reinit(cell);
290 }
291
292 break;
293 }
294 default:
295 {
297 break;
298 }
299 }
300 }
301
302
303
304 template <int dim>
305 const std::optional<::FEValues<dim>> &
307 {
308 if (current_cell_location == LocationToLevelSet::inside)
309 return fe_values_inside_full_quadrature.at(active_fe_index);
310 else
311 return fe_values_inside;
312 }
313
314
315
316 template <int dim>
317 const std::optional<::FEValues<dim>> &
319 {
320 if (current_cell_location == LocationToLevelSet::outside)
321 return fe_values_outside_full_quadrature.at(active_fe_index);
322 else
323 return fe_values_outside;
324 }
325
326
327
328 template <int dim>
329 const std::optional<FEImmersedSurfaceValues<dim>> &
331 {
332 return fe_values_surface;
333 }
334
335
336
337 template <int dim>
338 template <typename VectorType>
340 const hp::FECollection<dim> &fe_collection,
341 const Quadrature<1> &quadrature,
342 const RegionUpdateFlags region_update_flags,
343 const MeshClassifier<dim> &mesh_classifier,
344 const DoFHandler<dim> &dof_handler,
345 const VectorType &level_set,
346 const AdditionalData &additional_data)
347 : mapping_collection(&::hp::StaticMappingQ1<dim>::mapping_collection)
348 , fe_collection(&fe_collection)
349 , q_collection_1D(quadrature)
350 , region_update_flags(region_update_flags)
351 , mesh_classifier(&mesh_classifier)
352 , face_quadrature_generator(q_collection_1D,
353 dof_handler,
354 level_set,
355 additional_data)
356 {
357 // Tensor products of each quadrature in q_collection_1d. Used on the
358 // non-intersected cells.
359 hp::QCollection<dim - 1> q_collection;
360 for (const auto &quadrature : q_collection_1D)
361 q_collection.push_back(quadrature);
362
363 initialize(q_collection);
364 }
365
366
367
368 template <int dim>
369 template <typename VectorType>
371 const hp::MappingCollection<dim> &mapping_collection,
372 const hp::FECollection<dim> &fe_collection,
373 const hp::QCollection<dim - 1> &q_collection,
374 const hp::QCollection<1> &q_collection_1D,
375 const RegionUpdateFlags region_update_flags,
376 const MeshClassifier<dim> &mesh_classifier,
377 const DoFHandler<dim> &dof_handler,
378 const VectorType &level_set,
379 const AdditionalData &additional_data)
380 : mapping_collection(&mapping_collection)
381 , fe_collection(&fe_collection)
382 , q_collection_1D(q_collection_1D)
383 , region_update_flags(region_update_flags)
384 , mesh_classifier(&mesh_classifier)
385 , face_quadrature_generator(q_collection_1D,
386 dof_handler,
387 level_set,
388 additional_data)
389 {
390 initialize(q_collection);
391 }
392
393
394
395 template <int dim>
396 void
398 const hp::QCollection<dim - 1> &q_collection)
399 {
400 current_face_location = LocationToLevelSet::unassigned;
401
402 Assert(fe_collection->size() > 0,
403 ExcMessage("Incoming hp::FECollection can not be empty."));
404 Assert(
405 mapping_collection->size() == fe_collection->size() ||
406 mapping_collection->size() == 1,
408 "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
409 Assert(
410 q_collection.size() == fe_collection->size() || q_collection.size() == 1,
412 "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
413 Assert(
414 q_collection_1D.size() == fe_collection->size() ||
415 q_collection_1D.size() == 1,
417 "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
418
419 fe_values_inside_full_quadrature.emplace(*mapping_collection,
420 *fe_collection,
421 q_collection,
422 region_update_flags.inside);
423 fe_values_outside_full_quadrature.emplace(*mapping_collection,
424 *fe_collection,
425 q_collection,
426 region_update_flags.outside);
427 }
428
429
430
431 template <int dim>
432 template <typename CellAccessorType>
433 void
436 const unsigned int face_no,
437 const unsigned int q_index_in,
438 const unsigned int active_fe_index_in,
439 const std::function<void(::FEInterfaceValues<dim> &,
440 const unsigned int)> &call_reinit)
441 {
442 current_face_location =
443 mesh_classifier->location_to_level_set(cell, face_no);
444
445 // These objects were created with a quadrature based on the previous cell
446 // and are thus no longer valid.
447 fe_values_inside.reset();
448 fe_values_outside.reset();
449
450 switch (current_face_location)
451 {
453 {
454 call_reinit(*fe_values_inside_full_quadrature, q_index_in);
455 break;
456 }
458 {
459 call_reinit(*fe_values_outside_full_quadrature, q_index_in);
460 break;
461 }
463 {
464 unsigned int q_index = q_index_in;
465
466 if (q_index == numbers::invalid_unsigned_int)
467 {
468 unsigned int active_fe_index = active_fe_index_in;
469
470 if (active_fe_index == numbers::invalid_unsigned_int)
471 {
472 if constexpr (std::is_same_v<
474 CellAccessorType> ||
475 std::is_same_v<
477 CellAccessorType>)
478 active_fe_index = cell->active_fe_index();
479 else
480 active_fe_index = 0;
481 }
482
483 if (q_collection_1D.size() > 1)
484 q_index = active_fe_index;
485 else
486 q_index = 0;
487 }
488
489 AssertIndexRange(q_index, q_collection_1D.size());
490
491 face_quadrature_generator.set_1D_quadrature(q_index);
492 face_quadrature_generator.generate(cell, face_no);
493
494 const Quadrature<dim - 1> &inside_quadrature =
495 face_quadrature_generator.get_inside_quadrature();
496 const Quadrature<dim - 1> &outside_quadrature =
497 face_quadrature_generator.get_outside_quadrature();
498
499 // Even if a cell is formally intersected the number of created
500 // quadrature points can be 0. Avoid creating an FEInterfaceValues
501 // object if that is the case.
502 if (inside_quadrature.size() > 0)
503 {
504 fe_values_inside.emplace(*mapping_collection,
505 *fe_collection,
507 inside_quadrature),
508 region_update_flags.inside);
509
510 call_reinit(*fe_values_inside, /*q_index=*/0);
511 }
512
513 if (outside_quadrature.size() > 0)
514 {
515 fe_values_outside.emplace(*mapping_collection,
516 *fe_collection,
518 outside_quadrature),
519 region_update_flags.outside);
520
521 call_reinit(*fe_values_outside, /*q_index=*/0);
522 }
523 break;
524 }
525 default:
526 {
528 break;
529 }
530 }
531 }
532
533
534
535 template <int dim>
536 const std::optional<::FEInterfaceValues<dim>> &
538 {
539 if (current_face_location == LocationToLevelSet::inside)
540 return fe_values_inside_full_quadrature;
541 else
542 return fe_values_inside;
543 }
544
545
546
547 template <int dim>
548 const std::optional<::FEInterfaceValues<dim>> &
550 {
551 if (current_face_location == LocationToLevelSet::outside)
552 return fe_values_outside_full_quadrature;
553 else
554 return fe_values_outside;
555 }
556
557
558#include "fe_values.inst"
559
560} // namespace NonMatching
const std::optional<::FEInterfaceValues< dim > > & get_outside_fe_values() const
Definition fe_values.cc:549
void initialize(const hp::QCollection< dim - 1 > &q_collection)
Definition fe_values.cc:397
const hp::QCollection< 1 > q_collection_1D
Definition fe_values.h:664
void do_reinit(const TriaIterator< CellAccessorType > &cell, const unsigned int face_no, const unsigned int q_index, const unsigned int active_fe_index, const std::function< void(::FEInterfaceValues< dim > &, const unsigned int)> &call_reinit)
Definition fe_values.cc:434
const std::optional<::FEInterfaceValues< dim > > & get_inside_fe_values() const
Definition fe_values.cc:537
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:339
typename FaceQuadratureGenerator< dim >::AdditionalData AdditionalData
Definition fe_values.h:491
const std::optional<::FEValues< dim > > & get_outside_fe_values() const
Definition fe_values.cc:318
void initialize(const hp::QCollection< dim > &q_collection)
Definition fe_values.cc:101
const std::optional< FEImmersedSurfaceValues< dim > > & get_surface_fe_values() const
Definition fe_values.cc:330
void reinit_internal(const CellIteratorType &cell, const unsigned int q_index, const unsigned int mapping_index, const unsigned int fe_index)
Definition fe_values.cc:178
const hp::QCollection< 1 > q_collection_1D
Definition fe_values.h:336
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int)
Definition fe_values.cc:150
typename QuadratureGenerator< dim >::AdditionalData AdditionalData
Definition fe_values.h:146
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:45
const std::optional<::FEValues< dim > > & get_inside_fe_values() const
Definition fe_values.cc:306
unsigned int size() const
unsigned int size() const
Definition collection.h:308
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ update_default
No update.
#define DEAL_II_ASSERT_UNREACHABLE()
Definition hp.h:117
static const unsigned int invalid_unsigned_int
Definition types.h:220